Modelling sequential branching dynamics with a multivariate branching Gaussian process
Authors: Elvijs Sarkans
Sumon Ahmed
Magnus Rattray
Alexis Boukouvalas
Venue: Accepted by TMLR ()

Abstract: The Branching Gaussian Process (BGP) model is a modification of the Overlapping Mixture
of Gaussian Processes (OMGP) where latent functions branch in time. The BGP model
was introduced as a method to model bifurcations in single-cell gene expression data and
order genes by inferring their branching time parameter. A limitation of the current BGP
model is that the assignment of observations to latent functions is inferred independently
for each output dimension (gene). This leads to inconsistent assignments across outputs
and reduces the accuracy of branching time inference. Here, we propose a multivariate
branching Gaussian process (MBGP) model to perform joint branch assignment inference
across multiple output dimensions. This ensures that branch assignments are consistent and
leverages more data for branching time inference. Model inference is more challenging than
for the original BGP or OMGP models because assignment labels can switch from trunk to
branch lineages as branching times change during inference. To scale up inference to large
datasets we use sparse variational Bayesian inference. We examine the effectiveness of our
approach on synthetic data and a single-cell RNA-Seq dataset from mouse haematopoietic
stem cells (HSCs). Our approach ensures assignment consistency by design and achieves
improved accuracy in branching time inference and assignment accuracy.

Published in Transactions on Machine Learning Research (05/2023)
Modelling sequential branching dynamics with a multivariate branching Gaussian process
Elvijs Sarkans BIOS Health
ebujs.sarkons@@maalcoo com
Sumon Ahmed University of Dhaka
sumonSduu.ac.bl
Magnus Rattray University of Monchester
Alexis Boukouvalas PROWLER.Ko
aleris. Serrissboookkoovaae@@aaaooooom
Reviewed on OpenReview: https: // operreveeww.we.... net/ forum? id- 9KB0lstTg
Abstract
The Branching Gaussian Process (BGP) model is modification of the Overlapping Mixture of Gaussian Processes (OMGP) where latent functions branch in time. The BGP model was introduced a method to model bifurcations in single-cell gene expression data and order genes by inferring their branching time parameter. A limitation of the current BGP model is that the assignment of observations to latent functions is inferred independently for each output dimension (gene). This leads to inconsistent. assignments across outputs and reduces the accuracy of branching time inference. Here, we propose a multivariate branching Gaussian process (MBGP) model to perform joint. branch assignment inference across multiple output dimensions. This ensures that branch assignments leverages more data for branching time inference. Model inference is more challenging than for the original BGP or OMGP models because assignment labels can switch from trunk to branch lineages as branching times change during inference. To scale up inference to large datasets we use sparse variational Bayesian inference. We examine the effectiveness our approach on synthetic data and a single-cell RNA-Seq dataset from mouse haematopoictic stem cells (HSCs). Our approach ensures assignment. consistency by design and achieves improved accuracy in branching time inference and assignment accuracy.
consistent and
1 Introduction
Many algorithms have been developed to uncover trajectories from single-cell gene expression data (Saelens et al. 2019). Genome-scale single-cell experiments are destructive, so that time information lost, and tra- jectory inference provides a useful approach to associate cells with inferred time labels (pseudotimes) help recover gene expression dynamics. The resulting trajectorics are of great biological interest and branching trajectories are of particular interest in systems where cells differentiale into multiple alternative fates.
Pseudotime inference methods capture broad changes across high- hign-dimensional gene expression data. For example, popular approach is the Monocle algorithm (Qiu et al. 2017). which combines dimensionality reduction with minimal spanning tree estimation to order cells along trajectories. However, these algorithms do not model differences in branching dynamics across individual genes. In this paper, we consider the complementary problem of modelling geno-specific branching dynamics along pscudotime trajectorics. This is biologically important as it identifies the branching order of different genes, leading a better under- standing of gene networks involved in cellular differentiation. Two existing methods tackle the problem of

---

Published in Transactions on Machine Learning Research (05/2023)
inferring gene-specific branching times: Branch Expression Analysis Modelling (BEAM) QQu et al. 2017] and Branching Gaussian Process (BGP) (Boukouvalas et al. 2018). BEAM was the first method to tackle the branching time inference problem by using a relatively straightforward spline-based approach. Boukouvalas ct al. (2018) improved on the branching time estimation with BGP, espocially for carly branching genes and genes with low signal-to-noise. They achieved this improved performance by combining trajectory inference with branch label inference to better identify gene-specific branching times. The PscudotimeDE) (Song & Li, 2021) and Tradeseq (Van den Berge et al. 2020) methods have also been proposed to estimate which genes are branching. However, these methods do not infer the branching time and cannot be used in ordering the branching times of different genes.
BGP uses an approach inspired by the Overlapping Mixture of Gaussian Processes model (OMGP; Lázaro- Gredilla et al., 2012) to probabilistically label cells with the branches they belong to. OMGP is a mixture model developed for time-series data where each mixture component a Gaussian Process (GP) function and data points can be assigned to alternative mixture components. In the case of single-cell gene expression data, pseudotimes represent cells' progression through some developmental process and are analogous to time-series labels. Therefore, each cell can be assigned to the mixture components of OMGP based on pseudotime, an approach adopted by the GPfates method for trajectory inference [Lönnberg et al. 2017). However, in the standard OMGP, mixture components are independent and are represented by priori independent latent GP functions, whereas cellular trajectories are typically branching as cells begin in the same state but later differentiate into one multiple alternative states during a developmental process. Therefore, [Boukovvlas] et al. (2018) adapted the OMGP model by introducing dependence among the mixture components. Their BGP model uses a branching kernel that forces latent. GP functions representing mixture components to intersect at branching location and hence to become dependent on one another. Before the branching time cells are assigned to trunk function and after the branching time cells can be assigned to either of two child functions. The branching time is a parameter of the model that can be inferred to distinguish early, late and non-branching genes.
Trajectory inference methods such as Monocle can be used to provide cellular branch label information after some global cellular branching time inferred from the reduced dimension gene expression data. However, early-branching genes may be incorrectly assigned to the trunk according to this global assignment. The BGP model allows the branching time for specific genes to be carlier than the global branching time and new branch labels to be inferred after this time. [Boukouvalas et al. (2018) showed the BGP approach outperformed BEAM, especially with regards to such early branching genes for which no global branching information is available. However, in the original BGP model the branch assignments were made independently for cach gene and could potentially be inconsistent across genes, which makes little biological sense. Here, we extend the model to address this limitation. In the proposed multivariate BGP (MBGP) we perform joint inference across all genes in a way that ensures consistency in cell assignments. This increases the power of the method to infer branch labels, so that MBGP is better able to correct errors in the Monocle branch labelling and can assign cells to branches with much weaker prior label information than BGP. An added benefit is that single fit of the model on many genes (MBGP) is faster than many fits on a single gene (BGP) due to the label inference being shared.
In Section 2 we define the modelling problem. In Section 3 we demonstrate the limitations of the BGP model on single cell gene expression data. In Scction | we describe the MBGP model and in Section 5 we compare its performance to the original BGP model on both synthetic and single-cell gene expression data.
2 Problem definition
Consider observations Y E R associated with inputs t E R (without loss of generality, assumed to be monotonically increasing and in [0, 1I) with the following properties:
Inputs define a branching process for each dimension d. Prior to an unknown branching point by E [0, 1), i.c when t, < ba, data point Y.,d is associated with a "trunk" while after this point data points can be associated with one of two underlying bifurcation "branches" that generated it.
2

---

Published in Transactions on Machine Learning Research (05/2023)
The minimum branching point across dimensions defines the global branching point. After this point each input in is associated with one of two global bifurcation branches.
The assignment of data points to branches is consistent with the global branching across outputs (but unknown). That is, there exists global assignment of inputs la to branches, such that if ln is in branch A, then all outputs Yn,d with bn > by are also associated with branch A.
For datasets that satisfy the properties above, we want to solve two problems:
branching time inference: for each output [1,2,.., D} determine the branching time ba.
data point labeling: for each data point n E 2,,, N determine which branch belongs to. The formal definition is somewhat dense, but it tries capture simple idea illustrated in Figure
Y
Y
0.75
0.5
0.50
0.25
0.0
0.00
-0.5
0..5
-0.50
-1.0 bj
b2
0.5
1.0
0.0
Figure Illustaation of the problem definition. We plot a dataset with = 2 outputs over = 10 data points corresponding to equidistant inputs on [0, 1]: t 0, 0.11...,.., &10 1.0. The branching time inference problem is solved by b1 0.3 and by = 0.5. The data point labeling problem is solved via data point assignments indicated by the data point colours. The trunk is coloured green while the two branches (A and B) are coloured yellow and pink respectively. In a biological context, the data Y,, corresponds to a gene that branches early and whose expression rises in branch A and falls in branch B, while the data Y.2 corresponds to a gene that branches later, whose expression falls in branch A and rises in branch B. The cell assignments are consistent once both genes have started branch.
As discussed in the introduction, the abstract problem definition is motivated by the trajectory inference problem in single-cell gene expression data. In this application, Ya,d measures how much cell n expresses gene d at input fn (usually called "pseudotime"). The branching points represent when the expression of a gene begins differ between lincages (branches). The "consistent" data point assignment to branches captures the idea that different cells evolve via consistent trajectories.
We note that no "ground truth" branching points cell labels are available for real datasets, which makes the algorithm evaluation somewwat subjective. We will therefore augment the evaluation of MBGP to include synthetic datasets with known ground truth.
3 Branching Gaussian Process (BGP): Overview and limitations
We highlight some the key equations and techniques used in the original BGP method (Boukouvalas et al. 2018) that we will be building on in $ We show how a good model for solving the branching time inference, but makes errors the data point labeling problem and does not ensure data point labeling consistency across outputs.
3

---

Published in Transactions on Machine Learning Research (05/2023)
The outputs in BGP are modelled as independent. This makes the extension from to many outputs straightforward (see $.3 for details). However, as will see in $3.4 this mathematical convenience comes at cost to modelling capabilities.
3.1 Model specification
Let E Nx be single output data conforming to the full problem definition from We use latent GP functions to represent the trunk and branches of the data.
We first discuss the kernel construction that captures the requirement that the trunk and branch latent functions cross a single location the branching or bifurcation point. We then write down the model likelihood, and finally briefly discuss the inference.
3.1.1 Branching kernel
For clarity we restrict our attention to three latent GP functions f2, fa} that represent the trunk and branches of the bifurcation. We will sometimes also denote these three latent functions by (f.g. h. with f denoting the trunk. For convenience, an evaluation of latent functions N imputs (corresponding cells in y) is denoted by Fe E R3xN where fi corresponds to the i-th row (1 3))
f2.
Notation: we use Fe to describe the latent functions using OMGP-compatiblee matrix shapes. In Section 3.3 we will define reshaped but equivalent F that makes the variational inference equations simpler.
In order for the latent trunk and branch functions to capture bifurcations, we need to constrain them all to cross al a given input point. Yang al. (2016) described the joint probability density and the corresponding GP covariance for two latent functions constrained to cross each other given input. This was extended by Bookouvalas et al. (2018) to three latent functions crossing at single point. We bricfty ro-state the key governing equations.
We place the GP priors on f, 9: h, constrain them to all cross at the branching time b, and restrict the prior covariances to all have the same structure k:
f(t) ((,
(t) ((t), h(t) k(t,),
f(b) = g(b) = h(6).
Let us consider a vector of times and evaluate the covariance of the random variables ft), g(t), h(t). denoted by S. Since the individual covariance structures of f, g and h are identical, the diagonal blocks of S are given by Kig == k(t,, t)); whereas off off-diagonal blocks are given by
ke ktt, b) t2)
L
k(b, 6)
giving us the overall covariance
K
L
L
S
L
K
L
(1)
L
K
Figure 2 shows example of this covariance structure. The diagonal blocks show the covariance structure of kernel & (in this case, an exponentiated quadratic) whereas the off-diagonal blocks reflect the intersection of latent functions al time where the functions are constrained to be equal.
Note that the branching kernel can be extended to capture multiple bifurcations by extending the number of latent functions considered. This general case can be treated in much the same way as the single bifurcation case. However, here restrict our attention to modelling single bifurcation event.
4

---

Published in Transactions on Machine Learning Research (05/2023)
f
g
2.0
1.5
g
-0.5
h
Figure 2: Illustration of the covariance matrices used by Boukouvalas et al. 2018 The input contains 100 evenly distributed time points within the range [0, 10] for each latent function, and the global branching time is to 4. For the joint covariance structure latent functions f, and h are intersecting at the branching time b (left sub-panel) and the covariance matrix from which these three functions have been sampled (right sub-panel). The covariance matrix is evaluated at every point within the imput domain [0, 10] for all three functions.
have used an exponentiated quadratic. Three
3.1.2 Likelihood
Let Ze 10 be a binary indicator matrix determining the assignment N data points to precisely one the three latent functions. That is, each row of 2 has precisely one nonzero entry. The model likelihood given by
ply FC zd) = diag II)
(2)
We note that the extension multiple outputs is done by employing a different z for each gene, all of which are independent. This makes the likelihood factorise, that the extension is mathematically straightforward and we omit its discussion from this section.
As in dizaro-Gredille al. 2012
categorical prior placed on the indicator matrix Z
P P(Z9)
H 2 n..
where H E Nx3 is the data point assignment probability matrix. That is, H n,s represents the probability of cell being assigned to latent function i, so that Σ He n,i Ho.. outputs case there is different Ze and H for each gene.
for all n. We note that in the many-
3.2 Equivalence with OMGP
Whilst BGP was inspired by OMGP, the of somewhat. different notation means that the model likelihood equivalence is not obvious. Let us show it explicitly.
Equation
(4) from OMGP Lázaro Grodilla et al. 2012 defines
ply FC F Zo)
Ih Th (yn; Fin, Za.... Recall that Ze has precisely
nonzero entry in each
We

---

Published in Transactions on Machine Learning Research (05/2023)
denote its index for row to be in so that Znii. 1. Then
IIII (wns Pe 7 II N mi..::o)) i=i==
Pr o²)
πI
N
V/ni
Z Fi in
=N(y; diag 20
0²)
(4)
thus establishing our claim. Note that we are considering only three branches, but the equivalence would hold for more general branching kernels as well since the derivation above does not depend on the number of branches.
3.3 Inference
The BGP inference scheme is based on an important re-statement of the model likelihood in 2 Let the rows of Fe be given by {f. hf ht. Define a reshaped F := (f1, (11,91, h1, 92, 92, h2,...,
........N SN, WN,,,, hx) E R3Nx1,
Given Ze as above with Z = (as in $3.2,, define Z E {0, 1 to be a binary indicator matrix with everywhere except for one element on each row: 2n. 1)+, for n 1,2,........
We highlight that the important consequence of the above definitions is that diag 2 2= ZF since
3N
diag( (Z = Σ
Ze Fi, Fú = FF = Fs(n-1)+H = 2n, (
1=1
This gives an alternative formulation of
2
F : p F zo) = Nyy diag(Z" = Ny; 2F,00111
(5)
The log-likelihood of the model is not analytically tractable as it requires integrating out the indicator matrix Z. Boukouvalas et al.] (2018) used variational inference to derive an analytical solution for the lower bound as in Greaill et al. 2012 More precisely, they modelled the distributions of branch assignment Z and latent processes F as independent: q(,, = Q (F)); this sometimes referred to as the mean-field approximation. A standard variational lower bound on the marginal likelihood is derived using Jensen's 's inequality (see e.g. King & Lawrence 2006)
logp (B) 2 Eo(z) [logp (3) Z)] DEL IQ PP())
(6)
Boukouvalas et al. (2018) then chose to model data point assignments to latent branches Z as independent
3N Bernoulli-like variables for each data point. That is, Q == 11 - U.Z. where Un, i is the posterior
probability of data point being assigned to particular branch in 1,2,3 with == 3(n +) (probability
is 0 for all other outcomes). The branch can be either the trunk state when in S b, or one of the branches
after bifurcation when fn > b. By denoting == | [Z] they arrived at an important bound on the KL divergence for the above choice of Q(Z) (see Equations (18)-(19) from Boukouvalas el al. (018))
3N
Und Hn
DKL (Q(Z))P(() = Σ
Σ
Un log
where we have defined H as a counterpart to H from 3 in the obvious way. Next, F is integrated out derive an exact variational collapsed bound of the marginal log-likelihood logp(y) (Boukouvalas et al., 2018). The model's parameters can be estimated by maximising a variational bound on the log-likelihood.
Boukouvalas et al. 2018)) also derived a sparse approximation bound of their model using inducing points that allows the BGP model to scale up to larger datasots.
An important aspect of the BGP is that the branching time posterior probability is not available analytically and is calculated numerically using the approximate marginal likelihood evaluated at a uniformly spaced set of candidate branching points. This approach would scale exponentially with the number of branching times in the MBGP model that we introduce below and therefore we use an alternative gradient computation estimate only the mode of the branching posterior.

---

Published in Transactions on Machine Learning Research (05/2023)
3.4
Limitations of the method
BGP is a good model for solving the branching time inference problem as evidenced by comparisons with the relevant methods, ie. BEAM (Qiu et al. 2017) and the Mixture of Factor Analysers (Campbell & Yau 2017), see BBouoovvllas al. (2018).
Having said that, its applicability is hindered by its shortcomings in solving the data point labelling problem. The BGP model assigns data points to branches output-by-outph without ensuring consistency across outputs (see section $2). That is, a data point may be assigned to branch A for one output and to branch B for another output.
Notably, when applied to biological datasets analysing early branching genes that branch before the global branching time identified by a pseudotime algorithm, then no prior information on cell branch assignment prior to the global branching time is available. In this case it is likely to be beneficial to leverage data from multiple genes to consistently solve the cell labelling problem, since there will be high uncertainty given data from only one gene. We note also that the gene-by-gene inference is computationally inefficient.
To explore these issues empirically, have applied the BGP method on a single cell mouse hematopoietic stem coll dataset (Paul al.) 2015). The data contain 4423 cells and a global pscudotime and branching structure was inferred using Monocle (Qiu et al., 2016). Details of the experimental setup are given in the supplementary material.
(a) gene A
(b) gene B
Gene Expression
3
2
1
0
-1
0.2
0.4
0.8
1.0
0.2
0.6
0.8
0.5
0.5
0.0
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Pseudotime
Figure 3: Mouse haematopoietic stem cells (Paul et al. 2015 Inconsistent cell assignment between the pair of early branching genes MPO (gene A) and CTSG (gene B). The 1000 randomly sub-sampled cells that have been used for inference are coloured according to the posterior global branching pattern derived using the BGP algorithm (yellow means more likely be in branch A, blue means more likely to be in branch B). The points are used to represent the cells have been assigned to the same branches, i.e. consistent across these two genes. On the other hand, the triangles represent the inconsistent cell assignment, i.e. cells that have been assigned to different branches across these two genes.
Boukouvalas et al. (2018) tried to mitigate the risk inconsistent cell assignments by placing a strong global prior over the cell assignments, but the problem can remain even when strong prior is used. For instance, Figure X shows the posterior cell assignment of randomly sub-sampled 1000 cells for two early branching genes MPO and CTSG. To avoid the inconsistent cell assignment, very strong prior confidence 0.95

---

Published in Transactions on Machine Learning Research (05/2023)
used in both cases. This strong prior assumes that cells after the global branching point in pseudotime ordering are unlikely to switch their branching assignment away from the prior. Among the 1000 cells that are used in the inference, we have found that 88 cells (triangle markers) have inconsistent cell assignments, i.c. have been assigned to different branches in the posterior. The posterior cell assignment inconsistency on the eight biomarkers of HSCs that show very strong evidence of branching has been summarised in Table Posterior cell assignments have been calculated using a high confidence (0.80) for all eight genes. For each gene, different numbers of cells show strong evidence of being assigned to the different branches. The cells may actually belong to particular branch they have become biased towards that branch in case of the gene under consideration. For each pair of genes, we identify the number of cells that have been assigned to different branches in posterior by each gene. We find that among the 1000 cells used in the inference of BGP, 468 cells show inconsistency across all eight genes.
Table 1: Number of cells (among 1000 cells used in the inference) inconsistently assigned to different branches by BGP for each pair of the eight biomarkers of hematopoietic stem cells (HSCs).
MPO
PRTN3
CALR
CAR1
CAR2
GSTM1
ELANE
CTSG
88
99
75
230
156
158
76
MPO
81
127
288
201
213
115
PRTN3
147
299
219
226
136
CALR
223
146
147
72
CAR1
191
205
182
CAR2
159
GSTMI
105
The second limitation of the BGP model is computational complexity. Since BGP estimates the branching time posterior by using the normalised approximated likelihood calculated at a set of candidate branching points over the imput domain, the computational requirements are high when analysing many outputs. Although the model inference can be done in parallel (per output), this still requires signíficant compute. In our experiments reported in Section 5 BGP on 1000 data point dataset takes around one hour per output for total of hours across all 8 outputs whilst MBGP takes around 30 mins for all outputs" Bookouvalas et al. (2018) had to sub-sample the data points in order to minimise the per-output computational time, but the computational requirements still remain very high. This exceptionally high computational load makes the BGP model impractical in many instances.
4 Multivariate Branching Gaussian Process (MBGP)
We propose Multivariate Branching Gaussian Process (MBGP) model that performs inference jointly across outputs. This approach addresses BGP's data point label inconsistency issue, is computationally faster and improves the statistical power available to infer branch assignments by jointly considering all outputs.
The MBGP model builds on the BGP model discussed $3.1 It uses the same branching kernel (sce $3.1.1 It also uses the same likelihood from (5, but with one modification the Z and H assignment matrices and prior probabilities and hence develop new variational approximation for the marginal likelihood. This change to the model addresses the consistent data point labelling problem from $ by ensuring consistent data point labelings across all outputs.
shared across all outputs. This means the model no longer factorises across outputs
Another change in MBGP its treatment of branching times. Recall that $3.3 described how the branching times were obtained by comparing approximate marginal likelihoods on grid of candidate points. As this does not scale to many outputs due to the curse of dimensionality, we treat branching times as model parameters that are optimised along with other parameters.
SSstemm configuratio:: Intel(R) Core(TM) 17-1065G7 CPU 1.30Hz 1.50 GHz, 16 GB RAM, 64-bit operating system, x64-hased processor.
8

---

Published in Transactions on Machine Learning Research (05/2023)
4.1 Model and variational bound
The likelihood factorises for multiple independent outpuls so that the exlension of the single-output likelihood in to multiple independent outputs straightforward
p PY 2)
II N oF d=1
where D is the number of outputs, ya the dth column of the data matrix Y (y1,..., yn) and Fa denotes the tt column vector from the set of latent GP functions F = {Fi. FD}
Crucially, the association matrix Z is shared by all outputs with categorical prior (equation 30 placed over Z. We treat the branching time for each output by as an explicit parameter of this prior. The prior over the latent functions Fa for each output a branching GP also parameterised by the branching time and kernel parameters 0a,
N 3N
P (Z II II
H (ba) Zn.i 2n.)
Fabba GP (0, k(Oa))ba)) for d = 1.2...,..
The categorical prior is adjusted in an output-specific manner so that data points are either assigned to the trunk or the two branches depending on the branching time. That is, 3(x)) = t < h, while
[1))),.3mm (()))) = ih,, - hnd otherwise, with all other entries in row n set to zero. The variational approximation for the assignment probabilities Q(2) has the same structure
Q (Z
II II (ba Zm.i
(9)
with 1)41 = 1 if en < b, while [U(b,,, 3(n 1)42. U (b).. is the inferred probability that data point is assigned to the first branch.
[Mn, 1 - un]] otherwise, where un
The log likelihood terms also depend on the branching parameter. If we denote - (b1,...b)) e Rº and integrate over QZ b), then by using Jensen's incquality in the standard way bring log inside the integral we have
D
>
EQ(2)ba [logp (ya) Fa, Z, ba)ll
Σ
DKL IQ (Z(ba) (PP
D
2
ND
(2TT2)
Σ
1
(yl ya + F. A.Fa U.U.U
20
d=1
DKL IQ PP
(10)
1=1
where U₂ & Eo(z)h) (Z) and A. A EQ(Z)))) ZZ are the variational expectations conditioned on the output-specific branching time ba-

---

Published in Transactions on Machine Learning Research (05/2023)
We note that Ua is straightforward to compute given equation 9 but Aa is not. Computing it explicitly is computationally prohibitive, so establish key identity.
N
(Ad)as = Eo(z)ba)) [(m2))s]
= (Z(0)
Zk..
(separate out the =j case) = (b(()(ba))
Zu Zed
(1 (1- + (16b))
Zi, Zui
Oet
J=1
N
(Zz Zen as it's binary) (l())
Zk.iZus
(1-61g) ) (2(
Σ
Oeg
k=1
N
(Zx = when + (b()(ba)) (ba)
Σ
Zu.i
dig-
(11)
kk=1
Therefore, A is diagonal and the individual terms are now easy to compute using equation and equation I1
Finally, note that counterpart of equation 7 holds via straightforward computation
N 3N
Unii(ba) (n,((ba)
DKL (Q(Z (() =
Σ>
Un (ba)
(12)
By plugging equation 12 and equation 11 into equation T0 we finally oblain a computationally feasible variational lower bound on the likelihood logp (Y b). Finally, we obtain a sparse collapsed variational bound by integrating out the latent functions F and introducing inducing point approximation. This is relatively straightforward and well-known, so only provide the full details in the appendix.
4.2 Inference
Although it is possible to estimate the lengthscale and process variance parameters jointly with the branching points, we prefer to estimate them prior to inferring branching locations. We set the branching location of every output at the starting of time (i.c. b - 0) and estimate the lengthscale and the process variance parameters for each output by maximum likelihood. After estimating lengthscale and process variance, the only remaining paramelers to estimale are the branching time vector and data point assiguments (Z). Using this approach we assume that the bifurcating functions have similar smoothness and variance before and after the branching time. This allows to simplify inference by reducing the number of effective parameters that need to be estimated. We also employ multiple restarts to reduce the effect of local minima.
By doing gradient search we no longer infer the posterior over branching times as in the original BGP model. If a full Bayesian posterior over branching times is required then an MCMC approach can be adopted, although that would be computationally demanding.
5 Experimental results
We demonstrate our approach by fitting both MBGP and BGP models to noisy synthetic data. By fitting models to realistic noisy samples, will see that MBGP substantially improves the data point assignment and improves the branching time learning when compared to BGP. We will also fit the MBGP model to real data, namely, single-cell RNA-seq data from mouse hematopoietic stem cells (Paul el al. 2015) and see how it resolves the cell label inconsistency problem whilst retaining biologically sensible fits.
5.1 Synthetic noisy data
In many real-world applications, there is no ground truth for data point assignments and branching times. This can make it difficult compare mothods a quantitative manner. In order address this, we will draw samples from fixed MBGP model where both the ground truth branching points as well as data point assignments are known. This allows us to provide quantitative comparisons of MBGP and BGP.
10

---

Published in Transactions on Machine Learning Research (05/2023)
gene 0
gene 1
2
0.5
Synthetic Gene Expression
-1
-1.0
-2
gene
gene
1.5
-0.5
10
1.5
0.5
2.0
3.0
025 0.50
0.75
0.00
0.25
0.50
0.75
1.00
Pseudotime
Figure 4: An example of the MBGP model fit on noisy synthetic data with 4 outputs. The data is represented in a scatter plot with inputs along the x a-axis and colours indicating branch assignments. The model's mean predictions are the solid lines. The vertical blue dotted lines are the true branching times. The vertical magenta dashed lines are the learned branching times.
5.1.1 Setup
We sample noisy data from realisations of MBGP and reject samples where the latent. branches cross in than location. The MBGP that we sample from has branching points set at target locations distributed equidistantly on [0, 1]. Each sample consists of 100 data points. In all experiments that follow we employ the same data point label setup and use the same fitting procedure. This holds true for both MBGP as well as BGP.
We initialise data point label probabilities for latent branches by sampling uniformly from [0.5, 1] for the true label. We then set the data point label prior probabilities to 0.5 for inputs in [0, 0.8] (that is, no prior knowledge) and then 0.8 for the true label for inputs in [0.8, 1], which we will refer to as the "informative prior'. Note that even this "informative prior' is weaker than some of the priors explored in the original BGP paper, see 3..4 for details.
Unless explicitly stated otherwise, the fitting procedure is multiple restart from 4 different branching points, namely, 0.2, 0.4, 0.6, 0.8 applied uniformly to all outputs. For each restart, we perform gradient descent along the evidence lower bound surface using L.BBGS--B. We demonstrate a typical model fit in Figure
5.1.2 MBGP assigns data point labels more accurately
One of the key hypothesised advantages of using model that considers all outputs at the same time is that by avoiding assignment inconsistencies, namely the biologically impossible situation where cell A is assigned to one branch in gene 0 and a different branch in gene 1, should improve cell label assignment. accuracy as it uses information from genes.
In order to test the hypothesis we drew 10 random noisy synthetic samples as described above. We then constructed MBGP and BGP models with the ground truth branching times. In contrast. to the other
11

---

Published in Transactions on Machine Learning Research (05/2023)
synthetic data experiments we kept the branching times fixed to their ground truth locations and only optimised the other parameters, most notably the data point assignment variational parameters.
We found that MBGP achieved a better mean percentage of correct data point labels at 96. 0.77% compared to BGP's 0.5% (the = value is the standard error). For more detail, see Figure We note that achieving 100% accuracy on data point assignment challenging task given the high dataset. noise.
100
92
mBGP
BGP
Figure 5: Boxplots of data point assignment correctness averaged over all outputs for each of the 10 samples in MBGP vs BGP. Reported as percentage, so that 100% is the best score you can hope for. We can see MBGP outperforming BGP and in fact being close to the theoretical highest score.
We also looked at the distribution of incorrect data point labels over inputs. It showed that MBGP performs similarly to BGP at early and late inputs, but holds a considerable advantage in the range of [0.3, 0.8). see Figure G The equivalence on [0.8, 1] is not surprising as both models use the informative data point label prior for these inputs. The full per-sample results are reported in the appendix.
Overall, we note MBGP's impressive performance on the data point labeling problem. It suggests that alongside accurate model fitting procedure for branching points MBGP would constitute powerful tool for branching time and data point assignment analysis.
5.1.3 MBGP improves branching time estimates
The final synthetic data experiment looks the overall performance of MBGP vs BGP when we're learning both data point assignments and branching locations.
It also explores the inconsistently labelled data point problem in BGP. This happens when data point is assigned to different. branches in e.g. gene 0 and gene 1, which is biologically unrealistic.
We draw 20 random noisy samples and fit MBGP as well as BGP models using the setup described in Section 5.1.1 We measure (i) the root mean square error (RMSE) between the learned and ground truth branching points, (ii) the number of correctly learned data point labels, and (iii) the number of inconsistently labeled data points by BGP.
The overall RMSE and correct data point label performance for both models is summarised in Table which shows MBGP outperforming BGP on both measures. We further see in Table 3 that the difference in performance is significant for both RMSE as well as the percentage of correctly learned data point labels.
We note however that the absolute values of RMSE are quite high for MBGP and we hypothesise that this is the key factor the drop of correct data point label percentage from 96.7% in Section 5.1.2 (where we used the ground truth branching locations) to 86.4% here. It therefore seems that improvements learning branching locations would be worthwhile future research direction.
12

---

Published in Transactions on Machine Learning Research (05/2023)
All incorrect labels
All incorrect labels
14
12
2-
0.2
Figure 6: Histogram of all incorrect data point assignments across all samples and all outputs. The z-axis is the input of the incorrectly labeled data point. The y-axis is normalised to represent the incorrect data points in each bin a percentage of all data points in that bin. We see that MBGP significantly improves on the data point labeling task and hardly makes any mistakes after 0.5.
For each sample, compute (i) the root learned branching points (lower is better), and (ii) % the percentage of correctly learned data point labels (higher is better). We see that in both cases MBGP outperforms BGP, but we note that we see relatively high standard deviations. We believe this down to the high noise level in the synthetic data.
Table 2:
square error (RMSE) between the true and
Mean
Standard deviation
RMSE
0.178
0.074
RM RMSEEGP
0.222
0.070
BMP 86.4
5.1
%BP
80.4
4.1
Finally, we report that BGP on average has 32.0 + 9.8 inconsistently assigned data points out of total of 100.
We note that duc to the high noise level in random samples, some of them can be extremely difficult to model, see the appendix for some examples.
5.2 Hematopoiesis single-cell RNA-seq
We apply the MBGP model on the gene expression of mouse hematopoictic stem cells used in Paul et al. (2015). As in Boukouvalas et al. (2018) we have used the Wishbone algorithm /Stty al. (2016) to derive the pscudotime of cach cell as well the global branching structure. The model also readily applicable using pseudotime from other methods such as Monocle (Qiu et al. 2017
Boukouvalas et al. (2018 have given branching time network of eight genes (Figure (a) for which they found the highest evidence of branching. In the branching time network shown in Figure 7 (a), the most probable branching time is annotated with each gene and directed edges are used to represent pairwise ordering of genes based on the most likely gene-specific branching time. [oukouvalas et al. (0018) have grouped genes based on the branching time order relationships. For instance, both PRTN3 and CTSG are early branching genes and branch before ELANE, GSTM1, CAR2 and CAR1. On the other hand, MPO and CALR (yellow group) also early branching genes and branch before GSTM1, CAR2 and CARI, but
13

---

Published in Transactions on Machine Learning Research (05/2023)
Table For each sample, we compute (i) ARMSE: == RMSE RM where RMSE is the root mean square error between the true and learned branching points, and (ii) A% == ((MMGP) % (BGP) where % denoles the percentage of correctly learned data point labels. In both cases, values above mean that MBGP is outperforming BGP, whereas values below show the opposite. We report the mean and standard error over the 20 samples. As the reader can see, MBGP outperforms BGP.
Mean
Standard error
A RMSE
0.0442
0.0193
A% %
5.98
1.40
not necessarily before ELANE. Finally, the third group consists of genes ELANE and GSTMI (blue group) that branch before the gene CARI. The genes CARI and CAR2 have the latest branching time.
(b) Estimated branching time
PHINS(O.08)
CTSGIO
MPOIC.061
CALRIOLOS;
(a) Branching time gene network
CT50
2
BLANE
Gone Expresion
1-
2-
B4014
ELANEI0.211
CAR2IC 18)
GSTM2100201
CAR1(028))
MPQ
bat38
Brunching Time
Pscudotime
Figure 7: Mouse haematopoietic stem cells (Paul al. 2015 : (a): Estimated gene branching time network from Boukouvalas et al. (2018) the most probable branching time of cach gene is shown. The directed edges are used to represent the pairwise ordering of genes based on their branching time. The edge colours are representing group of genes having the same later branching genes. (b): The MBGP fit. The cells (points) are coloured according to the global branching pattern derived using the Wishbone (Setty et al. 2016) algorithm. The dashed line in each sub-figure represents the estimated branching time for cach gene, which is also depicted inside the parenthesis in each sub-figure caption.
Although BGP finds plausible ordering of the genes for this data, we showed in Table I there are large number of inconsistencies of cells to branches when using the original BGP algorithm. Therefore we have applied the MBGP model on the expression profiles of these eight genes which ensures consistency of cell to branch assignments. To speed up inference, we have randomly sub-sampled the data and have used 1000 cells for the inference and have used the sparse model with 60 inducing points.
Figure 7 (b) shows the MBGP model fit on this network of genes. In cach sub-figure, the dashed line represents the inferred gene-specific branching time. We reproduce the same pairwise ordering of genes published in BBououvalas et al. (2018) wherein it was found to be biologically meaningful. If we compare the inferred gene-specific branching times from Figure 7 (b), we see that our model has successfully reproduced the same branching time network of Figure (a).
5.2.1 MBGP improves cell assignment accurately
In Appendix C we show examples of cell label inconsistency when applying the BGP algorithm the mouse haematopoictic stem cell data. In contrast, MBGP ensures by construction no such inconsistencies can occur.
14

---

Published in Transactions on Machine Learning Research (05/2023)
Additional while MBGP ensures the cell assignment consistency, it also improves the posterior cell assign- ments. We have analysed the posterior cell assignment by MBGP which follows the global branching pattern provided as prior. As the model has learned the gene-specific branching points for a number of genes earlier than the global branching point (main paper Figure 6(b)), therefore, few cells from the trunk state become biased towards either of the branches to facilitate the identification of these pioneer genes. We have sum- marised this result in Table 4 where a posterior cell assignment confidence of 0.80 is used determine cell has been assigned to branch.
The number of cells switches states from prior to posterior as a result of the MBGP model inference for the mouse haematopoietic stem cell data (Paul et al. 2015). A posterior cell assignment confidence of 0.80 is used to determine a cell has been assigned to a branch. The model learns the genc-specific branching points for a number of genes earlier than the global branching point (main paper Figure 6(b)), therefore, few cells from the trunk state become biased towards either of the branches to facilitate the identification of these pioneer genes.
Table 4:
Number of cells in prior
Number of cells posterior
Trunk
278
187
Branch 1
423
466
Branch 2
299
347
6 Conclusion
The BGP model was developed to infer branching events from single-cell gene expression data. We uncovered its main limitations, namely inconsistent data point assignments and heavy computational requirements. To mitigate these limitations, we have presented the multivariate BGP (MBGP) model that extends BGP by sharing data point assignments across outputs (genes). On both simulated and mouse hematopoictic stem cell gene expression data we have shown that the MBGP model achieves similar or better branching time accuracy to the original BGP model whilst addressing the assignment inconsistency issue, making model inferences more biologically relevant.
In Section 5.1.2 we showed that when the MBGP model does not have to infer the branching locations, it performs the cell assignment task with impressive accuracy. However, we have also seen in Section 5.1.3 that it does not reach the same accuracy when it also has to learn the branching locations. Thus exploring improvements in the learming of branching locations seems like a valuable future research direction.
Although MBGP is faster than the original BGP formulation, inference is still time-consuming for large number of outputs since GP function has to be learned for each output. This also ignores dependencies between gene expression profiles from different genes. Future work could develop multiple output formula- tion where individual genes are functional combinations of a smaller set of latent GPs (Wilson et al. 201 1) This would allow our approach to scale better to many thousands of outputs, making the approach more practical for large scale single cell branching analysis.
Finally, we note that during the review process of this paper new method, scFates Faaure et al. 2022). released that also aims to identify carly-branching genes. Direct comparisons between MBGP and scFates are of obvious interest, but have not been evaluated in this paper.
References
Alexis Boukouvalas, James Hensman, and Magnus Rattray. BGP: identifying gene-specific branching dy- namics from singlo-cell data with branching gaussian process. Genome biology, 19(1):65, 2018.
Kieran R Campbell and Christopher Yau. Probabilistic modeling of bifurcations in single-cell gene expression data using bayesian mixture of factor analyzers. Wellcome open research, 2, 2017.
15

---

Published in Transactions on Machine Learning Research (05/2023)
Louis Faure, Ruslan Soldatov, Peter V Kharchenko, and Igor Adameyko. scFates: scalable python package for advanced pseudotime and bifurcation analysis from single-cell data. Bioinformatics, 39(1), 11 2022.
Nathaniel J King and Neil Lawrence. Fast variational inference for gaussian process models through kl-correction. In European Conference on Machine Learning, pp. 270-281. Springer, 2006.
Miguel Lázaro-Gredill., Steven Van Vaerenbergh, and Neil Lawrence. Overlapping mixtures of gaussian processes for the data association problem. Pattern Recognition, 45(4):1386-133) 2012.
Tapio Lönnberg, Valentine Svensson, Kylie R James, Daniel Fernandez- Ruiz, Ismail Sebina, Ruddy Mon- tandon, Megan SF Soon, Lily G Fogg, Arya Sheela Nair, Urijah Liligeto, et al. Single-cell rna-seq and computational analysis using temporal mixture modelling resolves th1 th1/tfh fate bifurcation in malaria. Science immunology, 2(9), 2017.
Franziska Paul, Ya'ara Arkin, Amir Giladi, Diego Adhemar Jaitin, Ephraim Kenigsberg, Hadas Keren-Shaul Deborah Winter, David Lara-Astiaso, Meital Gury, Assaf Weiner, el al. Transcriptional heterogeneity and lineage commitment in myeloid progenitors. Cell, 163(())1663 1666-77766111777 2015.
Xiaojic Qiu, Andrew Hill, Yi-An Ma, and Cole Trapnell. Single-cell mrna quantification and diffcrential analysis with census. Nature Methods, pp. 309 315, 2016.
Xiaojie Qiu, Qi Mao, Ying Tang, Li Wang, Raghav Chawla, Hannah Pliner, and Cole Trapnell. Reversed graph embodding resolves complex single-cell trajectories. Nature Methods, 14(10):979, 2017.
Wouter Saclens, Robrecht Cannoodt, Helena Todorov, and Yvan Sacys. A comparison of single-cell trajectory inference methods. Nature biotechnology, 37(5):547-5533 554, 2019.
Manu Setty, Michelle D Tadmor, Shlomit Reich-Zeliger, Omer Angel, Tomer Meir Salame, Pooja Kathail, Kristy Choi, Scan Bendall, Nir Friodman, and Dana Po'er. Wishbone identifics bifurcating developmental trajectories from single-cell data. Nature biotechnology, 34(6):637, 2016.
Dongyuan Song and Jingyi Jessica Li. Pseudotimede: inference of differential gene expression along cell pscudotime with well-calibrated p-values from singlo-cell ma soquencing data. Genome biology, 22(1): 1 25, 2021.
Koen Van den Berge, Hector Roux de Bézieux, Kelly Street, Wouter Saelens, Robrecht Cannoodt, Yvan Saeys, Sandrine Dudoit, and Lieven Clement. Trajectory-based differential expression analysis for single- cell sequencing data. Nature communications, 11(1):1 13, 2020.
Andrew Gordon Wilson, David Knowles, and Zoubin Ghahramani. Gaussian process regression networks. arXiv preprint Xic::11100444: 2011.
Jing Yang, Christopher A Penfold, Murray R Grant, and Magnus Rattray. Inferring the perturbation time from biological time course data. Bioinformatics, (2(19):2956-2964 2016.
A Sparse collapsed variational lower bound
We follow similar steps Sizaro-Gredilla et al. (2012) to derive a sparse variation bound from the full data variational bound derived in Section 4.1 and shown in equation ??
16

---

Published in Transactions on Machine Learning Research (05/2023)
Starting from the bound in equation ??, integrate out the latent GP functions F to derive the collapsed variational bound for the marginal log likelihood.
logp (Y)
2
ND
log log((2 (2TO²)
Y
2
20²
D
D
log
K s51bal
2
d=1
D
Σ
log Ado - + K 17//val
d=1
D
+
Yd Un AAro Kj(ba)
UJ yd
2
d=1
Σ
DRL. (2) |((b),)
(13)
d=1
where Kftha denotes the covariance matrix dependent on the output branching parameter ba.
The sparse variation lower bound for our model Ls is
ND
L A
log(2 ?) 20
Y
2
D
log Pal - col
2 d= |
1
[r (AaK tr (AcK fur
SucK Kunia Kufa)
20
Σ
DDKL
(14)
d=1
where we have defined
P. è 1 L K of A,K fs o0-2
ca 2 Ri L K Ksffba U -2
and 8 is the set of inducing points and the Choleksy factors La and Rg are defined as
Kuns Lab
Pa R.
B Detailed experimental results on synthetic data
B.1 Labelling performance
The full sam-s results are reported in Table 5 We have already reported the key summary statistics in the main part of the paper, so we omit further discussion. We also show in Figure g an example of the MBGP model fit on the noisy synthetic data. We note again MBGP's impressive performance on the data point labeling.
17

---

Published in Transactions on Machine Learning Research (05/2023)
gene
gene
gene 2
gene 3
15
10
0.5
0.0
0.5
1.0
-1
gene
gene
gene 6
gene
0.5
0.0
00.
-1.0
-1.5
0.5
-2.0
-2.
gene
gene
10
0.0
-0.5
0.0
-1.0
0.5
Figure 8: An example of the MBGP model fit noisy synthetic data with 10 outputs. The data is represented in a scatter plot with input along the axis and colours indicating branch assignments. The model's mean predictions are the solid lines. The vertical magenta dashed lines are the true branching times we do not learn branching locations in this example.
We draw the reader's attention the data points around the true branching point and suggest that it is a highly nontrivial problem to assign them to branches given the noise levels.
18

---

Published in Transactions on Machine Learning Research (05/2023)
Table 5: This table contains the per-sample results for MBGP and BGP. Each row corresponds a different 10 output noisy synthetic datasct. inCBGP only applies to the BGP model and describes how many data points have inconsistent assignments across the outputs, that is, for how many data points do assign them to one branch in output A and to a different branch in output B? We note that lower is better and 100 is the absolute worst result you could get while 0 is the best result. % the percentage of correctly learned data point labels (higher is better). In this mcasure, we ignore inconsistent data point assignments and simply average over putputs.
We can see that the MBGP model is able infer data point assignment with higher accuracy, although both models are able to achieve accuracy > 90%. We suspect MBGP outperforms BGP due to its multivariate nature allowing it to more accurately learn data point assignments around branching locations where there is significant uncertainty.
ineRGP
%(MBGP)
%(BGP)
0
22
95.3
91.3
1
24
94.4
91.5
2
18
99.6
93.8
3
26
95.2
91.3
4
23
94.7
90.3
97.6
93.0
6
18
97.2
92.7
7
16
99.3
95.1
8
27
99.2
92.6
9
31
94.3
90.1
B.2 Details on synthetic data branching performance
We provide more details the synthetic data experiments discussed in Section 5.1.3. We remind the reader that the final synthetic data experiment looks at the overall performance of MBGP vs BGP when we're learning both data point assignments and branching locations.
The per-sample RMSE and correct data point label performance for both models is provided in Table 6 We have already reported the key summary statistics in the main part of the paper, so we omit further discussion. However, we will now look two different regimes: "casy" and "hard' infer branching data, corresponding to samples 2 and 14 respectively in Table
The model fits on sample 2 from Table f are discussed in Figure 9 and Figure 10 for MBGP and BGP respectively. See the captions for a discussion.
The model fits on sample 14 from Table displayed in Figure which pertains to MBGP, and Figure 12 which pertains to BGP. See the captions for a discussion.
Finally, for completeness, we also provide (i) the distribution of incorrect data point labels across input, see Figure 13 and (ii) the mean percentage of correct data point labels, see Figure 14
C Mouse haematopoietic stem cells: Inconsistent cell assignment by BGP
In the BGP model, the inference is performed independently per gene which causes inconsistent cell assign- ments across different genes. To minimise cell assignment inconsistency prior probability of 0.95 as suggested by BBoukuvvlass et al. (2018), but the problems remains to great exlent. Figure 15 shows the expression proliles as well as BGP fits of eight genes where triangle markers are used to indicate the cells that have been assigned to different branches incousistently. Among the 1000 cells used the inference of BGP, 468 cells show inconsistency across all cight genes. would have included more genes this number would have been greater. Thus the cell assignment reported by the BGP is not informative and likely to be misleading.
very strong cell assignment
19

---

Published in Transactions on Machine Learning Research (05/2023)
gene 0
gene
gene
gene 3
2
3
2
0
-1
gene 4
gene
gene 6
gene 7
1.5
2
2.5
0.0
10
2.0
0.5
0.5
15
0.0
10
-05
11.5
-1
0.5
-1.0
0.0
2.0
gene 8
gene
1.0
1.5
2.0
2.5
-1.5
3300
0.0
0.0
0.5
Figure 9: This Figure shows the MBGP model fit on sample from Table G The data is represented a scatter plot with inputs along the æ axis and colours indicating branch assignments. The model's mean predictions are the solid lincs. The vertical blue dottod lines are the true branching times. The vertical magenta dashed lines are the learned branching times.
Overall the model does good job of learning the true branching times and thus does well on data point assignment: 4MBGP 93.7% vs %BGP 81.3%, see Table 6 We hypothesise that this could be due to the carly branching outputs, namely, output 0 and output 1 providing useful evidence for data point assignment that the her outputs then benefit. from. Note, in particular, how MBGP learns the correct branching time for output 7 whereas BGP learns incorrect dynamics, see Figure 10
20

---

Published in Transactions on Machine Learning Research (05/2023)
0.
Gene
Gene
2.0
Gene
0.5
Gene
0.0
15
0.0
10
0.5
1.0
115
-1.5
2.0
0.00
025
0.50
0.75
10
0.00
0.25
0.50
0.75
1.0
0.00
0.55
0.50
0.75
0.00
0.25
0.50
0.75
1.00
15
Gene 4
Gene 5
Gene
Gene 7
25
0.0
10
40
0.5
0.5
11.0
0.0
10
0.5
-1.5
0.5
0.0
0.00
0.25
0.50
0.75
0.00
0.25
0.50
0.75
100
0.50
0.75
0.00
0.25
0.50
100
0.5
Gene 8
Gene 9
-1.0
0.0
-1.0
-1.5
-3.0
0.00
025
0.50
0.75
10
0.00
0.55
0.50
0.75
10
This Figure shows the BGP model on sample 2 from Table 6 The data is represented scatter plot with input along the x axis and colours indicating branch assignments. The model's mean predictions are the solid lines. The vertical blue dashed lines are the true branching times. The vertical magenta dashed lincs are the learned branching times. See Figure for the key discussion points.
Figure 10:
21

---

Published in Transactions on Machine Learning Research (05/2023)
gene
gene 1
gene 2
gene
2
1
1
0
-1
-2
-2
-3
-1
gene
gene 5
gene
gene
1
1.
1.5
2.0
10
10
15
0
0.5
0.0
0.5
0.0
0.5
0.0
1.0
0.5
-2
1.0
1.5
gene
gene 9
10
0.5
0
0.0
-2
0.0
0.5
10
0.0
0.5 10
Figure 11: This Figure shows the MBGP model fit on sample 14 from Table The data represented in a scatter plot with input along the axis and colours indicating branch assignments. The model's mean predictions are the solid lincs. The vertical blue dottod lines are the true branching times. The vertical magenta dashed lines are the learned branching times.
Overall the model struggles with learning the true branching times and thus data point assignment suffers well: %MBCP = 79.4% vs %BGP = 79.8%, see Table 0 We hypothesise that this could be due to the carly branching output 0 being hard to model due to the branches staying close to each other. In particular, this means that output 0 contributes incorrect data point assignment suggestions to the other outputs thus making the model's task harder than it should be.
22

---

Published in Transactions on Machine Learning Research (05/2023)
Gene
2
Gene
0.0
Gene
2.0
Gene
1.5
1.0
1.0
1.5
0.5
2.0
10
0.5
-2
3.0
1.0
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
10
0.00
0.25
0.50
0.75
100
10
Gene 4
Gene 5
1.5
Gene
2.0
Gene 7
0.5
10
15
0.0
10
00
-0.5
-1.0
0.0
-1.5
).5
1.0
-2.0
1.0
.1.5
0.00
0.25
0.50
0.75
10
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
0.00
0.25
0.50
0.75
10
Gene
10
Gene
0.5
0.0
0000
0.5
0.0
0.5
1.00
0.00
0.25
0.50
0.75
1.00
0.25
0.50
0.75
10
Figure 12: This Figure shows the BGP model fit on sample 14 from Table 6 The data is represented in a scatter plot with input along the x axis and colours indicating branch assignments. The model's mean predictions are the solid lines. The vertical blue dashed lines are the true branching times. The vertical magenta dashed lines are the learned branching times. See Figure 11 for the key discussion points.
23

---

Published in Transactions on Machine Learning Research (05/2023)
All MMBGP incorrect labels
All BGP incorrect labels
300
250
200
150
100
50
0.0
0.2
0.4
0.6
0.8
10
0.0
0.2
0.4
0.6
0.8
10
Figure 13: Histogram of all incorrect data point assignments across all samples and all outputs. The z-axis is the input of the incorrectly labeled data, point. The y-axis is normalised to represent the incorrect data points in each bin percentage of all data points in that bin.
90
85
80
75
70
mBGP
BGP
Figure 14: Boxplots of data point assignment correctness averaged over all outputs for each of the 20 samples in MBGP vs BGP. Reported as a percentage, so that 100% is the best score you can hope for.
24

---

Published in Transactions on Machine Learning Research (05/2023)
(a) CTSG
(b) MPO
0.4
0.6
0.8
(c) PRTN3
(d) CALR
10.0
7.5
0.0
-2.5
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.5
0.8
1.0
0.5
0.0
0.2
0.1
0.6
0.8
1.0
0.0
3.
0.4
0.6
0.8
1.
(e) CAR1
CAR2
0.0
0.4
00
(g) GSTMI
(h) ELANE
1.0
0.0
0.4
1.0
Figure 15: Mouse haematopoictic stem cells (Paul et al.] 2015): Inconsistent cell assignment across the top eight branching genes. The 1000 randomly sub-sampled cells that have been used for inference coloured according to the posterior global branching pattern derived using the BGP algorithm. The points used to represent the cells have been assigned to the same branches, i.e. consistent across the eight genes. On the other hand, the triangles represent the inconsistent cell assignment, iee. cells that have been assigned to different branches across the eight genes.
25

---

Published in Transactions on Machine Learning Research (05/2023)
Table 6: This table contains the per-sample results for MBGP and BGP. Each row corresponds to a different 10 output noisy synthetic dataset. The RMSE columns mcasure the root mean square error between the learned and ground truth branching points. inenGP only applies to the BGP model and describes how many data points have inconsistent assignments across the outputs, that is, for how many data points do we assign them to one branch output A and to a different branch in output B? We note that lower is better and 100 is the absolute worst result you could get while 0 is the best result. % the percentage of correctly learned data point labels (higher is better).
this measure, we ignore inconsistent assignments and simply average over all outputs. Finally, 4x = %(M BGP) - (BGG)) above means MBGP is learning data point assignment better, below 0 mcans BGP is doing better.
As reported in the main body of the paper, MBGP outperforms BGP. In addition to that highlight that there's a correlation between how well MBGP learns the branching points and how well it assigns data point labels. When it learns the branching point locations better or equal to BGP, then it seems to consistently do better at data point assignment. As before, this motivates exploring schemes for learning more accurate branching locations.
RMSEMBGP
RMSE
iNCBGP
%(MBGP)
%(BGP)
A%
0
0.155656
0.161103
45
87.9
80.5
7.4
1
0.107528
0.172943
43
91.2
83.9
7.3
0.060403
0.248472
47
93.7
81.3
12.4
3
0.144460
0.209304
42
90.1
79.8
10.3
0.107687
0.122775
32
92.5
83.5
9.0
5
0.231522
0.171707
23
82.3
85.2
-2.9
6
0.145352
0.301008
43
87.7
75.0
12.7
7
0.125235
0.335639
34
89.2
77.8
11.4
8
0,277916
0.280737
27
80.4
78.1
2.3
9
0.149851
0,274390
37
89.7
76.9
12.8
10
0.249903
0.329334
41
81.2
69.4
11.8
0.239118 11
0.212270
25
79.8
81.2
-1.4
0.172359 12
0.123819
26
85.2
85.6
-044
13
0.149704
0.187880
13
87.1
85.0
2.1
14
0.329045
0.290695
20
79.4
79.8
-0.4
15
0.210808
0.230580
29
86.6
77.7
8.9
16
0.116043
0.114808
40
91.9
83.2
8.7
17
0.138261
0.296907
21
89.2
77.3
11.9
18
0.134857
0.165543
27
88.1
83.6
4.5
19
0.318348
0.218915
25
74.8
83.6
-8.8
Figure 16 shows an example of the BGP model applied on an early branching gene from mouse hematopoi- etic stem cell dataset (Paul et al.) 2015). The data contain 4423 cells. The log normalised data were transformed to wero mean for each individual gene. The DDRTree algorithm from Monocle 2 (Qiu et al. 2016) is used to infer the global cellular branching pattern as well as pseudotime for each cell. We have used a sub-sample of 870 cells and 30 inducing points to speed up the model inference. Figure 16 (a) shows the prior cell assignment of all 4423 cells which has been derived by running the Monocle 2 algorithm. This global branching structure used as the informative prior. Figure 16 (b) shows the posterior cell assignment of 870 cells that have been used in the inference (top sub-panel). The cells that are away from the global branching time (the black dashed line in Figure 16 (a)) have been assigned to either of the branches with high confidence. On the other hand, the cells closer to around the global branching point possess high level of uncertainty. This is also the case for the cells equidistant from both branches. The bottom sub-pancl of Figure 16 (b) shows the inferred gene-specific branching dynamics along with the posterior probability of the branching location. This uncertainty has been also reflected in Figure 16 (a) where the magenta background depicts the uncortainty associated with the geno-specific branching time (the blue solid line). The uncertain regions shown in both figures (magenta background in Figure 16 (a) and green dots near branching time in Figure
16 (b)) indicative of how cell assignment uncertainty is incorporated into the gene-specific
26

---

Published in Transactions on Machine Learning Research (05/2023)
branching-time posterior uncertainty. It emphasises one of the major benefits of developing probabilistic methods like the BGP for downstream analysis. As single-cell data are noisy, the cell assignment to different branches should be probabilistic while identifying the gene-specific branching dynamics.
(a)
(b)
Gene Expression
Gene Expression
00
15
Pseudotime
Pseudotime
Figure 16: Mouse haematopoietic stem cells (Paul al. 2015 BGP fit for the early branching gene MPO. The DDRIree algorithm of Monocle 2 is used estimate the pseudotime for each cell as well the global branching pattern. (a): Monocle global assignment of cells to the trunk state (purple) and two branches (green and red). The black dashed line is the global branching lime. The most probable gene branching time is shown by the blue solid line along with posterior uncertainty over the branching location (magenta background). (b): The top sub-panel shows the posterior cell assignment uncortainty for the cells used in the inference. Cell assignment to one of the branches is depicted. The bottom sub-panel represents the posterior probability distribution over the gene-specific branching location.
27


Reviews:
[{"review_id": "WochB0UQ32", "reviewer": "Reviewer_NqTz", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: \n1. the application of cell gene modeling is quite interesting \n2. utilization of multi-task \n3. clear improved performance\n\nWeakness:\n1. only compared against BGP, with many potentially simple multi-task baselines not considered. \n2. presentation needs improvements, in particular the early part of the paper. \n\n\nOther comments: \n-The introduction is not very clear in describing the problem, before discussing issues with current problems. What is the input/output of the problem, and what is exactly trajectories? Maybe an illustrative figure here would be helpful. In addition, discussion on why such a problem is important is needed. \n-  \" global branches\" what does this mean? if you have more genes, do they also belong to one of two branches?\n- notations are not consistent with each other, and the organization and presentaion is not quite clear: for example, just in problem definition\n1. both t and b are used to indicate time. \n2. bifurcation branch label should be clearly labeled with a letter and clarify it is an output\n3.  a property \"consistent cell labelling across genes\" shown in the problem definition before even defying the problem is strange. \n\nOther presentation issues: \n1. \"as we will see in §?? this mathematical convenience comes at a cost to modelling capabilities\".\n2. highlight the best performing methods in tables\n\n- \"are optimised along with other hyper-parameters\" generally only parameters are optimized. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "this paper proposes a Multivariate Branching Gaussian Process, an improvement over BGP, by considering multiple genres at the same time to achieve  label consistency and with better accuracy/efficiency shown. The main innovations are to use shared parameters (Z, H) and direct modeling of the branching time. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "\npresentation improvement and more empirical evaluations ", "broader_impact_concerns": "none", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "LgrlyNOUS5", "reviewer": "Reviewer_zLb8", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strong points: \n\n1) The proposed method considered joint inferences across all the genes.\n\n2) MBGP improved accuracy in branching time inference and assignment accuracy, and it is faster than the original BGP.\n\n## Weak points:\n\n1) The empirical experiments and comparisons are not sufficient;\n\n2) The existing work on which this study is based does not seem too new or current;\n\n3) Too many typos and format issues.\n\nIn addition, although there are empirical experiments in this paper, I did not find codes related to these experiments, so I am not sure the empirical description/experiments would be enough to reproduce.\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Modeling sequential branching dynamics is an interesting topic, especially, it is important for analyzing single-cell gene expression data. In this paper, the author(s) tried to fix the limitations in BGP model and proposed a multivariate BGP version, called MBGP, to jointly analyze branch assignment inference. And through the experiments on synthetic data and a real-world single-cell RNA-seq data, the author(s) demonstrated the advantages of the proposed method.\n\n\n", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Detailedly, I have the following comments/questions:\n\n1. By experiments with only two datasets (that is, one simulated dat set and one real dataset), then it concluded that the proposed method performed better than other existing method(s), I think it is not convincing. Moreover, the author(s) did not provide a comparison of the proposed method with existing method(s) on that real dataset in the main text.\n\n2. In the introduction, the author(s) mentioned several existing methods, why not compare the proposed method with these methods? In experiments, the author(s) should at lease explain the reason.\n\n3.The existing work on which this study is based does not seem too new or current. \n\nAfter looking at the references of this paper, it seems that only two citations are from 2020 onwards, namely 2020 and 2021. Since there is not a lot of overlap between my research and the work of this paper, I'm not sure whether the author(s) haven't done enough search on the literature, or this field has stagnated. Could the author(s) please explain that?\n\n4. Many issues/typos:\n\n(1) On page 3, Line 5, \"??\", what does it denote?\n\n(2) In the main text, sometimes mathematical formulas have punctuation at the end and sometimes they don't.\n\n(3) Giving the full name the first time \"BGP\" appears in the main text, not until section 3.\n\n(4) The caption of Table 1 should end with punctuation.\n\n(5) Please unify the table styles in the main text.\n\n(6) In Figure 5, the numbers on the axes are too small.\n\n(7) In Figure 5, what does MMBGP mean?\n\n(8) Many formatting inconsistencies in the references, please check and correct them carefully.\n\n\nI look forward to the response from the author(s). Thanks!\n\n\n\n\n", "broader_impact_concerns": "No.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "frLQebyt6Q", "reviewer": "Reviewer_TGmh", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Significance.** The proposed method (combination of model and inference procedure) seems very bespoke and application-specific. I am unsure if the framework would readily transfer beyond the analysis of gene expression data. Now, within this particular application, I am not well-versed in the domain, so I cannot assess with any confidence the importance or difficulty of solving these particular subproblems.\n\n**Novelty / originality.** The originality of the contribution is nontrivial but modest. Augmenting a model and its inference by going beyond the mean-field assumption is a conceptually simple improvement. Proposing a computationally feasible inference procedure for this model doesn’t appear to be too challenging since the necessary components are already in place. However, combining them into one coherent framework turns out as usual to be not so trivial and this work seems to have done that successfully, as evidenced by their qualitative and quantitative experimental results.\n\n**Clarity / Presentation.** The structure of this paper is well-organized and the writing is mostly clear and effective.\n\nSection 2 defines provides a formal definition of the two interdependent problems that this paper seeks to solve. Section 3 formalizes the existing BGP model and relates it to OMGP. Thereafter, they illustrate the failure modes of this class of models in solving the second subproblem of interest — namely, the cell labelling problem. In brief, it is possible for the same cell to be assigned to varying branches across different genes, which is biologically impossible. Section 4 summarizes their core contribution and Section 5 provides empirical results to validate their method.\n\nIn general, the presentation could be adapted more to the ML research community that primarily constitutes the target audience of TMLR. For example, I think there is room for improvement in the problem definition section, which is perhaps a bit too terse and assumes a certain degree of familiarity with gene expression data. Some other issues are also referred to in the Requested Changes section.\n\n**Quality / soundness.** I checked the equations and derivations, particularly those that describe the proposed, and these seem sound and correct. Some equations seem unnecessarily verbose, and could probably be deferred to the appendix. For example eq. 11 explicitly spells out a multivariate normal log density, eq. 12 is quite fleshed out. In terms of the proposed inference procedure, the authors propose a two-stage approach to estimating the GP hyperparameters and the branching points and cell assignments. While this is somewhat justified by intuition, it would be more rigorous to include an ablation study that compares against a joint treatment of these (hyper)parameters. Otherwise, the experimental setup appears to be sound and the empirical improvements over BGP appear to be significant.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The Overlapping Mixture of Gaussian Processes (OMGP) and its variants like the Branching Gaussian Process (BGP) are models designed to detect bifurcation points in single-cell gene expression data. However, the current BGP model has a limitation in that it independently infers the assignment of observations to latent functions for each output dimension (gene), which can lead to inconsistent assignments and reduced accuracy in branching time inference. In this paper, the authors propose a new model called the multivariate branching Gaussian process (MBGP) and perform joint branch assignment inference across multiple output dimensions to ensure assignment consistency and improve accuracy in branching time inference and assignment accuracy. The authors test the effectiveness of their so-called “multivariate branching Gaussian process” model on synthetic data and a real dataset of mouse stem cells and develop a variational inference procedure to scale the model to large datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Broadly speaking, I find the figures in this to be quite difficult to parse. Beyond simple things like adding axis labels (which are missing in almost all figures), there is generally some room for improvement in things like figure captions and legends and beyond.\n\nFor instance, it isn't sufficient to label subfigures \"MPO\" and \"CTSG\" (in Figure 2) even if it is stated explicitly somewhere in the text that these are genes, since this likely won’t immediately occur to the target audience of TMLR.\n\nFigure 2 is tough to parse. Not only are there many different components (markers with different colours and shapes, colour bars, vertical shaded regions, curves, etc.) there are no labels in the form of figure legends, let alone axis labels. Furthermore, it is not possible to understand the figure from the caption alone, one really has to parse through the main text.\n\n### Misc\n\n- On Page 3, a section reference did not render correctly\n- Page 9: \"Finally, obtain a sparse collapsed variational bound by [...]\"", "broader_impact_concerns": "No Broader Impact Statement section is included. However, in my opinion, there are no obvious concerns regarding the ethical implications of the work that necessitates the inclusion of this section.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The claims of the paper relate to an improvement of the BGP model and application to single-cell transcription experiments. This claim seems to be sufficiently interesting to the audience to warrant consideration for this venue. Some reviewers commented on lack of comparison to other competing methods and presentation issues that would improve the paper substantially. The authors responded to those comments and offered to update the paper with four specific changes including a more current review of literature. I encourage to make the changes that they suggested in response to the comments.
Decision: Accept (poster)