Scientific Model

Overview

phospho-velocity models the temporal dynamics of phosphosite intensities as:

\[ y_{s,c}(t) = f_{s,c}(t) + \varepsilon \]

where \(f_{s,c}(t)\) is a smooth latent trajectory for phosphosite \(s\) in cell line \(c\), and \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\) is measurement noise.

Gaussian Process Trajectory

The latent function \(f\) is modeled as a draw from a Gaussian Process:

\[ f \sim \mathcal{GP}(0, k(t, t')) \]

with a Matérn 5/2 kernel:

\[ k(t, t') = \sigma_f^2 \left(1 + \frac{\sqrt{5}\,r}{\ell} + \frac{5r^2}{3\ell^2}\right) \exp\!\left(-\frac{\sqrt{5}\,r}{\ell}\right), \quad r = |t - t'| \]

where \(\ell\) is the length-scale (default 30 min) and \(\sigma_f^2\) is the signal variance.

The Matérn 5/2 kernel is once-differentiable almost surely, making it the scientifically correct choice for computing a well-defined GP derivative (phosphosite velocity).

Phosphosite Velocity

Velocity is defined as the time derivative of the log2-intensity trajectory:

\[ v_{s,c}(t) = \frac{d}{dt} f_{s,c}(t) \]

In the GP-derivative approach (default), velocity is derived analytically from the joint GP over \((f, df/dt)\) using the derivative kernel. Let \(r = |t - t'| / \ell\):

\[ \frac{\partial k}{\partial t}(t, t') = -\sigma_f^2 \cdot \frac{5(t-t')}{3\ell^2} \cdot \left(1 + \sqrt{5}\,r\right) \exp\!\left(-\sqrt{5}\,r\right) \]

\[ \frac{\partial^2 k}{\partial t\,\partial t'}(t, t') = \sigma_f^2 \cdot \frac{5}{3\ell^2} \cdot \left(1 + \sqrt{5}\,r - 5r^2\right) \exp\!\left(-\sqrt{5}\,r\right) \]

The posterior mean and variance of \(df/dt\) at any \(t\) are obtained by conditioning on the observations via standard GP conditioning — yielding a full uncertainty estimate for velocity with no finite-difference approximation.

Observation Model

Following Robin et al. (2019), observation uncertainty scales inversely with phosphosite localization probability:

\[ \sigma_{\text{obs}} = \sigma_{\max} - (\sigma_{\max} - \sigma_{\min}) \cdot p_{\text{loc}} \]

where \(p_{\text{loc}} \in [0, 1]\) is the localization probability from MaxQuant.

An optional empirical variance model (ObservationModel.fit_variance_model) fits a Robin et al. 2019 exponential relationship between \(|\text{log-ratio}|\) and SD:

\[ \sigma(x) = A \cdot \exp(B \cdot \log|x|) \]

Bayesian Inference

In Bayesian mode (PyMC), priors are placed on the GP hyperparameters:

\[ \ell \sim \text{HalfNormal}(50), \quad \sigma_f \sim \text{HalfNormal}(\text{scale}), \quad \sigma_n \sim \text{HalfNormal}(1) \]

where scale defaults to 2.1 for orphan sites and is tightened for well-connected sites via the network prior (see below).

The NUTS sampler is used to obtain full posterior distributions over trajectories and velocities. Posterior velocity is computed analytically using the derivative kernel on each MCMC draw.

Network Prior

Kinase–substrate edges (from KinomeXplorer / NetworKIN / Creixell et al. 2015) are used as regularization weights:

\[ \lambda_s = \frac{1}{1 + |\mathcal{K}_s|} \]

where \(|\mathcal{K}_s|\) is the number of kinases targeting site \(s\). The prior scale on \(\sigma_f\) is set to \(2\lambda_s + 0.1\), so well-connected sites (small \(\lambda_s\)) receive a tighter prior and orphan sites a wider one.