$\displaystyle\int d\theta\,n{[\partial_{\theta}\varphi(\theta,q)]}^{\rm dr}$
$\displaystyle=$ $\displaystyle-\int
d\theta\,n{[\partial_{\theta}\varphi(\theta,-q)]}^{\rm dr}={1}^{\rm
dr}(q)[{\varphi}^{\rm dr}(q,q)-{\varphi}^{\rm dr}(q,-q)]={1}^{\rm dr}(q)V(q),$
$\displaystyle\int d\theta\,n(\theta){\varphi}^{\rm dr}(\theta,q)$
$\displaystyle=$ $\displaystyle({1}^{\rm dr}(q)-1),$ $\displaystyle\int
d\theta\,\theta n{\varphi(\theta,\pm q)}^{\rm dr}$ $\displaystyle=$
$\displaystyle({\theta}^{\rm dr}-\theta)(\pm q),$ $\displaystyle\int
d\theta\,\theta n{[\partial_{\theta}\varphi(\theta,\pm q)]}^{\rm dr}$
$\displaystyle=$ $\displaystyle(1-{1}^{\rm dr})+{\theta}^{\rm dr}(q)U(q).$
(SG.126)

As a consequence the even functions are given by

$\displaystyle\int d\theta\,n({g}^{\rm dr}-{g}^{\rm dr}_{0})$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}T^{2}{g}^{\rm dr}(q)}{3({\theta}^{\rm
dr})^{2}}\left(-2({1}^{\rm
dr})V(q)+V(q)-\left(\frac{1}{v_{F}}-2U(q)+f(\alpha)\right)({1}^{\rm
dr}-1)\right)+\frac{\pi^{2}T^{2}{(\partial g)}^{\rm dr}}{3({\theta}^{\rm
dr})^{2}}({1}^{\rm dr}-1),$ $\displaystyle\int d\theta\,\theta n({g}^{\rm
dr}-{g}^{\rm dr}_{0})$ $\displaystyle=$ $\displaystyle 0.$ (SG.127)

And the odd functions satisfy the similar relations

$\displaystyle\int d\theta\,n({g}^{\rm dr}-{g}^{\rm dr}_{0})$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\int d\theta\,\theta n({g}^{\rm dr}-{g}^{\rm
dr}_{0})$ $\displaystyle=$ $\displaystyle\frac{\pi^{2}T^{2}{g}^{\rm
dr}(q)}{3({\theta}^{\rm dr})^{2}}\left(({1}^{\rm dr}-1)-(2{\theta}^{\rm
dr}(q)-q)U(q)-\left(\frac{1}{v_{F}}-2U(q)+f(\alpha)\right)({\theta}^{\rm
dr}(q)-q)\right)$ $\displaystyle+\frac{\pi^{2}T^{2}{(\partial g)}^{\rm
dr}(q)}{3({\theta}^{\rm dr})^{2}}({\theta}^{\rm dr}(q)-q).$

Noting that a subscript ${}^{\prime}0$’ denotes that the zero temperature
value of a function is being considered. With the functions $f(\alpha)$
introduced as potential corrections due to additional Lagrange parameters in
the thermal energy, which are expected to take the form

$\displaystyle\frac{\epsilon^{\prime\prime}}{\epsilon^{\prime}}\sim\frac{1}{v_{F}(\theta)}-U(\theta)+f(\alpha)$
(SG.129)

#### General expression for conserved charges and currents at constant
$\mu(q)$ and low $T$

A general hydrodynamic charge field is given by $g(\theta)=\theta^{n}/n!$ as

$\displaystyle Q_{n}=\int d\theta\,{g}^{\rm dr}(\theta)n(\theta).$ (SG.130)

Such the Sommerfeld expansion (SG.118) can be applied by identifying

$\displaystyle f(\theta)$ $\displaystyle=$ $\displaystyle{g}^{\rm dr},$
$\displaystyle f^{\prime}(q)-f^{\prime}(-q)$ $\displaystyle=$ $\displaystyle
2{(\partial g)}^{\rm dr}(q)-2\left({\varphi}^{\rm
dr}(q,q)-(-1)^{n}{\varphi}^{\rm dr}(q,-q)\right){g}^{\rm dr}(q).$ (SG.131)

When the function $g(\theta)$ is antisymmetric the low-temperature corrections
to the corresponding conserved charge vanish. As a consequence consider only
the case where $g(\theta)$ is symmetric

$\displaystyle Q_{n}(T)=\int n_{0}{g}^{\rm
dr}+\frac{\pi^{2}T^{2}}{3({\theta}^{\rm dr})^{2}}\left({(\partial g)}^{\rm
dr}-{g}^{\rm dr}(q)V(q)-\left(\frac{1}{v_{F}}-2U(q)\right){g}^{\rm
dr}(q)\right).$ (SG.132)

This must be combined with the low $T$ expansion for the even dressed
functions of eq. (G.5) and simplified. The overall result for the low-
temperature even conserved charges is found to be

$\displaystyle Q_{n}(T)-Q_{n}(0)=\frac{\pi^{2}T^{2}{1}^{\rm
dr}(q)}{3(\epsilon^{\prime}(q))^{2}}\left[\left(2(U(q)-V(q))-\frac{1}{v_{F}}\right){g}^{\rm
dr}(q)+{(\partial g)}^{\rm dr}(q)\right].$ (SG.133)

In addition to the general conserved charge, the general conserved current can
be written out. Specifically

$\displaystyle\mathcal{J}_{n}=\int d\theta\,vg{1}^{\rm dr}n=\int
d\theta\,\theta{g}^{\rm dr}n.$ (SG.134)

So identify the following to insert into the approximation formula

$\displaystyle f(\theta)$ $\displaystyle=$ $\displaystyle\theta{g}^{\rm dr},$
$\displaystyle f^{\prime}(q)-f^{\prime}(-q)$ $\displaystyle=$ $\displaystyle
2{g}^{\rm dr}(q)+2q{(\partial g)}^{\rm dr}(q)-2qU(q){g}^{\rm dr}(q).$ (SG.135)

Unlike the care of the hydrodynamic charges in the case of the currents only
antisymmetric $g(\theta)=\theta^{n}/n!$ need be considered. From the general
expression we then find

$\displaystyle\mathcal{J}_{n}(T)=\int d\theta\,\theta n_{0}{g}^{\rm
dr}+\frac{\pi^{2}T^{2}}{3({\theta}^{\rm dr})^{2}}\left({g}^{\rm
dr}(q)+q{(\partial g)}^{\rm dr}(q)-qU(q){g}^{\rm
dr}(q)-\left(\frac{1}{v_{F}}-2U(q)\right)q{g}^{\rm dr}(q)\right).$ (SG.136)

This is combined with the identity for temperature corrections to the odd
dressed functions, eq. (G.5). With the final result being

$\displaystyle\mathcal{J}_{n}(T)-\mathcal{J}_{n}(0)$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}T^{2}{(\partial g)}^{\rm dr}{\theta}^{\rm
dr}(q)}{3(\epsilon^{\prime}(q))^{2}}.$ (SG.137)

These can be used to incorporate more complicated forms of the thermal energy
by restoring the $f(\alpha)$ factor to find the charge and current corrections
at constant chemical potential as

$\displaystyle Q_{n}(T)-Q_{n}(0)$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}T^{2}{1}^{\rm
dr}(q)}{3(\epsilon^{\prime}(q))^{2}}\left[\left(2(U(q)-V(q))-\frac{1}{v_{F}}-f(\alpha)\right){g}^{\rm
dr}(q)+{(\partial g)}^{\rm dr}(q)\right],$
$\displaystyle\mathcal{J}_{n}(T)-\mathcal{J}_{n}(0)$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}T^{2}{(\partial g)}^{\rm dr}{\theta}^{\rm
dr}(q)}{3(\epsilon^{\prime}(q))^{2}}.$ (SG.138)

It should also be noted that these describe the behaviour of a symmetric Fermi
sea, which can be generalized by shifting the both Fermi seas separately in
the multi-Fermi sea scenario.

#### Low-temperature pressure at constant density

Temperature corrections to the pressure $\mathcal{P}_{s}=\mathcal{J}_{1}$ can
be determined from these equations, noting that

$\displaystyle\mathcal{P}_{s}=\int_{-\infty}^{\infty}d\theta\,{\theta}^{\rm
dr}\theta n.$ (SG.139)

So that the result at low-temperature with constant chemical potential is

$\displaystyle\mathcal{P}_{s}-\mathcal{P}_{s}(0)$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}T^{2}}{3v_{F}}.$ (SG.140)

Likewise, the density at finite temperature is determined by the integral

$\displaystyle\rho=\int d\theta\,{1}^{\rm dr}(\theta)n(\theta).$ (SG.141)

By application of the even conserved charge identity the temperature
correction to density is given by

$\displaystyle\rho-\rho(0)$ $\displaystyle=$
$\displaystyle-\frac{\pi^{2}T^{2}}{3v_{F}{\theta}^{\rm
dr}}\left(2V(q)({1}^{\rm dr}(q))+\left[\frac{1}{v_{F}}-2U(q)\right]({1}^{\rm
dr}(q))\right).$ (SG.142)

With these results for the density and pressure at fixed $q$ consider the
pressure at a constant low-temperature reference state $\rho_{\rm ref}$. For
the density to be constant it must be the case that

$\displaystyle\rho_{\rm ref}$ $\displaystyle=$
$\displaystyle\rho(q(T))=\rho(q_{B}),$ (SG.143) $\displaystyle=$
$\displaystyle\rho(q(0))-T^{2}f(q(0))=\rho(q(0))+\frac{\delta\rho}{\delta
q}\bigg{|}_{q=q(0)}(q_{B}-q(0)).$

With $f(q(0))$ being the small temperature correction to density at constant
chemical potential $\mu(q)$. As a consequence the temperature correction for
the pressure in a fixed state $\rho_{\rm ref}$ is given by

$\displaystyle\mathcal{P}(\rho(q(T)))$ $\displaystyle-$
$\displaystyle\mathcal{P}(\rho(q_{B}))=\mathcal{P}(\rho(q))+\frac{\pi^{2}T^{2}}{3v_{F}}-\mathcal{P}(\rho(q))-\frac{\delta\mathcal{P}(\rho(q))}{\delta
q}\bigg{|}_{q=q(0)}(q_{B}-q(0)),$ (SG.144) $\displaystyle=$
$\displaystyle\frac{2\pi^{2}T^{2}}{3}\left(\frac{1}{v_{F}}+V(q)-U(q)\right),$
$\displaystyle=$
$\displaystyle\frac{2\pi^{2}T^{2}}{3v_{F}}\left(1-2{v_{F}}{\varphi}^{\rm
dr}(q,-q)\right)=\frac{\pi^{2}T^{2}}{3v_{F}}(1+\partial_{q}v_{F})=\frac{\pi^{2}T^{2}}{3v_{F}}(1+2K\partial_{\rho}v_{F}),$

with the note that $\partial_{\rho}P=v_{F}^{2}$ .

###  G.6 Dressed functions identities at low temperature

For considering the spatial derivative of the thermal energy it is possible to
construct a relation

$\displaystyle\frac{\partial_{x}\epsilon}{\partial_{\theta}\epsilon}=-\partial_{x}\eta-\frac{v_{F}}{v_{s}(\theta)}\partial_{x}q,$
(SG.145)

the utility of which is that the spatial derivative of a thermal sea with
arbitrary Fermi points is put into correspondence with variables relying only
on symmetric Fermi points.

At zero temperature the derivative of the dressed functions can be explicitly
written out in terms of other dressed functions. These identities are
collected here. These identities appear frequently at various intermediate
steps.

$\displaystyle\partial_{\theta}{f}^{\rm dr}\bigg{|}_{\theta=\pm q}$
$\displaystyle=$ $\displaystyle{(\partial_{\theta}f)}^{\rm dr}-{\varphi}^{\rm
dr}(\pm q,q){f}^{\rm dr}(q)+{\varphi}^{\rm dr}(\pm q,-q){f}^{\rm dr}(-q),$
$\displaystyle(\partial_{q}{f}^{\rm dr})\bigg{|}_{\theta=\pm q}$
$\displaystyle=$ $\displaystyle{\varphi}^{\rm dr}(\pm q,q){f}^{\rm
dr}(q)+{\varphi}^{\rm dr}(\pm q,-q){f}^{\rm dr}(-q),$
$\displaystyle\partial_{q}({f}^{\rm dr}(q))$ $\displaystyle=$
$\displaystyle{(\partial f)}^{\rm dr}(q)+2{\varphi}^{\rm dr}(q,-q){f}^{\rm
dr}(-q).$ (SG.146)

With these relations many identities can be written out for other composite
quantities like the effective velocity

$\displaystyle\partial_{\theta}v_{s}(\theta)\bigg{|}_{\theta=q}$
$\displaystyle=$ $\displaystyle 1-2v_{F}{\varphi}^{\rm dr}(q,-q),$
$\displaystyle\partial_{q}v_{s}(\theta)\bigg{|}_{\theta=q}$ $\displaystyle=$
$\displaystyle-2v_{F}{\varphi}^{\rm dr}(q,-q),$
$\displaystyle\partial_{q}(v_{s}(q))$ $\displaystyle=$ $\displaystyle
1-4v_{F}{\varphi}^{\rm dr}(q,-q),$ (SG.147)

as well as $v_{E}={(\theta^{2})}^{\rm dr}/(2({1}^{\rm dr}))$, where

$\displaystyle\partial_{x}(v_{E}(\theta))\bigg{|}_{\theta=q}$ $\displaystyle=$
$\displaystyle-v_{F}\partial_{x}\eta,$
$\displaystyle\partial_{\theta}v_{E}(\theta)\bigg{|}_{\theta=q}$
$\displaystyle=$ $\displaystyle v_{F}.$ (SG.148)

Finally, the dressed kernel derivative identities are given by

$\displaystyle\partial_{x}{\varphi}^{\rm dr}(q,-q)$ $\displaystyle=$
$\displaystyle 2(\partial_{x}q){\varphi}^{\rm dr}(q,q){\varphi}^{\rm
dr}(q,-q),$ $\displaystyle\partial_{q}{\varphi}^{\rm
dr}(\theta,\mu)\bigg{|}_{\theta,\mu=q,-q}$ $\displaystyle=$ $\displaystyle
2{\varphi}^{\rm dr}(q,q){\varphi}^{\rm dr}(q,-q).$ (SG.149)

Consistent with the more general identities previously noted.

###  G.7 Diffusion at large coupling

In the large $c$ limit, keeping the leading term, the diffusion can be more
directly computed. Many of the steps are similar to the interacting scenario,
however due to explicit form of the thermal energy,

$\beta\epsilon=\beta\frac{(\theta-\theta^{-})(\theta-\theta^{+})}{2}+O\left(\frac{1}{c}\right),$
(SG.150)

it is possible to write out the relations in a more explicit form. This
section is included for completeness, however the results are completely
consistent with the $c\to\infty$ limit keeping only the $O(1/c^{2})$ terms. In
the large $c$ limit ${\varphi}^{\rm dr}\sim 1/(\pi c)$ so that the leading
contribution to the diffusion kernel yields

$\displaystyle\frac{1}{2}\partial_{x}(\int
d\alpha[\tilde{D}_{\infty}]\partial_{x}n)$ $\displaystyle\overset{c\gg
1}{\sim}$ $\displaystyle\frac{1}{2\pi c^{2}}\int
d\alpha|\theta-\alpha|(\theta-\alpha)\partial_{x}\left[\frac{\partial_{\theta}n(\theta)}{\partial_{\theta}\epsilon}\frac{\partial_{\alpha}n(\alpha)}{\beta\partial_{\alpha}\epsilon}(\partial_{x}\theta^{+}+\partial_{x}\theta^{-})\right].$
(SG.151)

This may be more easily simplified by recalling that

$\displaystyle\frac{\partial_{\theta}n}{\partial_{\theta}\epsilon}=-\sum_{\sigma=\pm}\sigma\frac{\delta(\theta-\theta^{\sigma})}{\theta-\frac{(\theta^{+}+\theta^{-})}{2}}.$
(SG.152)

So find that the integral becomes

$\displaystyle\frac{1}{2}\partial_{x}\left(\int
d\alpha[\tilde{D}_{\infty}]\partial_{x}n\right)$ $\displaystyle\overset{c\gg
1}{\sim}$ $\displaystyle-\frac{1}{2\pi
c^{2}}\partial_{x}\left(\sum_{\sigma=\pm}\sigma|\theta-\theta^{\sigma}|\frac{(\theta-\theta^{\sigma})}{\theta^{\sigma}-\eta}\frac{\partial_{\theta}n(\theta)}{\beta\partial_{\theta}\epsilon}(\partial_{x}\theta^{+}+\partial_{x}\theta^{-})\right).$
(SG.153)

Note that the presence of the $\delta$-function under the $\partial_{x}$
complicates the picture, since now acting on the integral introduces a
$\delta^{\prime}$-function that cannot be matched by other terms. This
represents the fluctuations around the edge of the Fermi seas. Explicitly, the
low-temperature diffusion kernel in the large $c$ limit is given by

$\displaystyle\frac{1}{2}\partial_{x}\left(\int
d\alpha[\tilde{D}_{\infty}]\partial_{x}n\right)$ $\displaystyle=$
$\displaystyle\frac{4}{\pi
c^{2}\beta}\sum_{\sigma}\sigma\delta(\theta-\theta^{\sigma})\partial_{x}^{2}\eta+\frac{4}{\pi
c^{2}}\sum_{\sigma,\sigma^{\prime}}\left(\frac{|\theta-\theta^{\sigma^{\prime}}|(\theta-\theta^{\sigma^{\prime}})}{\beta(\theta^{+}-\theta^{-})^{2}}\right)\delta^{\prime}(\theta-\theta^{\sigma})\frac{\partial_{x}(\beta\epsilon(\theta))}{\beta\partial_{\theta}\epsilon}(\partial_{x}\eta),$
(SG.154) $\displaystyle-\frac{8}{\pi
c^{2}}\sum_{\sigma}\delta(\theta-\theta^{\sigma})\frac{(\partial_{x}\eta)^{2}}{\beta(\theta^{+}-\theta^{-})}-\frac{4q}{\pi
c^{2}\beta^{2}}(\partial_{x}\eta)(\partial_{x}\beta)\sum_{\sigma}\delta(\theta-\theta^{\sigma}).$

This an be combined with with Euler GHD in terms of a single evolution
equation for the Fermi weights

$\displaystyle\partial_{t}n+v\partial_{x}n=\frac{1}{2}\partial_{x}\left(\int
d\alpha[\tilde{D}_{\infty}]\partial_{x}n\right),$
$\displaystyle\sum_{\sigma}\sigma\delta(\theta-\theta^{\sigma})$
$\displaystyle\left(\partial_{t}\theta^{\sigma}+v^{\rm
eff}(\theta)\partial_{x}\theta^{\sigma}\right)=\sum_{\sigma}\sigma\delta(\theta-\theta^{\sigma})\left(\frac{4}{\pi
c^{2}\beta}\partial_{x}^{2}\eta-\frac{8\sigma}{\pi
c^{2}}\frac{(\partial_{x}\eta)^{2}}{\beta(\theta^{+}-\theta^{-})}-\frac{4q\sigma}{\pi
c^{2}\beta^{2}}(\partial_{x}\eta)(\partial_{x}\beta)\right)$
$\displaystyle+\frac{4}{\pi
c^{2}}\sum_{\sigma,\sigma^{\prime}}\left(\frac{|\theta-\theta^{\sigma^{\prime}}|(\theta-\theta^{\sigma^{\prime}})}{\beta(\theta^{+}-\theta^{-})^{2}}\right)\delta^{\prime}(\theta-\theta^{\sigma})\left(\frac{\partial_{x}\beta}{\beta}\frac{\epsilon(\theta)}{\partial_{\theta}\epsilon}-\partial_{x}\eta-\frac{q}{\partial_{\theta}\epsilon}\partial_{x}q\right)(\partial_{x}\eta).$
(SG.155)

In the second line the term containing $\partial_{x}\beta$ should be noted as
vanishing, since $\epsilon(\theta^{\sigma})=0$. Consequently, by inserting a
test function and simplifying it follows that the $\partial_{x}\beta$ term
from the first and second lines are either identically zero or exactly cancel.
Subsequent removal of the test function then implies

$\displaystyle\sum_{\sigma}\sigma\delta(\theta-\theta^{\sigma})$
$\displaystyle\left(\partial_{t}\theta^{\sigma}+v^{\rm
eff}(\theta)\partial_{x}\theta^{\sigma}\right)=\sum_{\sigma}\sigma\delta(\theta-\theta^{\sigma})\left(\frac{4}{\pi
c^{2}\beta}\partial_{x}^{2}\eta-\frac{8\sigma}{\pi
c^{2}}\frac{(\partial_{x}\eta)^{2}}{\beta(\theta^{+}-\theta^{-})}\right)$
$\displaystyle-\frac{4}{\pi
c^{2}}\sum_{\sigma,\sigma^{\prime}}\left(\frac{|\theta-\theta^{\sigma^{\prime}}|(\theta-\theta^{\sigma^{\prime}})}{\beta(\theta^{+}-\theta^{-})^{2}}\right)\delta^{\prime}(\theta-\theta^{\sigma})\left(\partial_{x}\eta+\frac{q}{\partial_{\theta}\epsilon}\partial_{x}q\right)(\partial_{x}\eta).$
(SG.156)

Now the $\delta^{\prime}$-function must be dealt with. By considering the
diffusive corrections to conserved charges, which allows us to make sense of
the $\delta^{\prime}$-function. When integrating over $\int d\theta\,1$ the
result is a vanishing correction. A less trivial result is obtained by
considering the integration of $\int d\theta\,\theta$, and $\int
d\theta\,\theta^{2}/2$, corresponding to the momentum and energy diffusion
respectively

$\displaystyle\frac{1}{2}\int d\theta\,\theta\partial_{x}\left(\int
d\alpha[\tilde{D}_{\infty}]\partial_{x}n\right)$ $\displaystyle=$
$\displaystyle\frac{8}{\pi
c^{2}\beta}\left(q\partial_{x}^{2}\eta+(\partial_{x}q)(\partial_{x}\eta)-\frac{\partial_{x}\beta}{\beta}q\partial_{x}\eta\right)=\frac{8}{\pi
c^{2}}\partial_{x}\left(\frac{q\partial_{x}\eta}{\beta}\right),$
$\displaystyle\frac{1}{2}\int
d\theta\,\frac{\theta^{2}}{2}\partial_{x}\left(\int
d\alpha[\tilde{D}_{\infty}]\partial_{x}n\right)$ $\displaystyle=$
$\displaystyle\frac{8\eta}{\pi
c^{2}\beta}\left(q\partial_{x}^{2}\eta+(\partial_{x}q)(\partial_{x}\eta)-\frac{\partial_{x}\beta}{\beta}q\partial_{x}\eta\right)+\frac{8q}{\pi
c^{2}\beta}\left(\partial_{x}\eta\right)^{2}.$ (SG.157)

Note that this implies the viscosity in the large $c$ limit, after
identification of $\rho=2q$, as being given by

$\displaystyle\mu(\rho,T)=\frac{4T}{\pi c^{2}}\rho,$ (SG.158)

this result is consistent with the $c\gg 1$ limit of the low-temperature
diffusion with interaction.

###  G.8 Commutativity of low-temperature limit and spatial derivative

Along with the above, the $T=0$ limit commutes with the spatial derivative.
That this commutes is demonstrated below by considering the zero temperature
Fermi weight in the large $c$ limit and an LDA of the Fermi weight with

$\displaystyle
n(\theta)=\Theta(\theta-\theta^{-})-\Theta(\theta-\theta^{+})=-\sum_{\sigma=\pm}\sigma\Theta(\theta-\theta^{\sigma}).$
(SG.159)

The spatial derivatives are then computed as

$\displaystyle\partial_{x}n(\theta)$ $\displaystyle=$
$\displaystyle\sum_{\sigma=\pm}\sigma\delta(\theta-\theta^{\sigma})\partial_{x}\theta^{\sigma},$
$\displaystyle\partial_{x}^{2}n(\theta)$ $\displaystyle=$
$\displaystyle\sum_{\sigma=\pm}\left(\sigma\delta(\theta-\theta^{\sigma})(\partial_{x}^{2}\theta^{\sigma})-\sigma\delta^{\prime}(\theta-\theta^{\sigma})(\partial_{x}\theta^{\sigma})^{2}\right).$
(SG.160)

This result will now be contrasted with the finite temperature case.
Application of the derivatives in this case makes use of the Eq. (SC.39),
$\eta=\frac{\theta^{+}+\theta^{-}}{2}$, and
$q=\frac{\theta^{+}-\theta^{-}}{2}$ to write the finite temperature second
spatial deriative as

$\displaystyle\partial^{2}_{x}n$ $\displaystyle=$
$\displaystyle-\sum_{\sigma=\pm}\sigma\delta(\theta-\theta^{\sigma})\frac{\partial_{x}^{2}\epsilon}{\partial_{\theta}\epsilon}-\sum_{\sigma=\pm}\frac{\sigma}{\partial_{\theta}\epsilon(\theta^{\sigma})}\delta^{\prime}(\theta-\theta^{\sigma})\frac{(\partial_{x}\epsilon)^{2}}{\partial_{\theta}\epsilon},$
(SG.161) $\displaystyle=$
$\displaystyle-\sum_{\sigma=\pm}\sigma\delta(\theta-\theta^{\sigma})\frac{[(\theta^{-}-\theta)\partial_{x}^{2}\theta^{+}+(\theta^{+}-\theta)\partial_{x}^{2}\theta^{-}]+2\partial_{x}\theta^{+}\partial_{x}\theta^{-}}{2(\theta-\eta)}$
$\displaystyle-\sum_{\sigma=\pm}\frac{1}{q}\delta^{\prime}(\theta-\theta^{\sigma})\frac{[(\theta-\theta^{-})\partial_{x}\theta^{+}+(\theta-\theta^{+})\partial_{x}\theta^{-}]^{2}}{4(\theta-\eta)}.$

For clarity consider the $\delta^{\prime}(\theta-\theta^{\sigma})$ term on its
own. Interpretation of the $\delta^{\prime}(\theta-\theta^{\sigma})$ motivates
the application of an overall integral over $\theta$. So inserting the
expressions for the derivatives of the energy $\epsilon$ have

$\displaystyle\int
d\theta\sum_{\sigma=\pm}\frac{1}{q}\delta^{\prime}(\theta-\theta^{\sigma})\frac{[(\theta-\theta^{-})\frac{\partial_{x}\theta^{+}}{2}+(\theta-\theta^{+})\frac{\partial_{x}\theta^{-}}{2}]^{2}}{4(\theta-\eta)}f(\theta)]$
$\displaystyle=-\left[(\partial_{x}\theta^{+}+\partial_{x}\theta^{-})\sum_{\sigma}\partial_{x}(\theta^{\sigma})\frac{f(\theta^{\sigma})}{q}-\sum_{\sigma}(\partial_{x}\theta^{\sigma})^{2}\frac{f(\theta^{\sigma})}{q}+[(\partial_{x}\theta^{+})^{2}-(\partial_{x}\theta^{-})^{2}]\partial_{\theta}f(\theta)\right],$
$\displaystyle=-\sum_{\sigma}\delta(\theta-\theta^{\sigma})\left[\sigma(\partial_{x}\theta^{\sigma})^{2}(\partial_{\theta}f)(\theta)+\frac{(\partial_{x}\theta^{+})(\partial_{x}\theta^{-})}{q}f(\theta)\right].$
(SG.162)

This is now recombined with the $\delta(\theta-\theta^{\sigma})$ term and the
equivalence between the two expressions follows.


Systems biology

∗To whom correspondence should be addressed.

# Inferring Signaling Pathways with Probabilistic Programming

Systems biology

David Merrell ${}^{\text{\sfb 1,2}}$ Anthony Gitter ${}^{\text{\sfb 1,2,3,}*}$
${}^{\text{\sf 1}}$Department of Computer Sciences, University of
Wisconsin–Madison, USA
${}^{\text{\sf 2}}$Morgridge Institute for Research, Madison, Wisconsin, USA
${}^{\text{\sf 3}}$Department of Biostatistics and Medical Informatics,
University of Wisconsin–Madison, USA.

(2020; 2020)

###### Abstract

Motivation: Cells regulate themselves via dizzyingly complex biochemical
processes called signaling pathways. These are usually depicted as a network,
where nodes represent proteins and edges indicate their influence on each
other. In order to understand diseases and therapies at the cellular level, it
is crucial to have an accurate understanding of the signaling pathways at
work. Since signaling pathways can be modified by disease, the ability to
infer signaling pathways from condition- or patient-specific data is highly
valuable.
A variety of techniques exist for inferring signaling pathways. We build on
past works that formulate signaling pathway inference as a Dynamic Bayesian
Network structure estimation problem on phosphoproteomic time course data. We
take a Bayesian approach, using Markov Chain Monte Carlo to estimate a
posterior distribution over possible Dynamic Bayesian Network structures. Our
primary contributions are (i) a novel proposal distribution that efficiently
samples sparse graphs and (ii) the relaxation of common restrictive modeling
assumptions.
Results: We implement our method, named Sparse Signaling Pathway Sampling, in
Julia using the Gen probabilistic programming language. Probabilistic
programming is a powerful methodology for building statistical models. The
resulting code is modular, extensible, and legible. The Gen language, in
particular, allows us to customize our inference procedure for biological
graphs and ensure efficient sampling.
We evaluate our algorithm on simulated data and the HPN-DREAM pathway
reconstruction challenge, comparing our performance against a variety of
baseline methods. Our results demonstrate the vast potential for probabilistic
programming, and Gen specifically, for biological network inference.
Availability: Find the full codebase at https://github.com/gitter-
lab/sspshttps://github.com/gitter-lab/ssps
Contact<EMAIL_ADDRESS>
##  1 Introduction

Signaling pathways enable cells to process information rapidly in response to
external environmental changes or intracellular cues. One of the core
signaling mechanisms is protein phosphorylation. Kinases add phosphate groups
to substrate proteins and phosphatases remove them. These changes in
phosphorylation state can act as switches, controlling proteins’ activity and
function. A protein’s phosphorylation status affects its localization,
stability, and interaction partners (Newman et al., 2014). Ultimately,
phosphorylation changes regulate important biological processes such as
transcription and cell growth, death, and differentiation (Hunter, 2009;
Kholodenko et al., 2010).

Pathway databases characterize the signaling relationships among groups of
proteins but are not tailored to individual biological contexts. Even for
well-studied pathways such as epidermal growth factor receptor-mediated
signaling, the proteins significantly phosphorylated during a biological
response can differ greatly from those in the curated pathway (Köksal et al.,
2018). The discrepancy can be due to context-specific signaling (Hill et al.,
2017), cell type-specific protein abundances, or signaling rewiring in disease
(Pawson and Warner, 2007). Therefore, there is a need to learn context-
specific signaling pathway representations from observed phosphorylation
changes. In the clinical setting, patient-specific signaling pathway
representations may eventually be able to guide therapeutic decisions (Drake
et al., 2016; Halasz et al., 2016; Eduati et al., 2020).

Diverse classes of techniques have been developed to model and infer signaling
pathways (Kholodenko et al., 2012). They take approaches including Granger
causality (Shojaie and Michailidis, 2010; Carlin et al., 2017), information
theory (Cheong et al., 2011; Krishnaswamy et al., 2014), logic models (Eker et
al., 2002; Guziolowski et al., 2013; Gjerga et al., 2020), differential
equations (Schoeberl et al., 2002; Molinelli et al., 2013; Henriques et al.,
2017), non-parametric statistical tests (Zhang and Song, 2013), and
probabilistic graphical models (Sachs et al., 2005) among others. Some
signaling pathway reconstruction algorithms take advantage of perturbations
such as receptor stimulation or kinase inhibition. Although perturbing
individual pathway members can causally link them to downstream
phosphorylation changes, characterizing a complex pathway can require a large
number of perturbation experiments. Inferring pathway structure from temporal
phosphorylation data presents an attractive alternative. A single time series
phosphorylation dataset can reveal important dynamics without perturbing
individual pathway members. For instance, a kinase cannot phosphorylate
substrates before it is activated.

An alternative approach to pathway reconstruction selects a context-specific
subnetwork from a general background network. These algorithms can use
phosphorylation data to assign scores to protein nodes in a protein-protein
interaction network. They then select edges that connect the high-scoring
nodes, generating a subnetwork that may explain how the induced
phosphorylation changes arise from the source of stimulation. Extensions
accommodate temporal scores on the nodes (Patil et al., 2013; Budak et al.,
2015; Köksal et al., 2018; Norman and Cicek, 2019).

Our present work builds on past techniques that formulate signaling pathway
inference as a Dynamic Bayesian Network (DBN) structure estimation problem.
This family of techniques relies on two core ideas: (i) we can use a DBN to
model phosphorylation time series data; and (ii) the DBN’s structure
translates directly to a directed graph representing the signaling pathway.
Rather than identifying a single DBN that best fits the data, these techniques
take a Bayesian approach—they yield a _posterior distribution_ over possible
DBN structures. Some techniques use Markov Chain Monte Carlo (MCMC) to sample
from the posterior (Werhli and Husmeier, 2007; Gregorczyk, 2010). Others use
exact, enumerative inference to compute posterior probabilities (Hill et al.,
2012; Oates et al., 2014; Spencer et al., 2015).

We present a new Bayesian DBN-based technique, Sparse Signaling Pathway
Sampling (SSPS). It improves on past MCMC methods by using a novel proposal
distribution specially tailored for the large, sparse graphs prevalent in
biological applications. Furthermore, SSPS makes weaker modeling assumptions
than other DBN approaches. As a result, SSPS scales to larger problem sizes
and yields superior predictions in comparison to other DBN techniques.

We implement SSPS using the Gen probabilistic programming language (PPL).
Probabilistic programming is a powerful methodology for building statistical
models. It enables the programmer to build models in a legible, modular,
reusable fashion. This flexibility was important for prototyping and
developing the current form of SSPS and readily supports future improvements
or extensions.

##  2 Materials and methods

###  2.1 Model formulation

SSPS makes specific modeling assumptions. We start with the DBN model of Hill
et al. (2012), relax some assumptions, and modify it in other ways to be
better-suited for MCMC inference.

#### Preliminary definitions.

We first define some notation for clarity’s sake. Let $G$ denote a _directed
graph_ with vertices $V$ and edges $E(G)$. Graph $G$ will represent a
signaling pathway, with vertices $V$ corresponding to proteins and edges
$E(G)$ indicating their influence relationships. We use $\text{pa}_{G}(i)$ to
denote the _parents_ of vertex $i$ in $G$.

Let $X$ denote our time series data, consisting of $|V|$ variables measured at
$T$ timepoints. $X$ is a $T{\times}|V|$ matrix where the $j$th column
corresponds to the $j$th variable and the $j$th graph vertex. As a convenient
shorthand, let $X_{+}$ denote the _latest_ $T{-}1$ timepoints in $X$, and let
$X_{-}$ denote the _earliest_ $T{-}1$ timepoints in $X$. Lastly, define
$B_{j}\equiv X_{-,\text{pa}_{G}(j)}$. In other words, $B_{j}$ contains the
values of variable $j$’s parents at the $T{-}1$ earliest timepoints. In
general, $B_{j}$ may also include columns of nonlinear interactions between
the parents. We will only include linear terms, unless stated otherwise.

#### Model derivation.

In our setting, we aim to infer $G$ from $X$. In particular, Bayesian
approaches seek a _posterior distribution_ $P(G|X)$ over possible graphs. From
Bayes’s rule we know $P(G|X)\propto P(X|G)\cdot P(G).$ That is, a Bayesian
model is fully specified by its choice of _prior distribution_ $P(G)$ and
_likelihood function_ $P(X|G)$.

We derive our model from the one used by Hill et al. (2012). They choose a
prior distribution of the form

$P(G~{}|~{}G^{\prime},\lambda)\propto\exp\left(-\lambda|E(G)\setminus
E(G^{\prime})|\right)$ (1)

parameterized by a _reference graph_ $G^{\prime}$ and _inverse temperature_
$\lambda$. This prior gives uniform probability to all subgraphs of
$G^{\prime}$ and ‘‘penalizes’’ edges not contained in $E(G^{\prime})$.
$\lambda$ controls the ‘‘importance’’ given to the reference graph.

Hill et al. choose a Gaussian DBN for their likelihood function. Intuitively,
they assume linear relationships between variables and their parents:

$X_{+,j}\sim\mathcal{N}(B_{j}\beta_{j},\sigma_{j}^{2})\hskip 40.00006pt\forall
j\in\\{1\ldots|V|\\}.$

A suitable prior over the regression coefficients $\beta_{j}$ and noise
parameters $\sigma_{j}^{2}$ (Figure 1) allows us to integrate them out,
yielding this _marginal likelihood function_ :

$P(X|G)\propto\prod\limits_{j=1}^{|V|}T^{-\frac{|\text{pa}_{G}(j)|}{2}}\left(X_{+,j}^{\top}X_{+,j}-\frac{T{-}1}{T}X_{+,j}^{\top}(B_{j}\hat{\beta}_{ols})\right)^{-\frac{T-1}{2}}$
(2)

where $\hat{\beta}_{ols}=(B_{j}^{\top}B_{j})^{-1}B_{j}^{\top}X_{+,j}$ is the
ordinary least squares estimate of $\beta_{j}$. For notational simplicity,
Equation 2 assumes we have a single time course of length $T$. In general,
there may be multiple time course replicates with differing lengths. The
marginal likelihood generalizes to that case in a straightforward way.

In SSPS we use the same marginal likelihood function (Equation 2), but a
different prior distribution $P(G)$. We obtain our prior distribution by
decomposing Equation 1 into a product of independent Bernoulli trials over
graph edges. This decomposition in turn allows us to make some useful
generalizations. Define _edge existence variables_
$z_{ij}\equiv\mathbbm{1}[(i,j)\in E(G)]$. Let $Z$ be the $|V|{\times}|V|$
matrix of all $z_{ij}$. Then we can rewrite Equation 1 as follows:

$P(G|G^{\prime},\lambda)~{}\equiv~{}P(Z|G^{\prime},\lambda)~{}\propto~{}\prod_{(i,j)\notin
E(G^{\prime})}\\!\\!\\!e^{-z_{ij}\lambda}$ $=\prod_{(i,j)\in
E(G^{\prime})}\\!\\!\\!\left(\frac{1}{2}\right)^{z_{ij}}\\!\left(\frac{1}{2}\right)^{1{-}z_{ij}}\\!\\!\\!\\!\\!\prod_{(i,j)\notin
E(G^{\prime})}\left(\frac{e^{-\lambda}}{1{+}e^{-\lambda}}\right)^{z_{ij}}\\!\left(\frac{1}{1{+}e^{-\lambda}}\right)^{1{-}z_{ij}}$

where the last line is a true equality—it gives a normalized probability
measure. We see that the original prior is simply a product of Bernoulli
variables parameterized by a shared inverse temperature, $\lambda$. See
Appendix A for a more detailed derivation.

Rewriting the prior in this form opens the door to generalizations. First, we
address a shortcoming in the way reference graph $G^{\prime}$ expresses prior
knowledge. The original prior assigns equal probability to every edge of
$G^{\prime}$. However, in practice we may have differing levels of prior
confidence in the edges. We address this by allowing a real-valued prior
confidence $c_{ij}$ for each edge:

$P(Z|C,\lambda)=\prod_{(i,j)}\left(\frac{e^{-\lambda}}{e^{-c_{ij}\lambda}{+}e^{-\lambda}}\right)^{z_{ij}}\\!\left(\frac{e^{-c_{ij}\lambda}}{e^{-c_{ij}\lambda}{+}e^{-\lambda}}\right)^{1{-}z_{ij}}$
(3)

where $C$ is the matrix of all prior confidences $c_{ij}$, replacing
$G^{\prime}$. Notice that if every $c_{ij}{\in}\\{0,1\\}$, then Equation 3 is
equivalent to the original prior. In effect, Equation 3 _interpolates_ the
original prior, permitting a continuum of confidences on the interval $[0,1]$.

We make one additional change to the prior by replacing the shared $\lambda$
inverse temperature variable with a collection of variables,
$\Lambda=\\{\lambda_{j}~{}|~{}j=1,{\ldots},|V|\\}$, one for each vertex of the
graph. Recall that the original $\lambda$ variable determined the importance
of the reference graph. In the new formulation, each $\lambda_{j}$ controls
the importance of the prior knowledge for vertex $j$ and its parents:

$P(Z|C,\Lambda)=\prod_{(i,j)}\left(\frac{e^{-\lambda_{j}}}{e^{-c_{ij}\lambda_{j}}{+}e^{-\lambda_{j}}}\right)^{z_{ij}}\\!\left(\frac{e^{-c_{ij}\lambda_{j}}}{e^{-c_{ij}\lambda_{j}}{+}e^{-\lambda_{j}}}\right)^{1{-}z_{ij}}$
(4)

We introduced $\Lambda$ primarily to help MCMC converge more efficiently.
Experiments with the shared $\lambda$ revealed a multimodal posterior that
tended to trap $\lambda$ in high or low values. The introduction of vertex-
specific $\lambda_{j}$ variables yielded faster convergence with weaker
modeling assumptions—an improvement in both respects.

We implicitly relax the model assumptions further via our inference procedure.
For sake of tractability, the original exact method of Hill et al. (2012)
imposes a hard constraint on the in-degree of each vertex. In contrast, we use
a MCMC inference strategy with no in-degree constraints.

In summary, our model departs from that of Hill et al. (2012) in three
important respects. It permits real-valued prior confidences $C$, introduces
vertex-specific inverse temperature variables $\Lambda$, and places no
constraints on vertices’ in-degrees. See the full model in Figure 1 and
Appendix A for additional details.

$\displaystyle\lambda_{j}$
$\displaystyle\sim\text{Uniform}(\lambda_{\text{min}},\lambda_{\text{max}})$
$\displaystyle\forall j\in\\{1\ldots|V|\\}$ $\displaystyle
z_{ij}~{}|~{}c_{ij},\lambda_{j}$
$\displaystyle\sim\text{Bernoulli}\left(\frac{e^{-\lambda_{j}}}{e^{-c_{ij}\lambda_{j}}+e^{-\lambda_{j}}}\right)$
$\displaystyle\forall i,j\in\\{1\ldots|V|\\}$ $\displaystyle\sigma_{j}^{2}$
$\displaystyle\propto\frac{1}{\sigma_{j}^{2}}$ $\displaystyle\forall
j\in\\{1\ldots|V|\\}$ $\displaystyle\beta_{j}~{}|~{}\sigma_{j}^{2}$
$\displaystyle\sim\mathcal{N}\left(0,T\sigma_{j}^{2}(B_{j}^{\top}B_{j})^{-1}\right)$
$\displaystyle\forall j\in\\{1\ldots|V|\\}$ $\displaystyle
X_{+,j}~{}|~{}B_{j},\beta_{j},\sigma_{j}^{2}$
$\displaystyle\sim\mathcal{N}\left(B_{j}\beta_{j},\sigma_{j}^{2}I\right)$
$\displaystyle\forall j\in\\{1\ldots|V|\\}$

Figure 1: Our generative model. (top) Plate notation. DBN parameters
$\beta_{j}$ and $\sigma_{j}^{2}$ have been marginalized out. (bottom) Full
probabilistic specification. We usually set $\lambda_{\text{min}}\simeq 3$ and
$\lambda_{\text{max}}{=}15$. If $\lambda_{\text{min}}{>}0$ is too small,
Markov chains will occasionally be initialized with very large numbers of
edges, causing computational issues. The method is insensitive to
$\lambda_{\text{max}}$ as long as it’s sufficiently large. Notice the improper
prior $1/\sigma_{j}^{2}$. In this specification $B_{j}$ denotes
$X_{-,\text{pa}_{Z}(j)}$; that is, the parents of vertex $j$ depend on edge
existence variables $Z$.

###  2.2 Inference procedure

Our method uses MCMC to infer posterior edge existence probabilities. As
described in Section 2.1, our model contains two classes of unobserved random
variables: (i) the edge existence variables $Z$ and (ii) the inverse
temperature variables $\Lambda$. For each step of MCMC, we loop through these
variables and update them in a Metropolis-Hastings fashion.

#### Main loop.

At a high level, our MCMC procedure consists of a loop over the graph
vertices, $V$. For each vertex $j$, we update its inverse temperature variable
$\lambda_{j}$ and then update its _parent set_ $\text{pa}_{G}(j)$. All of
these updates are Metropolis-Hastings steps; the proposal distributions are
described below. Each completion of this loop yields one iteration of the
Markov chain.

#### Proposal distributions.

For the inverse temperature variables we use a symmetric Gaussian proposal:
$\lambda_{j}^{\prime}\sim\mathcal{N}(\lambda_{j},\xi^{2})$. In practice the
method is insensitive to $\xi$; we typically set $\xi{=}3$.

The parent set proposal distribution is more complicated. There are two
principles at work when we design a graph proposal distribution: (i) the
proposal should efficiently traverse the space of directed graphs, and (ii) it
should favor graphs with higher posterior probability. The most widely used
graph proposal distribution selects a _neighboring_ graph uniformly from the
set of possible ‘‘add-edge,’’ ‘‘remove-edge,’’ and ‘‘reverse-edge’’ actions
(Werhli and Husmeier, 2007; Gregorczyk, 2010). We’ll refer to this traditional
proposal distribution as the _uniform graph proposal_. In our setting, we
expect sparse graphs to be much more probable than dense ones—notice how the
marginal likelihood function (Equation 2) strongly penalizes
$|\text{pa}_{G}(j)|$. However, the uniform graph proposal exhibits a
preference toward _dense graphs_. It proposes ‘‘add-edge’’ actions too often.
This motivates us to design a new proposal distribution tailored for sparse
graphs—one that operates on our sparse _parent set_ graph representation.

For a given graph vertex $j\in V$, the parent set proposal distribution
updates $\text{pa}_{G}(j)$ by choosing from the following actions:

  * •

add-parent. Select one of vertex $j$’s non-parents uniformly at random, and
add it to $\text{pa}_{G}(j)$.

  * •

remove-parent. Select one of vertex $j$’s parents uniformly at random, and
remove it from $\text{pa}_{G}(j)$.

  * •

swap-parent. A simultaneous application of add-parent and remove-parent.
Perhaps surprisingly, this action is not made redundant by the other two. It
plays an important role by yielding updates that maintain the size of the
parent set. Because the marginal likelihood (Equation 2) changes steeply with
$|\text{pa}_{G}(j)|$, Metropolis-Hastings acceptance probabilities will be
higher for actions that keep $|\text{pa}_{G}(j)|$ constant.

These three actions are sufficient to explore the space of directed graphs,
but we need another mechanism to bias the exploration toward _sparse_ graphs.
We introduce this preference via the _probability_ assigned to each action.
Intuitively, we craft the action probabilities so that when
$|\text{pa}_{G}(j)|$ is too small, add-parent moves are most probable. When
$|\text{pa}_{G}(j)|$ is too big, remove-parent moves are most probable. When
$|\text{pa}_{G}(j)|$ is about right, all moves are equally probable.

We formulate the action probabilities for vertex $j$ as follows. As a
shorthand, let $s_{j}=|\text{pa}_{G}(j)|$ and define the _reference size_
$\hat{s}_{j}=\sum_{i=1}^{|V|}c_{ij}$. That is, $\hat{s}_{j}$ uses the prior
edge confidences $C$ to estimate an appropriate reference size for the parent
set. Then, the action probabilities are

$\displaystyle p(\texttt{add-parent}|s_{j},\hat{s}_{j})$ $\displaystyle\propto
1-\left(\frac{s_{j}}{|V|}\right)^{\gamma(\hat{s}_{j})}$ $\displaystyle
p(\texttt{remove-parent}|s_{j},\hat{s}_{j})$
$\displaystyle\propto\left(\frac{s_{j}}{|V|}\right)^{\gamma(\hat{s}_{j})}$
$\displaystyle p(\texttt{swap-parent}|s_{j},\hat{s}_{j})$
$\displaystyle\propto
2\left(\frac{s_{j}}{|V|}\right)^{\gamma(\hat{s}_{j})}\cdot\left(1-\left(\frac{s_{j}}{|V|}\right)^{\gamma(\hat{s}_{j})}\right)$

where $\gamma(\hat{s}_{j})=1/\log_{2}(|V|/\hat{s}_{j})$. We use these
functional forms only because they have certain useful properties: (i) when
$s_{j}{=}0$, the probability of add-parent is 1; (ii) when $s_{j}{=}|V|$, the
probability of remove-parent is 1; and (iii) when $s_{j}{=}\hat{s}_{j}$, all
actions have equal probability (Figure 2). Beyond that, these probabilities
have no particular justification. We provide additional information about the
parent set proposal in Appendix B.

Figure 2: Action probabilities as a function of parent set size. The reference
size $\hat{s}$ is determined from prior knowledge. It approximates the size of
a ‘‘typical’’ parent set. When $s{<}\hat{s}$, add-parent is most probable;
when $s{>}\hat{s}$, remove-parent is most probable; and when $s{=}\hat{s}$,
all actions have equal probability.

Recall that Metropolis-Hastings requires us to compute the reverse transition
probability for any proposal we make. This could pose a challenge given our
relatively complicated parent set proposal distribution. However, Gen provides
a helpful interface for computing reverse probabilities. The user can provide
an _involution_ function that returns the reverse of a given action. Gen then
manages the reverse probabilities without further intervention. This makes it
relatively easy to implement Metropolis-Hastings updates with unusual proposal
distributions.

#### Termination, convergence, and inference.

We follow the basic MCMC protocols described by Gelman et al. (2014). This
entails running multiple (i.e., 4) Markov chains and discarding the first half
of each chain as burnin. In all of our analyses, we terminate each Markov
chain when it either (i) reaches a length of 100,000 iterations or (ii) the
execution time exceeds 12 hours. These termination conditions are arbitrary
but emulate a real-world setting where it may be acceptable to let the method
run overnight.

Upon termination, we assess convergence with two diagnostics: Potential Scale
Reduction Factor (PSRF) and effective number of samples ($N_{\text{eff}}$).
PSRF identifies cases where the Markov chains fail to mix or achieve
stationarity. $N_{\text{eff}}$ provides a sense of ‘‘sample size’’ for our
inferred quantities. It adjusts the number of MCMC samples by accounting for
autocorrelation in each chain. For our purposes, we say a quantity has _failed
to converge_ if its PSRF $\geq{1.01}$ or $N_{\text{eff}}{<}10$. Note that
satisfying these criteria does not guarantee convergence. However, failure to
satisfy them is a reliable flag for non-convergence.

Assuming a quantity hasn’t failed to converge, we estimate it by simply taking
its sample mean from all samples remaining after burnin. In our setting we are
primarily interested in _edge existence_ probabilities; i.e., we compute the
fraction of samples containing each edge.

###  2.3 Probabilistic programming implementation

We implemented SSPS using the Gen PPL. We briefly describe the probabilistic
programming methodology and its advantages in our setting.

#### Probabilistic programming.

Probabilistic programming is a methodology for building statistical models.
It’s based on the idea that statistical models are _generative processes_
—sequences of operations on random variables. In probabilistic programming, we
express the generative process as a program written in a PPL. This program is
then compiled to produce a log-probability function, which can be used in
inference tasks. Probabilistic programming systems typically provide a set of
generic inference methods for performing those tasks—e.g., MCMC or Variational
Bayes.

Compare this with a more traditional approach, where the user must (i) derive
and implement the log-probability function and (ii) implement an inference
method that operates on the log-probability function. This process of manual
derivation and implementation is error-prone and requires a high degree of
expertise from the user. In contrast, probabilistic programming only requires
the user to express their model in a PPL. The probabilistic programming system
manages other details.

Probabilistic programming also tends to promote good software engineering
principles such as abstraction, modularity, and legibility. Most PPLs organize
code into functions, which can be reused by multiple statistical models.

#### Probabilistic programming languages.

Several PPLs have emerged in recent years. Examples include Stan (Carpenter et
al., 2017), Edward2 (Dillon et al., 2017), Pyro (Bingham et al., 2018), PyMC3
(Salvatier et al., 2016), and Gen (Cusumano-Towner et al., 2019). PPLs differ
in how they balance _expressive power_ and _ease of use_. For example, Stan
makes it easy to build hierarchical statistical models with continuous
variables but caters poorly to other model classes. At the other extreme, Gen
can readily express a large class of models—discrete and continuous variables
with complex relationships—but requires the user to design a customized
inference procedure.

#### Implementation in Gen.

We use the Gen PPL precisely for its expressive power and customizable
inference. While implementing SSPS, the customizability of Gen allowed us to
begin with simple prototypes and then make successive improvements. For
example, our model initially used a dense _adjacency matrix_ representation
for $G$, but subsequent optimizations led us to use a sparse _parent set_
representation instead. Similarly, our MCMC method started with a naïve ‘‘add
or remove edge’’ proposal distribution; we arrived at our sparse proposal
distribution (Section 2.2) after multiple refinements. Other PPLs do not allow
this level of control (Table 1).

PPL |  Host language |  Primary
model class |  Primary inference method
---|---|---|---
Stan |  custom language |  hierarchical, cont’s vars |  Black-box HMC
Edward2 |  Python/ TensorFlow |  ‘‘deep’’, cont’s vars |  Black-box variational
PyMC3 |  Python/Theano |  ‘‘deep’’, cont’s vars |  Black-box HMC
Pyro |  Python/PyTorch |  ‘‘deep’’, cont’s vars |  Black-box variational
Gen |  Julia |  discrete and cont’s vars; highly flexible |  Customizable MCMC

Table 1: A coarse comparison of some noteworthy PPLs. Gen provides
expressiveness but requires the user to implement an inference program for
their model. Cont’s vars: continuous variables; HMC: Hamiltonian Monte Carlo.

###  2.4 Simulation study evaluation

We use a simulation study to answer important questions about SSPS: How does
its computational expense grow with problem size? Is it able to correctly
identify true edges? What is its sensitivity to errors in the prior knowledge?
Simulations allow us to answer these questions in a controlled setting where
we have access to ground truth.

#### Data simulation process.

We generate each simulated dataset as follows:

  1. 1.

Sample a random adjacency matrix $A\in\\{0,1\\}^{|V|{\times}|V|}$, where each
entry is the outcome of a $\text{Bernoulli}(p)$ trial. $A$ specifies the
_structure_ of a DBN. We choose $p{=}5/|V|$ so that each vertex has an average
of 5 parents. This approximates the sparsity we might see in signaling
pathways. We denote the size of the original edge set as $|E_{0}|$.

  2. 2.

Let the weights $\beta$ for this DBN be drawn from a normal distribution
$\mathcal{N}(0,1/\sqrt{|V|})$. We noticed empirically that the $1/\sqrt{|V|}$
scale prevented the simulated time series from diverging to infinity.

  3. 3.

Use the DBN defined by $A,\beta$ to simulate $M$ time courses of length $T$.
We imitate the real datasets in Section 2.5 by generating $M{=}4$ time
courses, each of length $T{=}8$.

  4. 4.

Corrupt the adjacency matrix $A$ in two steps: (i) remove $r\cdot|E_{0}|$ of
the edges from $A$; (ii) add $a\cdot|E_{0}|$ spurious edges to the adjacency
matrix. This corrupted graph simulates the _imperfect prior knowledge_
encountered in reality. The parameters $r$ and $a$ control the ‘‘false
negatives’’ and ‘‘false positives’’ in the prior knowledge, respectively.

We use a range of values for parameters $|V|,r,$ and $a$, yielding a grid of
simulations summarized in Table 2. See Appendix C and Figure 6 for additional
details.

Parameter |  Meaning |  Values
---|---|---
$|V|$ |  Number of variables |  40, 100, 200
$T$ |  Time course length |  8
$M$ |  Number of time courses |  4
$r$ |  Fraction of original edges removed |  0.1, 0.5, 0.75, 1.0
$a$ |  Fraction of spurious edges added |  0.1, 0.5, 0.75, 1.0

Table 2: These parameters define the grid of simulated datasets in our
simulation study. There are $3{\times}4{\times}4{=}48$ distinct grid points.
For each one, we generate $K{=}5$ replicates for a total of 240 simulated
datasets. The graph corruption parameters, $r$ and $a$, range from very little
error (0.1) to total corruption (1.0).

#### Performance metrics.

We are primarily interested in SSPS’s ability to correctly recover the
structure of the underlying signaling pathway. The simulation study allows us
to measure this in a setting where we have access to ground truth. We treat
this as a probabilistic binary classification task, where the method assigns
an _existence confidence_ to each possible edge. We measure classification
performance using area under the precision-recall curve (AUCPR). We use
_average precision_ to estimate AUCPR, as opposed to the trapezoidal rule
(which tends to be overly-optimistic, see Davis and Goadrich (2006); Flach and
Kull (2015)).

Our decision to use AUCPR is motivated by the sparseness of the graphs. For
sparse graphs the number of edges grows linearly with $|V|$ while the number
of possible edges grows quadratically. Hence, as $|V|$ grows, the proportion
of positive instances decreases and the classification task increasingly
becomes a ‘‘needle-in-haystack’’ scenario.

Performance measurements on simulated data come with many caveats. It’s most
instructive to think of simulated performance as a sanity check. Since our
data simulator closely follows our modeling assumptions, poor performance
would suggest serious shortcomings in our method.

###  2.5 HPN-DREAM network inference challenge evaluation

We measure SSPS’s performance on experimental data by following the evaluation
outlined by the HPN-DREAM Breast Cancer Network Inference Challenge (Hill et
al., 2016). Signaling pathways differ across contexts—e.g., cell type and
environmental conditions. The challenge is to infer these context-specific
signaling pathways from time course data.

#### Dataset.

The HPN-DREAM challenge provides phosphorylation time course data from 32
biological contexts. These contexts arise from exposing 4 cell lines (BT20,
BT549, MCF7, UACC812) to 8 stimuli. For each context, there are approximately
$M{=}4$ time courses, each about $T{=}7$ time points in length. Cell lines
have differing numbers of phosphosite measurements (i.e., differing $|V|$),
ranging from 39 (MCF7) to 46 (BT20).

#### Prior knowledge.

Participants in the original challenge were free to extract prior knowledge
from any existing data sources. As part of their analysis, the challenge
organizers combined participants’ prior graphs into a set of edge
probabilities. These _aggregate priors_ summarize the participants’ collective
knowledge. They were not available to participants in the original challenge,
but we use them in our analyses of HPN-DREAM data. We provide them to each of
the baseline methods (see Section 2.6), so the resulting performance
comparisons are fair. We do not compare any of our scores to those listed by
Hill et al. (2016) in the original challenge results.

#### Performance metrics.

The HPN-DREAM challenge aims to score methods by their ability to capture
causal relationships between pairs of variables. It estimates this by
comparing predicted _descendant sets_ against experimentally generated
descendant sets. More specifically, the challenge organizers exposed cells to
AZD8055, an mTOR inhibitor, and observed the effects on other phosphosites.
From this they determined a set of phosphosites _downstream_ of mTOR in the
signaling pathway. These include direct substrates of the mTOR kinase as well
as indirect targets.

Comparing predicted descendants of mTOR against experimentally generated
descendants of mTOR gives us a notion of _false positives_ and _false
negatives_. As we vary a threshold on edge probabilities, the predicted mTOR
descendants change, which allows us to make a receiver operating
characteristic (ROC) curve. We calculate the resulting area under the ROC
curve (AUCROC) with the trapezoidal rule to score methods’ performance on the
HPN-DREAM challenge. Hill et al. (2016) provide more details for this
descendant set AUCROC scoring metric. AUCROC is sensible for this setting
since each descendant set contains a large fraction of the variables. Sparsity
is not an issue.

Because SSPS is stochastic we run it $K{=}5$ times per context, yielding 5
AUCROC scores per context. Meanwhile the baseline methods are all
deterministic, requiring only one execution per context. We use a simple
terminology to compare SSPS’s scores against those of other methods. In a
given context, we say SSPS _dominates_ another method if its _minimum_ score
exceeds that of the other method. Conversely, we say the other method
dominates SSPS if its score exceeds SSPS’s _maximum_ score. This _dominance_
comparison has flaws—e.g., its results depend on the sample size $K$. However,
it errs on the side of strictness and suffices as an aid for summarizing the
HPN-DREAM evaluation results.

###  2.6 Baseline pathway inference algorithms

Our evaluations compare SSPS against a diverse set of baseline methods.

#### Exact DBN (Hill et al., 2012).

This method was an early inspiration for SSPS and is most similar to SSPS.
However, the exact DBN method encounters unique practical issues when we run
it on real or simulated data. The method’s computational expense increases
rapidly with problem size $|V|$ and becomes intractable unless the ‘‘max-
indegree’’ parameter is set to a small value. For example, we found that the
method used more than 32GB of RAM on problems of size $|V|{=}100$, unless max-
indegree was set ${\leq}3$. Furthermore, the exact DBN method only admits
prior knowledge in the form of Boolean _reference edges_ , rather than
continuous-valued edge confidences. We overcame this by using a threshold to
map edge confidences to 1 or 0. We chose a threshold of 0.25 for the HPN-DREAM
challenge evaluation because it yielded a reasonable number of prior edges. We
ran Hill et al.’s implementation using MATLAB 2018a.

####  FunChisq (Zhang and Song, 2013).

This method is based on the notion that two variables $X,Y$ have a causal
relationship if there exists a _functional dependence_ $Y{=}f(X)$ between
them. It detects these dependencies using a chi-square test against the ‘‘no
functional dependence’’ null hypothesis. FunChisq was a strong competitor in
the HPN-DREAM challenge, despite the fact that it uses no prior knowledge. In
order to use FunChisq, one must first discretize their time course data. We
followed Zhang and Song’s recommendation to use 1D $k$-means clustering for
discretization. Detailed instructions are given in the HPN-DREAM challenge
supplementary materials (Hill et al., 2016). We used the FunChisq (v2.4.9.1)
and Ckmeans.1d.dp (v4.3.0) R packages.

#### LASSO.

We included a variant of LASSO regression as a simple baseline. It
incorporates prior knowledge into the typical primal formulation:

$\hat{\beta_{j}}=\text{argmin}_{\beta}\left\\{\|X_{+,j}-B_{j}\beta\|_{2}^{2}~{}+~{}\alpha\sum_{i=1}^{V}e^{-c_{ij}}|\beta_{i}|\right\\}$

where $c_{ij}$ is the prior confidence (either Boolean or real-valued) for
edge $(i,j)$. That is, the method uses _penalty factors_ $e^{-c_{ij}}$ to
discourage edges with low prior confidence. The method selects LASSO
parameters, $\alpha$, using the Bayesian Information Criterion described by
Zou et al. (2007). We use GLMNet (Friedman et al., 2010) via the GLMNet.jl
Julia wrapper (v0.4.2).

#### Prior knowledge baseline.

Our most straightforward baseline simply reports the prior edge probabilities,
performing no inference at all. Ideally, a Bayesian method should do no worse
than the prior—new time course data should only _improve_ our knowledge of the
true graph. In reality, this improvement is subject to caveats about data
quality and model fit.

###  2.7 SSPS software availability

We provide the SSPS code, distributed under a MIT license, via GitHub
(https://github.com/gitter-lab/sspshttps://github.com/gitter-lab/ssps) and
archive it on Zenodo
(https://doi.org/10.5281/zenodo.3939287https://doi.org/10.5281/zenodo.3939287).
It includes a Snakemake workflow (Koster and Rahmann, 2012) for our full
evaluation pipeline, enabling the reader to reproduce our results. The code
used in this manuscript corresponds to SSPS v0.1.1.

##  3 Results

We describe evaluation results from the simulation study and HPN-DREAM network
inference challenge. SSPS competes well against the baselines, with superior
scalability to other DBN-based approaches.

###  3.1 Simulation study results

We compare our method to the baselines listed in Section 2.6. We focus
especially on the exact DBN method of Hill et al. (2012), as SSPS shares many
modeling assumptions with it.

$|V|$ |  $N/\text{cpu-hr}$ |  $N_{\text{eff}}/\text{cpu-hr}$ |  MB per chain
---|---|---|---
40 |  70000 |  400 |  500
100 |  9000 |  140 |  1200
200 |  3000 |  60 |  1000

Table 3: Computational expense of SSPS as a function of problem size $|V|$.
$N$ is the number of iterations completed by a Markov chain. $N_{\text{eff}}$
accounts for burnin and autocorrelation in the Markov chains, giving a more
accurate sense of the method’s progress. The last column gives the approximate
memory footprint of each chain. The non-monotonic memory usage is an artifact
of the chain termination conditions ($N{>}$100,000 or time ${>}12$ hours).

$|V|$ |  max indeg |  ‘‘linear’’ |  ‘‘full’’
---|---|---|---
40 |  4 |  66s |  210s
5 |  770s |  3900s
6 |  6700s |  TIMEOUT
7 |  OOM |  OOM
100 |  3 |  250s |  520s
4 |  OOM |  OOM
200 |  2 |  53s |  140s
3 |  OOM |  OOM

Table 4: Computational expense of the exact DBN method of Hill et al. (2012)
measured in CPU-seconds, as a function of problem size $|V|$ and various
parameter settings. The method imposes an in-degree constraint on each vertex,
shown in the ‘‘max indeg’’ column. The columns ‘‘linear’’ and ‘‘full’’
correspond to different _regression modes_ , i.e., which interaction terms are
included in the DBN’s conditional probability distributions. ‘‘OOM’’ (Out Of
Memory) indicates that the method exceeded a 32GB memory limit. ‘‘TIMEOUT’’
indicates that the method failed to complete within 12 hours.

#### Computational expense.

Because SSPS uses MCMC, the user may allow it to run for an arbitrary amount
of time. With this in mind, we summarize SSPS’s timing with two numbers: (i)
$N/\text{cpu-hr}$, the number of MCMC samples per CPU-hour; and (ii)
$N_{\text{eff}}/\text{cpu-hr}$, the _effective_ number of samples per CPU-
hour. We also measure the memory footprint per Markov chain, subject to our
termination conditions. We measured these numbers for each simulation in our
grid (see Table 2).

Table 3 reports average values of $N/\text{cpu-hr}$,
$N_{\text{eff}}/\text{cpu-hr}$, and memory footprint for each problem size. As
we expect, $N/\text{cpu-hr}$ and $N_{\text{eff}}/\text{cpu-hr}$ both decrease
approximately with the inverse of $|V|$. In contrast, the non-monotonic memory
usage requires more explanation. It results from two causes: (i) our
termination condition and (ii) the sparse data structures we use to store
samples. On small problems ($|V|{=}40$), the Markov chain terminates at a
length of 100,000—well within the 12-hour limit. On larger problems
($|V|{=}100$ or $200$) the Markov chain terminates at the 12-hour timeout.
This accounts for the 500MB gap between small and large problems. The
_decrease_ in memory usage between $|V|{=}100$ and $200$ results from our
sparse representations for samples. Roughly speaking, the sparse format only
stores _changes_ in the variables. So the memory consumption of a Markov chain
depends not only on $|V|$, but also on the _acceptance rate_ of the
Metropolis-Hastings proposals. The acceptance rate is smaller for $|V|{=}200$,
yielding a net decrease in memory usage.

Recall that SSPS differs from more traditional MCMC approaches by nature of
its parent set proposal distribution, which is specially designed for sparse
graphs (see Section 2.2). When we modify SSPS to instead use a naïve uniform
graph proposal, we see a striking difference in sampling efficiency. The
uniform graph proposal distribution attains $N_{\text{eff}}/\text{cpu-hr}$ of
100, 10, and 0.2 for $|V|{=}40,100,$ and $200$, respectively—drastically
smaller than those listed in Table 3 for the parent set proposal. It’s
possible that the traditional proposal could achieve higher
$N_{\text{eff}}/\text{cpu-hr}$ by simply running faster. However, the more
important consideration is how $N_{\text{eff}}/\text{cpu-hr}$ changes with
$|V|$. Our parent set proposal distribution’s $N_{\text{eff}}/\text{cpu-hr}$
decays approximately like $O(1/|V|)$. This is better than what we might expect
from a simple analysis (Appendix B). Meanwhile, the traditional proposal
distribution’s $N_{\text{eff}}/\text{cpu-hr}$ decays faster than
$O(1/|V|^{4})$. This gap between $O(1/|V|)$ and $O(1/|V|^{4})$ sampling
efficiencies makes an enormous difference on large problems.

Table 4 summarizes the computational expense of the exact DBN method (Hill et
al., 2012). The method quickly becomes impractical as the problem size grows,
unless we enforce increasingly strict in-degree restrictions. In particular,
the exact DBN method’s memory cost grows exponentially with its ‘‘max in-
degree’’ parameter. The growth becomes increasingly sharp with problem size.
When $|V|{=}200$, increasing the maximum in-degree from 2 to 3 makes the
difference between terminating in ${<}1$ minute and exceeding 32GB of memory.
Such low bounds on in-degree are unrealistic, and will likely result in poor
inference quality. In comparison, SSPS makes no constraints on in-degree, and
its memory usage scales well with problem size.

The other baseline methods—FunChisq and LASSO—are much less computationally
expensive. Both finish in seconds and require less than 100MB of memory for
each simulated task. This highlights the computationally intense nature of
Bayesian approaches. Not every scenario calls for Bayesian inference. However,
Bayesian inference is valuable in scientific settings where we’re concerned
with uncertainty quantification.

Figure 3: Heatmap of AUCPR values from the simulation study. Both DBN-based
techniques (SSPS and the exact method) score well on this, since the data is
generated by a DBN. On large problems the exact DBN method needs strict in-
degree constraints, leading to poor prediction quality. LASSO and FunChisq
both perform relatively weakly. See Figure 7 for representative ROC and PR
curves.

#### Predictive performance.

The simulation study provides a setting where we have access to ‘‘ground
truth’’—the true simulated graph. We use AUCPR to score each method’s ability
to recover the true graph’s edges.

Figure 3 shows the AUCPR scores for our grid of simulations. Each heat map
shows AUCPR as a function of graph corruption parameters, $r$ and $a$. The
heat maps are arranged by method and problem size $|V|$. Each AUCPR value is
an average over 5 replicates. SSPS maintains fairly consistent performance
across problem sizes. In contrast, the other methods’ scores decrease with
problem size. For the exact DBN method, this is partially due to the small in-
degree constraints imposed on the large problems. It is forced to trade model
accuracy for tractability.

Figure 4 reveals further insights into these results. It plots _differential_
performance with respect to the prior knowledge, in a layout analogous to
Figure 3. Specifically, it plots the $t$-statistic of each method’s AUCPR,
paired with the prior baseline’s AUCPR. Whenever the prior graph has some
informative edges, SSPS outperforms the prior. On the other hand, SSPS’s
performance deteriorates whenever the prior contains _no_ true edges (i.e.,
$r{=}1$). Under those circumstances FunChisq may be a better choice. Since it
doesn’t rely on prior knowledge at all, it outperforms the other methods when
the prior is totally corrupted. However, we expect that in most realistic
settings there exists partially-accurate prior knowledge, in which case we
expect FunChisq to perform worse than SSPS.

These results confirm SSPS’s ability to identify the true network, given
partially-accurate prior knowledge and time series data consistent with the
modeling assumptions. SSPS is fairly robust with respect to the prior’s
quality and has consistent performance across different problem sizes.

Figure 4: Heatmap of differential performance against the prior knowledge,
measured by AUCPR paired $t$-statistics. SSPS consistently outperforms the
prior knowledge across problem sizes and shows robustness to errors in the
prior knowledge.

###  3.2 HPN-DREAM challenge results

We evaluated SSPS on experimental data from the HPN-DREAM challenge. The
challenge includes time series phosphorylation data from 32 biological
contexts: 8 stimuli applied to 4 breast cancer cell lines. Methods are scored
on their ability to correctly identify the experimentally derived descendants
of mTOR. Figure 5 shows bar charts comparing the methods’ AUCROC scores in
each context. Appendix D provides additional details.

Figure 5: Methods’ performances across contexts in the HPN-DREAM Challenge.
MCMC is stochastic, so we run SSPS 5 times; the error bars show the range of
AUCROC scores. The other methods are all deterministic and require no error
bars. See Figure 8 for example predicted networks, Figure 9 for AUCPR scores,
and Figure 10 for representative ROC and PR curves.

SSPS performs satisfactorily on this task overall. Employing terminology from
Section 2.5, SSPS dominates the exact DBN method in 18 of the 32 contexts,
whereas the exact DBN method dominates SSPS in only 9 contexts. Meanwhile,
SSPS dominates FunChisq in 11 contexts and is dominated by FunChisq in 15.
This is not surprising because FunChisq was among the top competitors in the
original challenge. LASSO, on the other hand, performs poorly. SSPS dominates
LASSO in 18 contexts and is dominated in only 6.

More puzzling is the strong performance of the prior knowledge baseline. SSPS
dominates the aggregate prior in only 9 contexts and is dominated in 21. This
is not isolated to our method. FunChisq outperforms and is outperformed by the
prior knowledge in 11 and 21 contexts, respectively. The aggregate prior’s
strong performance is consistent with the results from the original HPN-DREAM
challenge; this prior outperformed all individual challenge submissions (Hill
et al., 2016). Even though the aggregate prior gives identical predictions for
each context and totally ignores the time course data, it still attains better
performance than the other methods. This suggests either (i) the data is
relatively uninformative or (ii) the evaluation metric based on mTOR’s
descendants isn’t sufficiently precise to measure context-specific
performance. We suspect the latter, because FunChisq uses no prior knowledge
but was the top performer in the HPN-DREAM challenge’s _in silico_ tasks. An
evaluation based on one node’s descendants is not as discriminative as an
evaluation of the directed edges. Many different directed graphs can have
equivalent or similar mTOR descendants. However, it is experimentally
impractical to generate the context-specific gold standard networks that would
be required for a more precise edge-based evaluation.

##  4 Discussion

We presented SSPS, a signaling pathway reconstruction technique based on DBN
structure estimation. It uses MCMC to estimate the posterior probabilities of
directed edges, employing a parent set proposal distribution specially
designed for sparse graphs. SSPS is a Bayesian approach. It takes advantage of
prior knowledge with edge-specific confidence scores and can provide
uncertainty estimates on the predicted pathway relationships, which are
valuable for prioritizing experimental validation.

SSPS scales to large problems more efficiently than past DBN-based techniques.
On simulated data, SSPS yields superior edge predictions with robustness to
flaws in the prior knowledge. Our HPN-DREAM evaluation shows SSPS performs
comparably to established techniques on a community standard task. It is
difficult to make stronger statements in the HPN-DREAM setting because the
prior knowledge baseline performs so well and we can only evaluate the
predicted mTOR descendants, not the entire pathway. However, SSPS’s
scalability among Bayesian methods, strong results in the simulation, and
competitive performance in the HPN-DREAM challenge make it an attractive
option for further investigation of real phosphorylation datasets.

There are several potential limitations of SSPS relative to alternative
pathway signaling models. Prior knowledge is not available in some organisms
or biological conditions, reducing one advantage of our Bayesian approach.
Although SSPS is more scalable than related DBN techniques, it would struggle
to scale to proteome-wide phosphoproteomic data measuring thousands of
phosphosites. For large datasets, we recommend running SSPS on a pruned
version that includes only the highest intensity or most variable
phosphosites. SSPS, like most DBN techniques, models only observed variables.
It will erroneously exclude important pathway members, such as scaffold
proteins, that are not phosphorylated. Latent variable models or background
network-based algorithms are better suited for including unphosphorylated
proteins in the pathway. Background network methods can also impose global
constraints on the predicted pathway structure, such as controlling the number
of connected components or proteins’ reachability from relevant receptors
(Köksal et al., 2018).

There are many possible ways to improve SSPS. For example, it could be
extended to jointly model related pathways in a hierarchical fashion, similar
to Oates et al. (2014) and Hill et al. (2017). Alternatively, SSPS could be
modified to accommodate causal assumptions via Pearl’s intervention operators;
see the model of Spencer et al. (2015) for a relevant example. Combining
temporal and interventional data (Cardner et al., 2019) is another rich area
for future work. On the algorithmic side, we could improve our MCMC procedure
by adaptively tuning the parameters of its proposal distributions, as
described by Gelman et al. (2014). Because SSPS is a probabilistic program, it
is naturally extensible.

## Acknowledgements

We thank UW-Madison’s Biomedical Computing Group for generously providing
compute resources; the teams developing Snakemake and HTCondor (Thain et al.,
2005) for empowering us to use those resources effectively; the Gen team
(Cusumano-Towner et al., 2019) for designing a uniquely powerful probabilistic
programming language; and the HPN-DREAM challenge organizers for providing
experimental data for our evaluation.

## Funding

This work was funded by the National Institutes of Health (award T32LM012413)
and National Science Foundation (award DBI 1553206).

## References

  * Bingham et al. (2018) Bingham, E., Chen, J. P., Jankowiak, M., Obermeyer, F., Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P., Horsfall, P., and Goodman, N. D. (2018).  Pyro: Deep Universal Probabilistic Programming.  Journal of Machine Learning Research.
  * Budak et al. (2015) Budak, G., Ozsoy, O. E., Son, Y. A., Can, T., and Tuncbag, N. (2015).  Reconstruction of the temporal signaling network in Salmonella-infected human cells.  Frontiers in Microbiology, 6, 730.
  * Cardner et al. (2019) Cardner, M., Meyer-Schaller, N., Christofori, G., and Beerenwinkel, N. (2019).  Inferring signalling dynamics by integrating interventional with observational data.  Bioinformatics, 35(14), i577–i585.
  * Carlin et al. (2017) Carlin, D. E., Paull, E. O., Graim, K., Wong, C. K., Bivol, A., Ryabinin, P., Ellrott, K., Sokolov, A., and Stuart, J. M. (2017).  Prophetic granger causality to infer gene regulatory networks.  PLOS ONE, 12(12), e0170340.
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of Statistical Software, 76(1).
  * Cheong et al. (2011) Cheong, R., Rhee, A., Wang, C. J., Nemenman, I., and Levchenko, A. (2011).  Information Transduction Capacity of Noisy Biochemical Signaling Networks.  Science, 334(6054), 354–358.
  * Cusumano-Towner et al. (2019) Cusumano-Towner, M. F., Saad, F. A., Lew, A. K., and Mansinghka, V. K. (2019).  Gen: A general-purpose probabilistic programming system with programmable inference.  In Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI 2019, New York, NY, USA. ACM.
  * Davis and Goadrich (2006) Davis, J. and Goadrich, M. (2006).  The relationship between precision-recall and ROC curves.  In Proceedings of the 23rd International Conference on Machine Learning.
  * Dillon et al. (2017) Dillon, J. V., Langmore, I., Tran, D., Brevdo, E., Vasudevan, S., Moore, D., Patton, B., Alemi, A., Hoffman, M., and Saurous, R. A. (2017).  Tensorflow distributions.  arXiv, page arXiv:1711.10604.
  * Drake et al. (2016) Drake, J. M., Paull, E. O., Graham, N. A., Lee, J. K., Smith, B. A., Titz, B., Stoyanova, T., Faltermeier, C. M., Uzunangelov, V., Carlin, D. E., Fleming, D. T., Wong, C. K., Newton, Y., Sudha, S., Vashisht, A. A., Huang, J., Wohlschlegel, J. A., Graeber, T. G., Witte, O. N., and Stuart, J. M. (2016).  Phosphoproteome Integration Reveals Patient-Specific Networks in Prostate Cancer.  Cell, 166(4), 1041–1054.
  * Eduati et al. (2020) Eduati, F., Jaaks, P., Wappler, J., Cramer, T., Merten, C. A., Garnett, M. J., and Saez-Rodriguez, J. (2020).  Patient-specific logic models of signaling pathways from screenings on cancer biopsies to prioritize personalized combination therapies.  Molecular Systems Biology, 16(2), e8664.
  * Eker et al. (2002) Eker, S., Knapp, M., Laderoute, K., Lincoln, P., Meseguer, J., and Sonmez, K. (2002).  Pathway logic: symbolic analysis of biological signaling.  Pacific Symposium on Biocomputing, pages 400–412.
  * Flach and Kull (2015) Flach, P. A. and Kull, M. (2015).  Precision-recall-gain curves: PR analysis done right.  In Advances in Neural Information Processing Systems.
  * Friedman et al. (2010) Friedman, J., Hastie, T., and Tibshirani, R. (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33.
  * Gelman et al. (2014) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2014).  Bayesian Data Analysis.  CRC Press, 3 edition.
  * Gjerga et al. (2020) Gjerga, E., Trairatphisan, P., Gabor, A., Koch, H., Chevalier, C., Ceccarelli, F., Dugourd, A., Mitsos, A., and Saez-Rodriguez, J. (2020).  Converting networks to predictive logic models from perturbation signalling data with CellNOpt.  bioRxiv, page 2020.03.04.976852.
  * Gregorczyk (2010) Gregorczyk, M. (2010).  An introduction to Gaussian Bayesian networks.  In Q. Yan, editor, Systems Biology in Drug Discovery and Development: Methods and Protocols, chapter 6, pages 121–147. Springer.
  * Guziolowski et al. (2013) Guziolowski, C., Videla, S., Eduati, F., Thiele, S., Cokelaer, T., Siegel, A., and Saez-Rodriguez, J. (2013).  Exhaustively characterizing feasible logic models of a signaling network using Answer Set Programming.  Bioinformatics, 29(18), 2320–2326.
  * Halasz et al. (2016) Halasz, M., Kholodenko, B. N., Kolch, W., and Santra, T. (2016).  Integrating network reconstruction with mechanistic modeling to predict cancer therapies.  Science Signaling, 9(455), ra114.
  * Henriques et al. (2017) Henriques, D., Villaverde, A. F., Rocha, M., Saez-Rodriguez, J., and Banga, J. R. (2017).  Data-driven reverse engineering of signaling pathways using ensembles of dynamic models.  PLOS Computational Biology, 13(2), e1005379.
  * Hill et al. (2012) Hill, S. M., Lu, Y., Molina, J., Heiser, L. M., Spellman, P. T., Speed, T. P., Gray, J. W., Mills, G. B., and Mukherjee, S. (2012).  Bayesian inference of signaling network topology in a cancer cell line.  Bioinformatics, 28, 2804–2810.
  * Hill et al. (2016) Hill, S. M., Heiser, L. M., Cokelaer, T., Unger, M., Nesser, N. K., Carlin, D. E., Zhang, Y., Sokolov, A., Paull, E. O., Wong, C. K., Graim, K., Bivol, A., Wang, H., Zhu, F., Afsari, B., Danilova, L. V., Favorov, A. V., Lee, W. S., Taylor, D., Hu, C. W., Long, B. L., Noren, D. P., Bisberg, A. J., HPN-DREAM Consortium, Gray, J. W., Kellen, M., Norman, T., Friend, S., Qutub, A. A., Fertig, E. J., Guan, Y., Song, M., Stuart, J. M., Spellman, P. T., Koeppl, H., Stovolitzky, G., Saez-Rodriguez, J., and Mukherjee, S. (2016).  Inferring causal molecular networks: empirical assessment through a community-based effort.  Nature Methods.
  * Hill et al. (2017) Hill, S. M., Nesser, N. K., Johnson-Camacho, K., Jeffress, M., Johnson, A., Boniface, C., Spencer, S. E., Lu, Y., Heiser, L. M., Lawrence, Y., Pande, N. T., Korkola, J. E., Gray, J. W., Mills, G. B., Mukherjee, S., and Spellman, P. T. (2017).  Context specificity in causal signaling networks revealed by phosphoprotein profiling.  Cell Systems, pages 73–83.
  * Hunter (2009) Hunter, T. (2009).  Tyrosine phosphorylation: thirty years and counting.  Current Opinion in Cell Biology, 21(2), 140–146.
  * Kholodenko et al. (2012) Kholodenko, B., Yaffe, M. B., and Kolch, W. (2012).  Computational Approaches for Analyzing Information Flow in Biological Networks.  Science Signaling, 5(220), re1.
  * Kholodenko et al. (2010) Kholodenko, B. N., Hancock, J. F., and Kolch, W. (2010).  Signalling ballet in space and time.  Nature Reviews Molecular Cell Biology, 11(6), 414–426.
  * Köksal et al. (2018) Köksal, A. S., Beck, K., Cronin, D. R., McKenna, A., Camp, N. D., Srivastava, S., MacGilvray, M. E., Bodik, R., Wolf-Yadlin, A., Fraenkel, E., Fisher, J., and Gitter, A. (2018).  Synthesizing signaling pathways from temporal phosphoproteomic data.  Cell Reports, 24, 3607–3618.
  * Koster and Rahmann (2012) Koster, J. and Rahmann, S. (2012).  Snakemake–a scalable bioinformatics workflow engine.  Bioinformatics, 28(19), 2520–2522.
  * Krishnaswamy et al. (2014) Krishnaswamy, S., Spitzer, M. H., Mingueneau, M., Bendall, S. C., Litvin, O., Stone, E., Pe’er, D., and Nolan, G. P. (2014).  Conditional density-based analysis of T cell signaling in single-cell data.  Science, 346(6213), 1250689.
  * Molinelli et al. (2013) Molinelli, E. J., Korkut, A., Wang, W., Miller, M. L., Gauthier, N. P., Jing, X., Kaushik, P., He, Q., Mills, G., Solit, D. B., Pratilas, C. A., Weigt, M., Braunstein, A., Pagnani, A., Zecchina, R., and Sander, C. (2013).  Perturbation biology: inferring signaling networks in cellular systems.  PLOS Computational Biology, 9.
  * Newman et al. (2014) Newman, R. H., Zhang, J., and Zhu, H. (2014).  Toward a systems-level view of dynamic phosphorylation networks.  Frontiers in Genetics, 5, 263.
  * Norman and Cicek (2019) Norman, U. and Cicek, A. E. (2019).  ST-Steiner: a spatio-temporal gene discovery algorithm.  Bioinformatics, 35(18), 3433–3440.
  * Oates et al. (2014) Oates, C. J., Korkola, J., Gray, J. W., and Mukherjee, S. (2014).  Joint estimation of multiple related biological networks.  The Annals of Applied Statistics, 8, 1892–1919.
  * Patil et al. (2013) Patil, A., Kumagai, Y., Liang, K.-c., Suzuki, Y., and Nakai, K. (2013).  Linking Transcriptional Changes over Time in Stimulated Dendritic Cells to Identify Gene Networks Activated during the Innate Immune Response.  PLOS Computational Biology, 9(11), e1003323.
  * Pawson and Warner (2007) Pawson, T. and Warner, N. (2007).  Oncogenic re-wiring of cellular signaling pathways.  Oncogene, 26(9), 1268–1275.
  * Sachs et al. (2005) Sachs, K., Perez, O., Pe’er, D., Lauffenburger, D. A., and Nolan, G. P. (2005).  Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data.  Science, 308(5721), 523–529.
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., and Fonnesbeck, C. (2016).  Probabilistic programming in Python using PyMC3.  PeerJ Computer Science, 2, e55.
  * Schoeberl et al. (2002) Schoeberl, B., Eichler-Jonsson, C., Gilles, E. D., and Müller, G. (2002).  Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors.  Nature Biotechnology, 20(4), 370–375.
  * Shojaie and Michailidis (2010) Shojaie, A. and Michailidis, G. (2010).  Discovering graphical Granger causality using the truncating lasso penalty.  Bioinformatics, 26(18), i517–i523.
  * Spencer et al. (2015) Spencer, S. E., Hill, S. M., and Mukherjee, S. (2015).  Inferring network structure from interventional time-course experiments.  The Annals of Applied Statistics, 9, 507–524.
  * Thain et al. (2005) Thain, D., Tannenbaum, T., and Livny, M. (2005).  Distributed computing in practice: the Condor experience.  Concurrency - Practice and Experience, 17(2-4), 323–356.
  * Werhli and Husmeier (2007) Werhli, A. V. and Husmeier, D. (2007).  Reconstructing gene regulatory networks with Bayesian networks by combining expression data with multiple sources of prior knowledge.  Statistical Applications in Genetics and Molecular Biology, 6.
  * Zhang and Song (2013) Zhang, Y. and Song, M. (2013).  Deciphering interactions in causal networks without parametric assumptions.  arXiv, page arXiv:1311.2707.
  * Zou et al. (2007) Zou, H., Hastie, T., and Tibshirani, R. (2007).  On the ‘‘degrees of freedom’’ of the lasso.  The Annals of Statistics, 35, 2173–2192.

## Appendix

###  A Model formulation details

We provide additional information about our graph prior and marginal
likelihood function. We also describe some implications of SSPS’s model
assumptions.

#### Derivation of graph prior (Equation 4).

We step through a more detailed derivation of SSPS’s new graph prior. We begin
with the original graph prior (Equation 1) and rewrite it in terms of the edge
existence variables $Z$:

$\displaystyle P(G|G^{\prime},\lambda)$
$\displaystyle\propto\exp\left(-\lambda|E(G)\setminus E(G^{\prime})|\right)$
$\displaystyle=\exp\left(-\lambda\sum_{(i,j)\notin
E(G^{\prime})}z_{ij}\right)$ (5) $\displaystyle=\prod_{(i,j)\notin
E(G^{\prime})}\\!\\!\\!e^{-\lambda z_{ij}}$ $\displaystyle=\prod_{(i,j)\notin
E(G^{\prime})}\\!\\!\\!\left(e^{-\lambda}\right)^{z_{ij}}$
$\displaystyle\propto\left(\frac{1}{1+e^{-\lambda}}\right)^{V^{2}-|E(G^{\prime})|}\cdot\prod_{(i,j)\notin
E(G^{\prime})}\\!\\!\\!\left(e^{-\lambda}\right)^{z_{ij}}$
$\displaystyle=\prod_{(i,j)\notin
E(G^{\prime})}\\!\\!\\!\left(\frac{1}{1+e^{-\lambda}}\right)\left(e^{-\lambda}\right)^{z_{ij}}$
$\displaystyle=\prod_{(i,j)\notin
E(G^{\prime})}\\!\\!\\!\left(\frac{1}{1+e^{-\lambda}}\right)^{1-z_{ij}}\left(\frac{e^{-\lambda}}{1+e^{-\lambda}}\right)^{z_{ij}}$
(6)

Equation 6 shows the original prior is in fact a product of independent
Bernoulli variables—the edge existence variables $z_{ij}$. Equation 6
explicitly assigns probability to the edges _not_ contained in
$E(G^{\prime})$. However, it also implicitly assigns uniform probability to
every edge _contained_ in $E(G^{\prime})$. We deduce that they are
$\text{Bernoulli}(0.5)$ variables, allowing us to write the prior
$P(Z~{}|~{}G^{\prime},\lambda)$ in the following form:

$\prod_{(i,j)\in
E(G^{\prime})}\\!\\!\\!\left(\frac{1}{2}\right)^{z_{ij}}\\!\left(\frac{1}{2}\right)^{1{-}z_{ij}}\\!\\!\\!\\!\\!\prod_{(i,j)\notin
E(G^{\prime})}\left(\frac{e^{-\lambda}}{1{+}e^{-\lambda}}\right)^{z_{ij}}\\!\left(\frac{1}{1{+}e^{-\lambda}}\right)^{1{-}z_{ij}}$
(7)

just as shown in the main text.

Now we modify the prior to use continuous-valued edge confidences $c_{ij}$
instead of Boolean reference edges $E(G^{\prime})$. Intuitively, we want to
restate Equation 7 as a single product over all $Z$ variables, rather than two
separate products. Our goal is to find a function $q(c_{ij})$ such that

$P(Z~{}|~{}C,\lambda)=\prod_{(i,j)}q(c_{ij})^{z_{ij}}(1-q(c_{ij}))^{1-z_{ij}}.$

However, in order to remain consistent with the original prior $q(c_{ij})$
ought to be monotone-increasing and satisfy these criteria:

$q(0)=e^{-\lambda}/(1+e^{-\lambda})\hskip 50.00008pt\text{and}\hskip
50.00008ptq(1)=1/2.$

It turns out that choosing

$q(c_{ij})=\frac{e^{-\lambda}}{e^{-c_{ij}\lambda}+e^{-\lambda}}$

satisfies these requirements. This brings us to Equation 3 of the main text.

From there, it is straightforward to replace the single shared $\lambda$
variable with a set of vertex-specific $\Lambda$ variables and arrive at
Equation 4.

#### Marginal likelihood function details.

Equation 2 is obtained by (i) using a Gaussian DBN as the likelihood function
for $G$, (ii) assuming certain prior distributions for the DBN parameters, and
(iii) integrating the DBN parameters out. Specifically, let $\beta_{j}$ and
$\sigma_{j}^{2}~{}~{}\forall j\in\\{1\ldots|V|\\}$ be the DBN’s weight and
noise parameters, respectively. We assume an improper prior
$\sigma_{j}^{2}\propto 1/\sigma_{j}^{2}$ for the noise and a Gaussian prior
for the weights:

$\beta_{j}|\sigma_{j}^{2}\sim\mathcal{N}\left(0,T\sigma_{j}^{2}(B_{j}^{\top}B_{j})^{-1}\right).$

In other words, SSPS uses an improper joint prior
$P(\beta_{j},\sigma_{j}^{2})=P(\beta_{j}|\sigma_{j}^{2})P(\sigma_{j}^{2})$
with $P(\sigma_{j}^{2}){\propto}1/\sigma_{j}^{2}$. This choice allows
$\beta_{j}$ and $\sigma_{j}^{2}$ to be marginalized, yielding Equation 2.

The power $-|\text{pa}_{G}(j)|/2$ in Equation 2 is correct when the DBN only
uses linear terms. Recall that $B_{j}$ may in general contain columns of
nonlinear interactions between parent variables. When that is true, the
quantity $|\text{pa}_{G}(j)|$ should be replaced by the _width of_ $B_{j}$. We
elide this detail in the main text for brevity. Our implementation uses the
correct exponent.

Our implementation of the marginal likelihood function employs least recently
used caching to reduce redundant computation. Code profiling shows that this
yields a substantial improvement to efficiency. For additional in-depth
discussion of Equation 2, we recommend the supplementary materials of Hill et
al. (2012).

#### Additional insights about SSPS’s model assumptions.

SSPS’s model has interesting properties that could lead to method
improvements. For example, when we replace the shared $\lambda$ variable with
vertex-specific $\Lambda$ variables, the model effectively becomes a set of
$|V|$ independent models. The plate notation in Figure 1 makes this clear;
$X_{-}$ is the only shared variable, and it’s fully observed. This has
algorithmic implications. For example, future versions of SSPS could
parallelize inference at the vertex level, allocating more resources to the
parent sets that converge slowly.

In the course of deriving Equation 6, we showed that our prior is a log-linear
model over edge features. Equation 5 shows this most clearly. Future versions
of SSPS could use the expressiveness of log-linear densities over higher-order
graph features to capture richer forms of prior knowledge.

###  B Parent set proposal details

A key component of SSPS is its novel _parent set proposal distribution_. We
motivate its design and discuss its computational complexity in greater
detail.

#### Parent sets instead of edges.

The marginal likelihood (Equation 2) is a function of the graph $G$. However,
it depends on $G$ only via its _parent sets_ , which are encoded in the
matrices $B_{j}$. Accordingly, SSPS represents $G$ by storing a list of
parents for each vertex.

It makes sense to use a proposal distribution that operates directly on SSPS’s
internal parent set representation. This motivates our choice of the add-
parent, remove-parent, and swap-parent proposals listed in Section 2.2. There
is a natural correspondence between (i) likelihood function, (ii) data
structure, and (iii) proposal distribution.

#### Sampling efficiency.

We provide some intuition for the parent set proposal’s superior sampling
efficiency. Let $z_{ij}$ be a particular edge existence variable. The estimate
for $z_{ij}$ converges quickly if MCMC updates $z_{ij}$ frequently. Hence, as
a proxy for sampling efficiency, consider the number of times $z_{ij}$ gets
updated per unit time. We decompose this quantity into three factors:

$\frac{z_{ij}\text{ updates}}{\text{unit time}}=\epsilon\cdot\tau\cdot\alpha$

where

$\epsilon=\frac{\text{graph proposals}}{\text{unit time}}\hskip
40.00006pt\tau=\frac{z_{ij}\text{ proposals}}{\text{graph proposal}}$
$\alpha=z_{ij}\text{ acceptance probability}$

In other words, $\epsilon$ is the time efficiency of the proposal
distribution. The factor $\tau$ is the probability that a given proposal
_touches_ $z_{ij}$. Lastly, $\alpha$ is the proposal’s Metropolis-Hastings
acceptance probability.

For a given proposal distribution, we’re interested in how these factors
depend on $|V|$. For simplicity of analysis, assume the Markov chain is in a
typical state where the graph is sparse: $|E(G)|=O(|V|)$.

For the parent set proposal, execution time has no dependence on $|V|$ and
hence $\epsilon~{}{=}~{}O(1)$. Recall that the parent set proposal resides in
an outer loop, which iterates through all $|V|$ vertices. It follows that for
any particular proposal there is a $1/|V|$ chance that it acts on vertex $j$.
After choosing vertex $j$, there is on average a $O(1/|V|)$ chance that the
proposal affects $z_{ij}$. This follows from the sparsity of the graph: vertex
$i$ is typically a non-parent of $j$ and the probability of choosing it via an
add-parent or swap-parent action is $O(1/|V|)$. Hence, the parent set proposal
has a probability $\tau{=}O(1/|V|^{2})$ of choosing $z_{ij}$. Lastly, the
acceptance probability $\alpha$ has no dependence on $|V|$ and therefore
$\alpha=O(1)$. The product of these factors gives an overall sampling
efficiency of $O(1/|V|^{2})$ for the parent set proposal.

For the uniform graph proposal, $\epsilon$’s complexity depends on the
particular implementation. For sake of generosity we assume an efficient
implementation with $\epsilon~{}{=}~{}O(1)$. The proposal chooses uniformly
from $O(|V|^{2})$ actions: add-, remove-, or reverse-edge. The probability of
choosing one that affects $z_{ij}$ is $\tau=O(1/|V|^{2})$. Recall that the
marginal likelihood decreases steeply with parent set size. It follows that
add-edge actions will typically have low acceptance probability. Since the
graph is sparse, add-edge actions are overwhelmingly probable; the probability
of _not_ landing on one is $O(1/|V|^{2})$. If we assume the acceptance
probability is high for remove-edge and reverse-edge actions, (i.e., they are
accepted whenever they’re proposed), then this suggests $\alpha=O(1/|V|^{2})$,
averaged over many proposals. The product of these factors suggests a sampling
efficiency that decays like $O(1/|V|^{4})$.

This gap between $O(1/|V|^{2})$ and $O(1/|V|^{4})$ sampling efficiencies
explains most of the difference that we saw in Section 3.1. A more detailed
analysis may reveal why the parent set proposal attains sampling efficiencies
closer to $O(1/|V|)$ in practice.

###  C Simulation study details

Figure 6: A schematic of the simulation study. We randomly generate a true
network (upper left) and use it to simulate a time series dataset (upper
right). We corrupt the true network by adding and removing edges (lower left);
solid red edges have been added, dashed red edges have been removed, and black
edges are original. This corrupted network serves as partially inaccurate
prior knowledge for the inference techniques. Each technique produces a
predicted network (lower right) by assigning a score to each possible edge.
The predicted network is evaluated with respect to the true network.

ROC Curves

PR Curves

Figure 7: Representative ROC curves (top) and PR curves (bottom) from the
simulation study. We show curves for three different simulations:
$|V|=40,100,$ and $200$ (left, middle, right respectively). Each of these
simulations used corruption parameters $r=a=0.5$.

We give additional details about the simulation study’s methodology and
results.

#### Simulation process.

The simulation process described in Section 2.4 differs from SSPS’s modeling
assumptions in several ways. Recall that the simulator constructs a DBN to
generate time series data. This simulated DBN employs nonlinear interaction
terms. The simulator assumes that the data at each timepoint is a _cubic_
function of the data at the previous timestep. In contrast, all of our
analyses ran SSPS with an assumption of _linear_ dependencies. In other words,
the data contained complexities that SSPS was unable to capture. SSPS’s
performance in the simulation study suggests that it has some robustness to
modeling assumption mismatches.

We provide an illustration of the simulation process in Figure 6. It is
interesting to notice that the simulated networks do not resemble directed
acyclic graphs (DAGs) in any way. They do not have any sense of
directionality. Contrast this with the biological graphs shown in Figure 8.
Strictly speaking these are not DAGs, but they do have an overall direction.
Some vertices are source-like, and others are sink-like. Future simulations
and models could be more biologically realistic if they incorporated this kind
of structure.

#### Simulation study results.

Figure 7 gives some representative ROC and PR curves from the simulation
study. On problem sizes up to $|V|=100$, SSPS and the exact DBN method yield
similar curves in both ROC and PR space—though SSPS’s curves clearly dominate.
On larger problems the exact DBN method’s performance quickly deteriorates.
Computational tractability requires the exact method to impose highly
restrictive in-degree constraints. These observations are consistent with the
heatmaps of Figures 3 and 4 in the main text.

###  D HPN-DREAM challenge details

We provide additional details for the methodology and results of the HPN-DREAM
challenge evaluation.

#### Data preprocessing.

The HPN-DREAM challenge data needed to be preprocessed before it could be used
by the inference methods. The choices we made during preprocessing most likely
affected the inference results.

Many of the time series contain duplicate measurements. We managed this by
simply averaging the duplicates. We log-transformed the time series since they
were strictly positive and some methods (SSPS and exact DBN) assume normality.
This probably made little difference for FunChisq, which discretizes the data
as part of its own preprocessing.

#### Predicted networks.

Figure 8 visualizes networks from two biological contexts in the HPN-DREAM
challenge evaluation. This gives a sense of how the different inference
methods’ predictions differ from each other. All of the predicted networks are
fairly different, though the SSPS and exact DBN predictions are more similar
to each other than they are to FunChisq. FunChisq predicts more self-edges
than the other methods.

In the BT549 cell line, the experimentally detected mTOR descendants include
receptor proteins that would traditionally be considered upstream of mTOR in
the pathway. The experimental results are reasonable due to the influence of
feedback loops in signaling pathways. However, the number and positioning of
the mTOR descendants highlights the differences between the coarse HPN-DREAM
challenge evaluation, which is based on reachability in a directed graph, and
the more precise evaluation in our simulation study, where we have the edges
in the ground truth network.

Figure 8: Prior and predicted pathways from the HPN-DREAM challenge. We show
pathways from two contexts: cell lines BT549 (top row) and MCF7 (bottom row).
The stimulus is EGF for both contexts. SSPS attained the best AUCROC of all
methods in the (BT549, EGF) context and the worst in the (MCF7, EGF) context.
The yellow node is mTOR; red nodes are the experimentally generated (‘‘ground
truth’’) descendants of mTOR.

#### HPN-DREAM AUCPR.

For completeness, we complement the AUCROC results of Section 3.2 with the
corresponding AUCPR results. Figure 9 shows AUCPR in bar charts, with an
identical layout to Figure 5.

AUCPR leads us to similar conclusions as those from AUCROC. SSPS dominates the
exact DBN method in 19 contexts and is dominated in 10. Both SSPS and FunChisq
dominate each other in 14 contexts. However, SSPS dominates the prior
knowledge in only 9 contexts, and is dominated in 21. As before, we conclude
that SSPS attains similar performance to established methods on this task.

Figure 9: A bar chart similar to Figure 5 except that it shows AUCPR rather
than AUCROC. See Figure 5 for details about the layout.

#### ROC and PR curves.

Figure 10 shows ROC and PR curves from our HPN-DREAM evaluation. We focus on
two representative contexts: cell lines BT549 and MCF7, with EGF as the
stimulus.

The bar charts in Figure 9 tell us that SSPS was the top performer in the
(BT549, EGF) context. The ROC and PR curves are consistent with this. SSPS
dominates the other methods in ROC and PR space. In contrast, SSPS was the
worst performer in the (MCF7, EGF) context. The curves show SSPS performing
worse than random.

The LASSO ROC and PR curves are interesting. Its ROC curves show nearly random
performance. Its PR curves are straight lines. Manually inspecting its
predictions yields an explanation: (i) LASSO gives nonzero probability to a
very small number of edges; (ii) that small set of edges results in a very
small descendant set for mTOR; (iii) that small descendant set is incorrect.

ROC Curves

PR Curves

Figure 10: ROC curves (top) and PR curves (bottom) from the HPN-DREAM
challenge. We show results for two contexts: cell line BT549 (left) and MCF7
(right). The stimulus is EGF for both contexts. Since SSPS is stochastic, we
show all 5 of its curves in each plot. The other methods are all
deterministic, and therefore only have one curve in each plot.