Scaling up Continuous-Time Markov Chains Helps Resolve Underspecification
Authors: Alkis Gotovos
Rebekka Burkholz
John Quackenbush
Stefanie Jegelka
Venue: NeurIPS 2021 Poster (2021)
Keywords: discrete probabilistic models
continuous-time Markov chains
underspecification
approximate inference
cancer genomics

Abstract: Modeling the time evolution of discrete sets of items (e.g., genetic mutations) is a fundamental problem in many biomedical applications. We approach this problem through the lens of continuous-time Markov chains, and show that the resulting learning task is generally underspecified in the usual setting of cross-sectional data. We explore a perhaps surprising remedy: including a number of additional independent items can help determine time order, and hence resolve underspecification. This is in sharp contrast to the common practice of limiting the analysis to a small subset of relevant items, which is followed largely due to poor scaling of existing methods. To put our theoretical insight into practice, we develop an approximate likelihood maximization method for learning continuous-time Markov chains, which can scale to hundreds of items and is orders of magnitude faster than previous methods. We demonstrate the effectiveness of our approach on synthetic and real cancer data.

Scaling up Continuous-Time Markov Chains Helps Resolve Underspecification
Alkis Gotovos MIT alkisg@mit. edu
Rebekka Burkholz
Harvard University
rburkholz@hsph. harvard. edu
John Quackenbush
Stefanie Jegelka
Harvard University johng@hsph. harvard. edu
MIT stefje@mit.e
Abstract
Modeling the time evolution of discrete sets of items (c.g., genetic mutations) a fundamental problem in many biomedical applications. We approach this problem through the lens of continuous-time Markov chains, and show that the resulting learning task is generally underspecified in the usual setting of cross-sectional data. We explore a perhaps surprising remedy: including a number of additional independent items can help determine time order, and hence resolve underspecifi- cation. This is in sharp contrast to the common practice of limiting the analysis to small subset of relevant items, which is followed largely due to poor scaling of existing methods. To put our theoretical insight into practice, we develop an ap- proximate likelihood maximization method for learning continuous- chains, which can scale to hundreds of items and is orders of magnitude faster than previous methods. We demonstrate the effectiveness our approach on synthetic and real cancer data.
Markov
Introduction
Modeling the time evolution of physical processes with discrete states is an important machine learn- ing problem that spans wide range of application domains, including molecular dynamics (Crom- melin & Vanden-Eijnden, 2006), phylogenetics (Suchard et al., 2001), and computational medicine (Liu et al., 2015). Continuous-time markov chains have found success in such modeling tasks, sup- ported by extensive research into learning and inference in such models (Opper & Sanguinetti, 2008; Perkins, 2009; Archambeau & Opper, 2011; Rao Tch, 2012).
While the majority of previous work has focused on modeling a relatively small number of discrete states, many interesting applications involve the interaction between sets of items, which results in exponentially large state spaces. For example, when modeling the time progression of cancer, we want consider number of genetic alterations (e.g., mutations, copy number variations, etc.) that exhibit complex time-related dependencies (Beerenwinkel et al., 2014). To model how sets of n such alterations evolve over time, we need to consider state space of size 2n. What makes the problem even more challenging is the fact that the available data sets are typically cross-sectional; that they consist of (unordered) sets of items observed unknown time points without any further information about the history of the underlying process.
In this paper, focus the problem of learning particular parametric family of continuous-tima Markov chains from such cross-sectional data, and show that the resulting problem is general underspecifted. The issue of underspecification is faceeed and has been shown to permeate
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

large number of machine learning systems, often leading to poor generalization, lack of robustness, and spurious relationships when interpreting the resulting models (D' Amour et al., 2020). In our setting, we explore perhaps surprising remedy: including number of additional (approximately) independent items can help determine the time order of process events, and hence resolve underspec- ification. We theoretically show that these extra items act as "background clock", since counting the number of their occurrences in a data sample can help us estimate the time at which that sample was observed. In sharp contrast to common practice, which limits the analysis to small subset of relevant, highly-interacting items, our insight suggests that scaling up the learning procedure can be erucial for the robustness of the inferred models. We thus show that items deemed priori unimpor- tant to the application at hand may in fact be particularly valuable for recovering the time properties of the underlying physical process.
Existing learning approaches are, unfortunately, not well-suited for practically applying this insight; for example, the lfte-offtheearr method by Schill et al. (2019) for learning continuous-tin Markov chains to model cancer progression, scales exponentially in the number of items considered, thus limiting the analysis to around 20 items. To alleviate this issue, we propose an approximate likeli- hood maximization method that relies on a fast gradient approximation. On problems of 20 items our approach runs almost 1000 times faster that the state of the art, while it can also scale to prob- lems involving hundreds of items. In experiments on real cancer data, we demonstrate how some previous results may have been artifacts of underspecification, and how scaling up the analysis can result in more robust models of cancer progression.
2 Problem setup
Given a ground set V fl....,,n,, we consider a continuous-time Markov chain (CTMC) (22) on state space 2 (Grimmett & Stirzaker, 2001). Thus, the states of the Markov chain are subsets of V. and can equivalently be identified as binary vectors of size n. CTMCs are com- monly represented by a generator matrix Q E R="x2". Since the rows and columns of this matrix correspond to states in 2", we will use qs-R to denote the entry of Q at row indexed by S and column indexed by R. For S + R. the entry q--R represents the infinitesimal rate of transitioning from state S to state R, that is,
P(X,+8
R X = S) for $ R
Q--R = lim
81-+++
86
with gs- n 2 0. The quantity Øs == -S- represents the total infinitesimal rate of leaving state S, with Øs 2 O. Since the total rate of leaving S is equal to the sum of the rates of transitioning to any other state R. we have ûs = ERey Qs-R- for all S E 2V.
We are interested in modeling physical processes that induce an increasing accumulation of items over time, for example, the accumulation of genetic alterations when modeling cancer progression. To this end, we impose the following two constraints on the CTMC. First, we assume that the chain always starts from the empty state, that is, X() = V. Second, we only allow transitions that add a single item to the current set, that is, for S o R. qs-R # 0 only if R SU {?) for some V. Since items are never removed, this implies that lim X, V.
In this paper, we focus learning the CTMC from a given data set. Since learning directly the exponentially sized generator Q is a hopeless endeavor, we impose a particular parametric form on this matrix, recently introduced by Schill et al. (2019). We define a parameter matrix O pax= (also denoted by 0 as vector), and parameterize the off-diagonal entries of the generator matrix as
(s-su(n exp Ozj
Σ 0is YES
for all S C V, j E V\.S. The off-diagonal elements of θ encode positive (attractive) negative (repulsive) pairwise interactions between items in V. The presence of item increases (6iz > 0) or decreases (0., < 0) the rate of adding item j by a multiplicative factor of Wij := c.u,, The diagonal elements 0jj encode the intrinsic rate of adding element j when ignoring all other interactions.
Obtaining a sample from a CTMC under all aforementioned constraints can be done using the fol- lowing sequential procedure. At each step, given the current state S, we draw "holding time" h

---

Algorithm 1: Sampling a set from the marginal CTMC
Input: Parameters 0
Draw laha Exp(1)
t + 0, S 0
while < Lax do
Draw ~ (0)
Compute p. = su() st
Draw Cat(V \ (pi)ae)s))
1 ( - h. S S (x)
return S
from an exponential distribution with parameter Ús- This represents how much time will pass until we add the next item. Then, we draw the next item E V \ S according categorical distribution with probabilities x (S-Su(=)- The result a sequence of items (o)....,on) or E V, together with sequence of holding times (hi,..., f.,), h, R.
Ideally, our data set would consist of a number of such item sequences and holding times; we would then proceed learn the parameters 0 by maximizing the likelihood of our CTMC model. Unfortunately, such detailed data is typically not available in practice. In particular, it is rather common in biomedical applications only have access to cross-sectional data, i.e., to only observe the current state (in the form of an unordered subset of V) at a particular point in time without knowing how the process reached that state. Furthermore, is often the case that the time point of each observation can only be roughly estimated or, worse, is completely unknown.
Similarly to previous work, we assume that we are given data D { ). C V. with each 51()) representing state observed at potentially different unknown time Tats More specifically, we assume that Tats is random variable following an exponential distribution with density ptt) = e (Gerstung et al., 2009; Schill al., 2019), and investigate the problem of maxi- mizing the marginal likelihood,
p(D; 0)
II (SS)
II P(S()) ( (10)((t)d
Algorithm shows how to sample a set from p(S:0). The difference from our previous description of sampling from CTMC lies in the fact that the procedure only run until we reach the observation time tobs, and the result is an unordered set without any information about the holding times.
3 The value of "unimportant" items
Our ultimate goal is retrieve the ordered interaction structure present in the data by learning the parameter vector 9. Before attempting to maximize the marginal likelihood of eq. (1), though, natural first ask whether the limited information available in the data is sufficient infer such ordered interactions in the first place. We start with a simple example to showcase this issue.
3.1 Warmup example
Consider ground set V = (1, and true model parameters 0f = 03 O. and 0i2 = - for some a > O. When a is sufficiently large, the probability of item occurring before item goes to zero, therefore the model encodes the fact that item 1 is prerequisite (and therefore, always appears before) item 2. As a consequence, we will also have p([2); 0*) s 0 for the resulting marginal distribution of sets. One may intuitively think that this could be enough information to infer the true time order; the following proposition shows that this is not the case.
Proposition There is a one-dimensional family of parameters f ((s), and 81, 82 E R, such that p(:; (()) p(:: 0²), ). for all
< 82-
The proof follows by simple algebra, and can be found Appendix together with an illustration of the parameter family. Interestingly, for + $1, the model (s)) encodes the fact that, when

---

observe both items, item 2 always appears before item 1, which is the opposite of what is encoded in our true model.
3.2 Independent items background clock
While the above example paints a pessimistic picture for inferring the correct time order, we show that this is still possible given some additional side information. Suppose that we given an- other ground set V+ whose items have no interaction with the items in V, that is, the parameter matrix rul of the resulting model over Viu V U V+ has block diagonal structure with blocks θ R|VXXVI and O+ R Then, it is easy to that the distributions over V and V+ are conditionally independent given time t. That is, for any S C V, S+ S V+, we can write U(S t; Oun) p(5 t; 0) pS + t; 0+). Note that the same statement does not hold when con- sidering the marginal distributions, i.e., PSSUS+: P(un) & p(S: 0) p(S+: ++;+). Using this conditional independence property, we can rewrite the marginal probability of S U S+ as follows:
00
p(S U S4; (
p(S U US
t:: O dt))dt
p(s 0) p(5 t1: (p)(()d
(by cond. ind.)
p(s 1; 0) p() S+ 0 (p( )+ (by Bayes' rule)
= S4: 0mn)
p(S t: 0) (t) S+ )dt
(dividing with p(S+;0+)) (2)
Comparing equations (1) and (2), we see that the information about S "ginned by observing S4 can be explained via a posterior observation time distribution p(t S4; 04) that refines the prior p(t) by taking into account S4.
To gain further insight into this time posterior, we analyze in more detail a simplified setup, in which V+ consists of m identically distributed independent items parameterized by 0+. that is, O+ 0 In this case, |S+1| It follows a binomial distribution over m variables with success parameter 1 --e-wtt where W4 := e++ Intuitively, we expect these independent items to act as a "background clock": counting the number of items | should give us an indication about the observation time t. Moreover, as we increase the number of items m we should be able to get increasingly accurate estimates of the observation time. The following theorem formalizes these ideas.
Theorem 1. Let S C V be randomly drawn according the CTMC with parameter matrix O +. and let to be the true observation time of S4. Then, for any 8 E (0, 1), there exists an mo, such that for all m Z ma, the mean and variance of the posterior time distribution p(t S= 0) can be bounded as follows with probability at least 8,
log
log m m
Mppost t- tC((0,,t")
+ O
m
1 m
1 m
Vpost S O2(0+,t")
The quantities C₁ and C2 are constant w.r.t. m and encode how suitable a particular rate (quantified through 0_)) of an item is for estimating the specific observation time Intuitively, larger rates are better suited for estimating smaller times and vice versa. We further illustrate this effect in Appendix together with the detailed proof of the theorem. Note that for m + oo we get p(t S+; -> 6(t t²) in distribution, and therefore P(S S4: P + p(S ()) This means that observing an arbitrarily large amount of independent items equivalent to knowing the true observation time. In Appendix C we show that similar asymptotic convergence behavior can be observed even when the items in V+ are not identically distributed or exhibit some dependencies.
3.3 Discussion
To illustrate how this result applies to our two-item example discussed before, we define O Im (0 ..... of o(m))TT where each ol0 drawn uniformly in [-4, -2]. This corresponds individual

---

2
4
(a) m =
(b) T
4
5
2
25
(d) 100
Figure 1: Posterior time densities for N = 500 samples. Posteriors of samples with S = {1} are colored red, while those of samples with S {1, 2} are colored blue. The dashed line highlights the time posterior of specific sample, while the black square denotes its true observation time.
frequencies of items in S+ that roughly range from 0.01 to 0.1. We then draw N = 500 samples from the CTMC defined by Oun Figure shows the posterior time density corresponding each sample with S = {1} (red) and S = 2} (blue), for different numbers m of independent items. do not show the densities for samples with S 0. As expected, when m all samples follow the prior time density, and we have no way to distinguish their time order. However, we keep increasing m we notice that the two colors keep separating from each other in time, and it becomes clear that the red posteriors tend to concentrate around earlier observation times than the blue ones. This is a clear indication that item comes before item 2, and the additional time information obtained by the posterior time estimates should help pinpoint the true model within the family discussed in Proposition 1.
To avoid clutter,
Due to computational considerations, it is common practice to constrain the ground set to a small subset of items that interact strongly and deemed interesting the application at hand; for ample, genetic alterations that are known to have role in cancer development (Raphael & Vandin, 2015; Schill et al., 2019). Our results this section challenge this practice. We show that items be- lieved to be independent of the "interesting items", and thus considered unimportant to the analysis, can in fact be particularly valuable in estimating the observation time of data samples and recovering the time order properties of the underlying physical process.
For the remainder of this paper, will not assume any prior knowledge about the block-diagonal structure of @ull or the values the parameters 0.. Furthermore, we will make no distinction between "items of interest" and independent items. If there is indeed block-diagonal parameter structure in the data, the addition of regularization term that promotes sparsity can help recover this structure the presence of noise. For this reason, similarly to what was used by Schill et al. (2019), the final objective we maximize in practice is the -erruullrrzed marginal log-likelihood,
N
F(D:0)
1 N
Σ
log log Σ 10
(3)
d=1
5

---

Efficient approximate likelihood maximization
To put the previously discussed theoretical insight into practice, we need to be able to efficiently max- imize the above objective for ground sets V containing potentially hundreds of items. When using a first-order method for optimization, the bulk of the required computation is devoted to obtaining the gradient of the marginal log-likelihood p(S; 0) with respect to the parameters O. Computing the exact gradient requires ((2") computation; we propose here a method to compute a gradient approximation in a much more efficient manner. We start by deriving expressions for the likelihood in some simplified setups, and then use these results as building blocks for our setup of interest.
4.1 Full sequences
First, we consider setup in which we observe the arrival order of the items in the CTMC, which corresponds to building up an ordered sequence instead a set S in Algorithm 1. Furthermore, we assume that lah- +oo, which means that the returned sequence will be permutation of the ground set V. By definition of the CTMC, the probability of observing a sequence o (o1...., on) E Sv, where Sy denotes the permutation group over V, can be written
4a 1- (0)
p(a 0)
II P PX. X = =10=0(-1)) Ob-1) = II
(4)
go (
1=1
We denote subsequences using the bracket notation, [ = (oI...., ;;), with d 0. Also, with a slight abuse notation we directly use sequences in some places that require sets; therefore, the notation drr is to be understood as Bel(()) where set(ar := {a; i (1,...,/ol)} |ol)].
4.2 Partial sequences given time
As next step, instead of tubs +00, we assume that we are given the observation time tobs, and would like to derive the probability of observing a partial sequence time tobs- We can express this as the probability of the intersection of two events in the space of outcomes of the continuous-time Markov chain. Event A identifies the outcomes in which the first ordered elements agree with our partial sequence a. Event B identifies the outcomes for which the observation time labs falls between the addition of the k-tth and k + 1-th elements; that is, Tr < lobs < T We can then express the desired probability as p(a Lobs; 0) P (A) IP (B A).
(dI..... Tk), with h S n,
The probability of event A can be directly written using (4), except that, in this case the product taken from 1 k, instead of Regarding event B, first note that the sequence of events 1=1 )m>
((T: Lobs
AA)
is decreasing. As result, can write
( A) = P (Tk < tobs < Tx+1 A) = P (T& < tabs A) P (T&+1 < fats A
(5)
By definition of the CTMC, each holding time H₃ is an exponentially distributed random variable whose parameter depends on the i - items that have been added up to that point. We can also define the jump time Tę Σ H,, which represents the time at which the k-th change of state occurs. When we condition on knowing the first & items of the sequence, i.e., event A, the jump times Tp and T distributed sums of k and + independent exponential random variables respectively. The CDF of the sum of r independent exponential variables with rates Aj,-.., "ay...;Ay, Ap, also known as the hypoexponential distribution (Bibinger, 2013), is given by
r
1
FHEXP ,I....,>pp)
|I A
Σ
A (1
In our case, the rates A; of each exponential distribution are q (0), thus we compute the terms of (5) as IP (B) A) /HHxp(fatsi ...,, (()))))) FHHxP 9),1). ......... (ou))))
4.3 Marginal partial sequences
Next we consider the probability of observing a partial sequence a (ai,..., RR), with & S n, without knowledge of the observation time tobs We define events A and B as before, except that in

---

this case the observation time Tots is a random variable, and event B is defined by Ta < Tobs < Tk+1- In particular, we assume that Ta is an exponential random variable with rate Xats = 1. Rather than naively integrating over time the probability derived before, we can obtain the following greatly simplifed expression by making use of the memoryless property of exponential random variables (Bertsekas & Tsitsiklis, 2008) when computing P (B A).
Proposition 2. The marginal probability of a purtial sequence (0i,.... Th) can be written as
(1 (0) II L0 + 80 n-1 (
p(a; 0)
+ 1er (0)
The proof can be found in Appendix D, but the form of this expression allows for an intuitive in- terpretation. We can modify the original continuous-time Markov chain by adding an extra dummy state Xohs, such that for all X C V. -x (0) = 1, and ax implies that there
O. The first property
fixed rate to transition from any state to ohs: while the second property im- plies that Xahs is terminal state. It is easy to see then, that eq. (6) expresses exactly the probability of observing a "full" sequence a in this modified chain. Note that a "full" sequence in this chain is no longer permutation of V, but rather any sequence that ends with Xobs.
4.4 Marginal sets
Finally, we consider the probability of observing a S C V. without knowledge of the order in which the individual items arrived, and still under the assumption that the observation time Tats is an exponential random variable with rate Xobs 1. Making use of the result in the previous section, we can compute this probability by summing over all partial sequences that are permutations of S, that is, p(S; 0) Eoess p(o;; 0). However, this computation is infeasible for all but very small set sizes, since it requires summing over |S|! terms. To alleviate this problem, we propose a method to approximate the gradient of log p( S; 0). This gradient can be written as
Vo log ((: = 0) Vep(S:
V (: 0) p(S:0)
1
Σ oESs
Vop(o:0)
p(a: 0) aESs p(S:0)
Σ
Vo logp(o;; 0).
In the last expression, we observe that the gradient corresponding to each permutation a is weighed by the probability of observing that permutation of the given set S. This suggests stochastic approximation of the desired gradient by first sampling m permutations g(1), (MM) according to (((1):0)= := :p::0)/::00)00)) O):p(S;(pp 0), and then computing the average of the resulting gradients.
M
Ve logp(S;0))
Vo (o
M
Conceptually similar approaches have been used to approximate the gradient of the log normalizeer when maximizing the likelihood of energy-based probabilistic models (Song & Kingma, 2021).
remains to show how to obtain samples from p(- SS 0), which we do by using a Markov chain Monte Carlo method on Sg. More concretely, we employ a Metropolis-Hastingg chain (Levin et al., 2009); at each time step, given the current permutation a,, the chain proposes a new permutation Tnew according proposal distribution Q (), and transitions Onew with probability
min (1, pl new S: 5:0) Q( Onew: 0)
Paccopt
p S 0) Q (Tnew 0;
A simple choice for Q is the uniform distribution over all permutations regardless of the current state. We instead employ a more sophisticated proposal that based on the form of the marginal probability (6), and leads to faster and more stable learning in practice. While detailed discussion and theoretical analysis of mixing is beyond the scope of this paper, in Appendix we present our proposal in detail, and compare it the uniform.
5 Further related work
There has been a long line of research focused on understanding and reconstructing the evolution- ary history of tumors (Nik-Zainal et al., 2012; Welch et al., 2012; Turajlic et al., 2015; Mitchell
7

---

2
0.4
0.4
KL. divergence
3
5
4
0.2
0.2
10
20
10
30
Figure 2: The KI. divergence of recovered vs. true model as function of the number of extra in- dependent items, evaluated on two synthetic examples with known true parameters. The dashed lines represent the learned models given the true observations times. Error bars are over 96 repetitions.
al., 2018; Jolly& & Van Loo, 2018). However, these papers do not attempt to systematically model the timed interaction of genetic alterations across data set. More recently, Gerstung et al. (2020) showed that modern data resulting from deep sequencing procedures may contain partial order in- formation about specific genetic events. Incorporating this additional information into analysis like ours, could potentially help with further reducing underspecification. Finally, there is evidence that specific classes of mutational processes generate mutations at constant rate, and are largely unaffected by other factors (Alexandrov al., 2015; Campbell al., 2020; Alexandrov et al., 2020). The resulting so-called mutational signatures have been likened to biological clock, and could be candidates for explicitly defining the set V+ discussed earlier (cf. Section 3.2).
Experiments
The following experimental setup is common all our experiments. To optimize the objective (3), we use proximal AdaGrad method (Duchi et al., 201 We fix the initial step size to 1) 1, and the regularization weight to À 0.01. To initialize the parameter matrix O we train a diagonal model for 50 epochs, and then draw each off-diagonal entry from Unif([-0.2, 0.2]). For the gradient approximation (7), we use M 50 samples addition 10 burn-in samples that are discarded. The code used to run our experiments can be found at https: //tpthugbbttt.t//// cwm coe//rcoooooo /time/ tree/clean.
6.1 Synthetic data
We start by evaluating our approach on two synthetic data sets drawn from known true models. The first data set identical to the one discussed in Figure 1, and contains observations of two items that tend to follow a specific time order (item comes before item 2). The second data set contains observations of five items following the graphical structure shown on the upper right of Figure Edges with arrows +- represent positive interactions (0,, > 0), while edges with vertical lines + represent negative interactions (B., < 0). This type of structure, containing groups of items with bidirectional negative interactions, is common in genomics data; in fact, this exact structure can be found in the model learned from real cancer data in our following experiment (see Appendix F).
Our goal is to evaluate how well we can recover the true model by optimizing the regularized marginal likelihood; in particular, we want examine the effect of adding number m of inde- pendent items to the data. To quantify the quality of recovery, we compare the true distribution Price of marginal sequences, computed using eq. (6), to the distribution of marginal sequences p" of the learned model, approximated by drawing 10t marginal sequences from that model. In Figure we depict the KL divergence (pL(p" PPrue) for increasing numbers of extra items. To compare against the ideal case, we also train a model by maximizing the regularized marginal likelihood given the true observation times (cf. Section 4.2). The KL divergence between this model and the true one is shown as dashed line in the figure. We repeated each experiment 96 times; sources of randomness include the choice of independent items, parameter initialization, and the gradient approximation. We see that the recovery quality improves for both data sets with increasing m, and rapidly approaches the ideal. This shows that observing large enough number of additional items is practically equivalent to knowing the true observation times, and thus validates Theorem 1.

---

Table Learning run times (Intel Core i9 CPU)
Method
10
15
20
50
100
Schill et al. (2019)
2s
43 s
121 m
Ours
3
4
8s
1m 5s
33 43
EGFR(A)
PDGFRA(A)
first
EGFR
PDGFRA
first
first
20
40
60
80
100
20
40
60
80
100
TP53
MDM2(A)
first
first
IDH1
CDK4/A) first
first
20
40
60
80
100
0
20
40
60
80
100
11
11
Results on the TCGA glioblastoma data set indicating the tendencies the inferred time order of four selected pairs of genetic alterations. Error bars over 24 repetitions.
Figure 3:
6.2 Real cancer data
We next evaluate our approach on a data set that contains multiforme, aggressive type of brain cancer. The data is part of the TCGA PanCancer Atlas project (https:/ / //t///ww..r//// cancer wwwcccr.goo//ccoe gov/ttga), and we obtained a preprocessed version via cBioPortal (Cerami et al., 2012). We followed some further filtering procedures from previous work (Leiserson et al., 2013) and ended up with ground set n 410 items, each of which represents a point mutation, amplification, such genetic alterations is of fundamental importance to cancer research, as can help illuminate the processes that are involved in cancer initiation and progression.
378 tumor samples of glioblastoma
deletion of specific gene. Analyzing the interaction structure between
For the remainder of this section, we assume that the items in V are ordered by decreasing frequency, and when selecting subset we keep the most frequent items. In Table 1, we show that our learning approach runs about 1000 times faster on ground set of size n 20 compared to the exponentially scaling exact approach proposed by Schill et al. (2019).
In the absence of ground truth, we use the following procedure evaluate how some of the learned parameters behave as we increase the size of the ground set. Given pair of items b E V, and n-2 other (most frequent) items, we first learn CTMC model on the n total items using our approach. Next approximate the sequence distribution of just and b (marginally over time and over all other items) by sampling 10th marginal sequences from the learned model. Finally, we compute the proportion of sequences in which occurred before 6, given that both where observed; when this proportion is greater than 0.5, we infer that a is more likely to occur before b, and vice versa. In Figure we present these computed proportions for four chosen pairs of genetic alterations that are interesting both from a biological standpoint, as well as in terms of learning behavior.
First, we discover two significant ordered interactions, shown in the top row of the figure, which were not reported previous work analyzing glioblastoma data (Raphael & Vandin, 2015; Cristea et al., 2017: Schill et al., 2019). Our results suggest that amplification events for genes EGFR and PDGFRA tend to occur before their respective point mutation events. Interestingly, the co- occurrence of amplifications and mutations for EGFR has been observed before (Leiserson et al.,

---

2013: Sanchez- Vega et al., 2018), but there is little known about their time order. Note that the time order the PDGFRA alterations can be robustly inferred only after including more than 70 additional items in the analysis.
Second, confirm some previously reported interactions, for example, the tendency of IDH1 muta- tions to occur before TP53 ones, as shown in the bottom left of the figure. This was observed before by Schill al. (2019), and is also supported by some biological evidence (Watanabe al., 2009). Again, we can see that our result becomes robust only after n = 50. This indicates that the previous observation of this time order may have been in part due setup, e.g.. fixed initialization (see Appendix G). The final result in the bottom right of the figure seems support this claim: the amplifications of MDM2 and CDK4 do seem have clear time order, even though previous work has inferred both that CDK4(A) occurs before MDM2(A) (Cristea al., 2017); as well the opposite order (Schill et al., 2019).
fortuitous choice of the optimization
Finally, the parameter matrix θ learned from this data set (see Appendix G) has indeed an approxi- mately block-diagonal structure, which confirms the practical relevance of Theorem
7 Conclusion
Our experiments demonstrated that, in the presence of underspecification, special care is required when interpreting results that may be significantly affected by particular choices the optimization setup. In this paper, proposed effective way to mitigate such effects, namely scaling up the problem and taking into account additional items that are commonly discarded as insignificant. We hope that our approach will inspire further work into combating underspecifieation and learning robust time evolution models.
Acknowledgments and Disclosure of Funding
This work was partially supported by a Swiss National Science Foundation Early Postdoc. Postdoc.Mobility fellowship, by NSF CAREER award 1553284, NSF BIGDATA IIS-1741341, and NSF award 1717610, as well as a grant from the US National Cancer Institute (IR35CA220523).
References
Alexandrov, L. B., Jones, P. H., Wedge, D. C., Sale, E., Campbell, P. J., Nik-Zainal, S., and Stratton, M. R. Clock-like mutational processes in human somatic cells. Nature Genetics, 2015.
Alexandrov, B., Kim, J., Haradhvala, N. J., Huang, M. N. Tian Ng, A. W., Wu. Y., Boot, A.. Covington, K. R., Gordenin, D. A., Bergstrom, E. N., Islam, S. M. A., Lopez-Bigas, N., Klimczak, L. J., McPherson, J. R., Morganella, S., Sabarinathan, R., Wheeler, D. A., Mustonen, V., Getz, G., Rozen, S. G., and Stratton, M. R. The repertoire of mutational signatures in human cancer. Nature, 2020.
Archambeau, C. and Opper, M. Approximate inference for continuous-time markov processes. In Barber, D., Cemgil, A. T., and Chiappa, S. (eds.), Bayesian Time Series Models. 2011.
Beerenwinkel, N., Schwarz, R. F., Gerstung, M., and Markowetz, F. Cancer evolution: Mathematical models and computational inference. Systematic Biology, 2014.
Bertsekas, D. and Tsitsiklis, N. Introduction probability. Athena Scientific, 2008
Bibinger, M. Notes on the sum and maximum of independent exponentially distributed random variables with different scale parameters. arXiv, art. 1307.3945 2013.
Campbell, P. J., Getz, G., Korbel, J. O., Stuart, J. M., Jennings, J. L., Stein, L. D., Perry, M. D., Nahal-Bose, H. K., Ouellette, B. F. F., Li, C. H., Rheinbay, E., Nielsen, G. P., Sgroi, D. C., Wu, C.-L., Faquin, W. C., Deshpande, V., Boutros, P. C., Lazar, A. J., Hoadley, K. A., Louis, D. N., Dursi, L. J., Yung, C. K., Bailey, M. H., Saksena, G., Raine, K. M., Buchhalter, I., Kleinheinz, K., Schlesner, M., Zhang, J., Wang, W. Wheeler, D. A., Ding, I.., Simpson, J. T. OConnor, B. D., Yakneen, S., Ellrott, K. Miyoshi, N., Butler, A. P., Royo, R., Shorser, S. I., Vazquez, M.,
10

---

Rausch, T., Tiao, G., Waszak. S. M., Rodriguez-Martin, B., Shringarpure, S., Wu, D.-Y., Demidov,
G. M., Delaneau, O., Hayashi, S., Imoto, S., Habermann, N., Segre, A. V., Garrison, E.. Cafferkey. A., Alvarez, H. G., Hereeia- J. M., Muyas, F., Drechsel, O., Bruzos, A. L.,, Temes, J., Zamora, J., Baez- Baez-Ortega, A., Kim, H.-I.., Mashl, R. J., Ye, K., DiBiase, A., Huang, K.-I., Letunic, I. McIellan, M. D., Newhouse, S. J., Shmaya, T., Kumar, S., Wedge, D. C., Wright, M. H., Yellapantula, V. D., Gerstein, M., Khurana, E., Marques-Bonet, T., Navarro, A., Bustamante, C. D., Siebert, R., Nakagawa, H., Easton, D. F., Ossowski, S., Tubio, J. M. C., De La Vega, F. M., Estivill, X., Yuen, D., Mihaiescu, G. L., Omberg, L., Ferretti, V., Sabarinathan, R., Pich, O., Gonzalez-Perez, A., Taylor-Weinese A., Fittall, M. W. Demeulemeester.. J., Tarabichi, M., Roberts, N. D., Van Loo, P., Cortés-Ciriano, I., Urban, L., Park, P., Zhu, B., Pitkänen, E., Li, Y., Saini, N., Klimezak, L. J.. Weischenfeldt, J., Sidiropoulos, N., Alexandrov, L. B., Rabionet, R., Escaramis, G., Bosio, M., Holik, A. Z., Susak, H., Prasad, A., Erkek, S., Calabrese, C., Raeder, B., Harrington, F.., Mayes, S., Turner, D., Juul, S.,, Roberts, S. A., Song, I.., Koster, R., Mirabello, L.,, Hua, X., Tanskanen, T. J., Tojo, M., Chen, J., Aaltonen, L. A., Rätsch, G., Schwarz, R. F., Butte, A. J., Brazma, A., Chanock, S. J., Chatterjee, N., Stegle, O., Harismendy, O., Bova, G. S., Gordenin, D. A., Haan, D., Sieverling, L., Feuerbach, L., Chalmers, D., Joly, Y., Knoppers, B., Molnár-Gábor, F., Phillips, M., Thorogood, A., Townend, D., Goldman, M., Fonseca, N. A.. Xiang, Q. Craft, B., Piñieiro Yáñez, E., Muñoz, A., Petryszak, R., Füllgrabe, A., Al-Shahrour, F,, Keays, M., Haussler, D., Weinstein, J., Huber, W. Valencia, A.. Papatheodorou. I., Zhu, J., Fan, Y., Torrents. D.. Bieg, M., Chen, K., Chong, Z., Cibulskis, K., Eils. R., Fulton, R. S., Gelpi. J. L., Gonzalez, S., Gut, I. G., Hach, F., Heinold, M.. Hu, T., Huang, V., Hutter, B., Jäger, N., Jung. J., Kumar, Y., Lalansingh, C., Leshchiner, I., Livitz, D., Ma, E. Z.,, Maruvka, Y. E., Milovanovic, A., Nielsen, M. M., Paramasivam, N., Pedersen, S., Puiggrös, M., Sahinalp, S. C., Sarrafi, I., Stewart, C., Stobbe, M. D., Wala, J. A., Wang. J., Wendl, M., Werner, J., Wu, Z., Xue, H., Yamaguchi, T. N., Yellapantula, V., Davis-Dusenbery,
N., Grossman, R. L.,, Kim, Y., Heinold, M. C., Hinton, J., Jones, D. R., Menzies, A., Stebbings. L., Hess, J. M., Rosenberg, M.. Dunford, A. J., Gupta, M., Imielinski, M., Meyerson, M., Beroukhim, R., Reimand, J., Dhingra, P., Favero, F. Dentro, S., Wintersinger, J., Rudneva, V., Park, J. W., Hong, E. P., Heo, S. G., Kahles, A., Lehmann, K.-V. Soulette, C. M., Shiraishi, Y.. Liu, F., He, Y., Demirciolu, D., Davidson, N. R., Greger. L., Li, S., Liu, D., Stark, S. G., Zhang. F., Amin, S. B., Bailey, P., Chateigner, A., Frenkel- Morgenstem, M., Hou, Y., Huska, M. R., Kilpinen. H., Lamaze, F. C., Li, C., Li, X., Li, X., Liu, X., Marin, M. G., Markowski, J., Nandi, T., Ojesina, A. I., Pan Hammarrmrr,,r, Q., Park, P. J., Pedamallu, C. S., Su, H., Tan, P., Teh, B. T., Wang. J., Xiong, H., Ye, C., Yung, C., Zhang, X., Zheng, L., Zhu, S., Awadalla, P., Creighton, J., Wu, K., Yang, H., Göke, J., Zhang. Z., Brooks, A. N., Fittall, M. W., Martincorena, I., Rubio-Perez, C., Juul, M., Schumacher, S., Shapira, O., Tamborero, D., Mularoni, L., Hornshoj, H., Deu-Pons, J., Muifios, F., Bertl, J., Guo, Q., Gonzalez- Perez, A., Xiang, Q., and The ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium. Pan-cancer analysis of whole genomes. Nature, 2020.
Cerami, E., Gao, J., Dogrusoz, U., Gross, B. E., Sumer, S. O., Aksoy, B. A., Jacobsen, A., Byrne, C. J., Heuer, M. L.., Larsson, E., Antipin, Y., Reva, B., Goldberg, A. P., Sander, C., and Schultz, N. The cBio cancer genomics portal: An open platform for exploring multidimensional cancer genomics data. Cancer Discovery, 2012.
Cristea, S., Kuipers, J., and Beerenwinkel, N. pathTiMEx: Joint inference of mutually exclusive cancer pathways and their progression dynamics. Journal of Computational Biology, 2017.
Crommelin, D. T. and Vanden-Fijnden, Fitting timeseries by continuous-time markov chains: quadratic programming approach. Journal of Computational Physics, 2006.
D' D'Amour, A., Heller, K., Moldovan, D., Adlam, B., Alipanahi, B., Beutel, A., Chen, C., Deaton, J., Eisenstein, J., Hoffman, M. D., Hormozdiari, F., Houlsby, N., Hou, S., Jerfel, G., Karthike- salingam, A., Lucic, M., Ma, Y., McLean, C., Mincu, D., Mitani, A., Montanari, A., Nado, Z., Natarajan, V., Nielson, C., Osborne, T. F., Raman, R., Ramasamy, K., Sayres, R., Schrouff, J., Seneviratne, M., Sequeira. S., Suresh, H., Veitch, V., Vladymyrov, M., Wang, X., Webster, K., Yadlowsky, S., Yun, T., Zhai, X., and Sculley, D. Underspecification presents challenges for credibility in modern machine learning. arXiv, art. 201 2011.03395, 2020.
Duchi, J., Hazan, E., and Singer, Y. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research (JMLR), 2011.
11

---

Gerstung, M., Baudis, M., Moch, H., and Beerenwinkel, N. Quantifying cancer progression with conjunctive bayesian networks. Bioinformatics, 2009.
Gerstung, M., Jolly, C., Leshchiner, I., Dentro, S. C., Gonzalez, S., Rosebrock, D., Mitchell, J., Rubanova, Y., Anur, P., Yu, K., Tarabichi, M., Deshwar, A., Wintersinger, J., Kleinheinz, K., Vázquez-García, I., Haase, K., Jerman, L., Sengupta, S., Macintyre, G., Malikic, S., Donmez, N. Livitz, D. G., Cmero, M., Demeulemeester, J., Schumacher, S., Fan, Y., Yao, X., Lee, J., Schlesner, M., Boutros, C., Bowtell, D., Zhu, H., Getz, G., Imielinski, M., Beroukhim, R., Sahinalp, S. C.,, Ji, Y., Peifer, M., Markowetz, F., Mustonen, V., Yuan, K., Wang, W., Morris, Q. D., Spellman, P. T., Wedge, D. C., and Loo, P. V. The cvolutionary history of 2,658 cancers. Nature, 2020.
Grimmett, G. and Stirzaker, D. Probability and random processes. Oxford University Press, 2001.
Jolly, C. and Van Loo, P. Timing somatic events in the evolution of cancer. Genome Biology, 2018.
Leiserson, M. D. M., Blokh, D., Sharan, R., and Raphael, B. J. Simultancous identification of multiple driver pathways in cancer. PLoS Computational Biology, 2013.
Levin, D. A., Peres, Y., and Wilmer, E, L. Markov chains and mixing times. American Mathematical Society, 2009.
Liu, Y.Y.,, Li, S., Li, F. Song, L., and Rehg, J. M. Efficient learning of continuous-time hidden markov models for disease progression. In Neural Information Processing Systems (NeurIPS), 2015.
Mengersen, K. L. and Tweedie, R. L. Rates of convergence of the Hastings and Metropolis algo- rithms. The Annals of Statistics, 1996.
Mitchell, T. J., Turajlic, S., Rowan, A., Nicol, D., Farmery, J. H. R., O'Brien, T., Martincorena, I., Tarpey, P., Angelopoulos, N., Yates, L. R., Butler, A. P., Raine, K., Stewart, G. D., Challacombe, B., Fernando, A., Lopez, I., Hazell, S., Chandra, A., Chowdhury, S., Rudman, S., Soultati, A., Stamp, G., Fotiadis, N., Pickering. L., Au, L., Spain, I.., Lynch, J., Stares, M., Teague, J., Maura, F. Wedge, D. C., Horswell, S., Chambers, T., Litchfield, K., Xu, H., Stewart, A., Elaidi, R., Oudard, S., McGranahan, N., Csabai, I., Gore, M., Futreal, P. A., Larkin, J., Lynch, A. G., Szallasi, Z., Swanton, C., Campbell, P. J., and TRACERx Renal Consortium. Timing the landmark in the evolution of clear cell renal cell cancer: TRACERx renal. Cell, 2018.
Nik-Zainal, S., VanãLoo, P., Wedge, C., Alexandrov, L. B., Greenman, C. D., Lau, K. W., Raine, K. Jones, D., Marshall, J., Ramakrishna, M., Shlien, A., Cooke, S. L., Hinton, J., Menzies, A., Stebbings, L. A., Leroy, C., Jia, M., Rance, R., Mudie, L. J., Gamble, S. J., Stephens, McLaren, S., Tarpey, P. S., Papaemmanuil, E., Davies, H. R., Varela, I., McBride, D. J., Bignell, G. R., Leung, K., Butler, A. P., Teague, J. W., Martin, S., Jönsson, G., Mariani, O., Boyault, S., Miron, P., Fatima, A., Langerød, A., Aparicio. S. A. J. R., Tutt, A., Sieuwerts, A. M.. Borg, k., Thomas, G., Salomon, A. V., Richardson, A. L., Børresen- Dale, A.-L., Futreal, A., Stratton, M. R., and Campbell, P. J. The life history of 21 breast cancers. Cell, 2012.
J..
Opper, M. and Sanguinetti, G. Variational inference for markov jump processes. In Neural Informa- tion Processing Systems (NeurIPS), 2008
Perkins, T. Maximum likelihood trajectories for continuous-time markov chains. In Neural Infor- mation Processing Systems (NeurIPS), 2009.
Rao, V. and Teh, Y. Memc for continuous-time discrete-state systems. In Neural Information Processing Systems (NeurIPS), 2012.
Raphael, B. J. and Vandin, F. Simultaneous inference of cancer pathways and tumor progression from cross-sectional mutation data. Journal of Compulational Biology, 2015.
Sanchez- Vega,, F., Mina, M., and more. Oncogenic signaling pathways in the cancer genome atlas. Cell, 2018.
Schill, R., Solbrig, S., Wettig, T., and Spang, R. Modelling cancer progression using mutual hazard networks. Bioinformatics, 2019.
12

---

Song, Y. and Kingma, D. P. How to train your energy energy-based based models. arXiv, art. 2101.03288, 2021.
Suchard, M. A., Weiss, R. E., and Sinsheimer, J. S. Bayesian selection of continuous-time markov chain evolutionary models. Molecular Biology and Evolution, 2001.
Turajlic, S., McGranahan, N., and Swanton, Inferring mutational timing and reconstructing tumour evolutionary histories. Biochimica et Biophysica Acta (BBA) Reviews Cancer, 2015.
Watanabe, T., Nobusawa, S., Kleihues, P., and Ohgaki, H. IDHI mutations early events in the development of astrocytomas and oligodendroolionas. The American Journal of Pathology, 2009.
Welch, J. S., Ley, T. J., Link, D. C., Miller, C. A., Larson, D. E., Koboldt, D. C., Wartman, L. D., Lamprecht, T. L., Liu, Xia, J., Kandoth, C., Fulton, R. McLellan, M. D.. Dooling, D. J., Wallis, J. W., Chen, K., Harris. C. C.. Schmidt, H. K., Kalicki- Veizer, J. M., Lu, C., Zhang, Q., Lin, I.., O'Lauuglin, M. D., McMichael, J. F., Delehaunty, K. D., Fulton, A., Magrini, J., McGrath, S. D., Demeter, R. T., Vickery, T. L.., Hundal, J., Cook, L. I.., Swift, G. W. Reed, P., Alldredge, P. A., Wylie, T.. N., Walker, R., Watson, M. A., Heath, E., Shannon, W. D., Varghese, N., Nagarajan, R., Payton, J. E., Baty, J. D., Kulkarni, S., Klco. J. M., Tomasson, M. H., Westervelt, P., Walter, M. J., Graubert, T. A., DiPersio, J. F., Ding, L., Mardis, R., and Wilson, R. K. The origin and evolution of mutations in acute myeloid leukemia. Cell, 2012.
13


Reviews:
[{"review_id": "kFv6GzdRgbm", "reviewer": "Reviewer_ZrQu", "summary": "This paper deals with learning a particular class of multi-component continuous time Markov chains (CTMC) when time stamps are unobserved in the data. In this class binary components change their state irreversibly, rates follow a pairwise interaction model and one component changes at a time. The authors show that modeling independent elements provides information about the posterior distribution of observation time and provide an efficient stochastic approximation to the likelihood’s gradient.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "**Originality**: \nInferring temporal relations between sequences using independent mutations and a CTMC model is already a fundamental technique in phylogenetics. Therefore, I suggest mentioning the relation of the proposed approach to this field. Yet, as far as I know, the explicit quantification of the added value of independent components (Theorem 1) is novel.\nThe approximation technique is a novel combination of formulas for inferring the likelihood of partial sequences, stochastic gradient approximation and MCMC sampling.\n\n**Quality**: \nThe added value of additional independent components is supported by a sound derivation in Equation 2. This added value is quantified in Theorem 1 which is proved in the appendix.\nThe approximation algorithm contains a series of well-established steps.\nThe experimental section gives a good illustration of actual gain of leveraging additional independent components and illustrate a very nice application of modeling cancer tumor progression.\nWhile the title describes the method as a general CTMC, the actual model is a specific type that assume irreversible changes of components with pairwise influences on the rates. This subclass has an interesting important application and therefore it does not diminish the significance of this work. However, this specificity should be more clearly stated, at least in the abstract.\n\n**Clarity**: \nThe paper is clearly written. Here are some suggestions to improve readability:\n1.\t‘Underspecified’ is a broad term and makes the abstract a bit vague. It may be good to state the sentence from the introduction -- ‘consist of unordered sets of item observed at unknown time point’ -- already in the abstract.\n2.\tThe cancer example is very illustrative and provides a strong motivation. Furthermore, the paragraph in line 162 motivates the usage of independent mutations. It may be good to move it to the introduction, to make the motivation clearer at this stage.\nMinor questions:\n3.\tTheorem 1: should V be replaced by V_{+}?\n4.\tFigure 1: do the different curves represent posteriors given different values of S_{+}? \n5.\tWhere the N tumor samples taken from N distinct patients, or were multiple samples taken from some patients? This may influence the choice of modeling approach.\n\n**Significance**: \nThis work may be an important step in enabling modelling of cancer tumor progression, and may have also applications in other domains in which systems accumulate irreversible changes over time with dependencies among some of the components.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I do not see a potential negative social impact that is specific to this work.\n\nThe paper lacks a discussion of limitations of the method. For example, the experimental section used a fixed value of lambda. This seems to have a major effect on the sparsity of the inferred interaction matrix. Is there a reasonable way to tune this hyperparameter as done in prediction tasks? \n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6 hours"}, "extra_text": {}}, {"review_id": "YCs23nxBXWn", "reviewer": "Reviewer_SqMU", "summary": "This paper describes a continuous-time Markov Chain method for modeling the time evolution of events. By including independent items the authors are able to increase the limit on the number of items modeled. This research is particularly relevant as larger single-cell sequencing data sets are being collected.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "  The underlying biomedical problem is well-motivated and the authors briefly mention other possible applications. I think they could have provided more clarity on the specific problem --- the use of the term \"item\" is confusing in paragraph 2. The reader wonders if an \"item\" is a cell or a genetic locus, though this is clarified later in the paper.\n\n  Going from cross-section data to longitudinal inferences requires assumptions about the distribution of the samples. It would be helpful to be more clear about the assumptions that enable the underlying inferential goal of this work.\n\n  The \"warm-up\" example clearly illustrates an important point. But the title \"warm-up\" seems pedantic. Perhaps a more descriptive title such as \"A counter-example for inferring time order\" would help the reader to understand the motivation of the example. The discussion in section 3.3 is helpful in that it links back to the counterexample.\n\n  There are some standards for genetic nomenclature and it does not seem the case that the authors have used that nomenclature here. This leads to some confusion as to what exactly is meant by \"amplification\" in figure 3. Is this an increase in copy number or increase in transcription or translation due to an epigenetic change or something else. It would be helpful to clearly state what the data can and can't inform about the mechanisms.\n\n  The meaning of the error bars in figure 3 is not clearly stated.\n\n  The development of the method is predicated on a set of \"interesting items\" $V$ and a set of \"independent items\" that are used to help inform the clock (lines 1-168). Then, in the paragraph on lines 169-174, the authors state, \"we will make no distinction between 'items of interest' and independent items\". It is not clear, theoretically, what motivates or justifies this shift in perspective. Prior to this change, I was left considering how one might select \"interesting\" and \"independent\" genetic loci, but it seems that has been rendered unnecessary. It would be helpful for the authors to clarify what motivates and theoretically justifies the elimination of the need to specify $V_+$ - the \"independent\" set of items.\n\n  The description of the experiments makes it impossible to reproduce the results of the paper. On line 280, the authors state that they evaluate on a data set of tumor samples from glioblastoma from TCGA, but TCGA contains expression, mutation, and methylation data. It is not clear what data the authors used and the url to the entire project does not enable reproduction of the analysis. Perhaps a more detailed description with code can be provided in supplementary information.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "  The authors only compared to Schill (2019). There are many tumor evolution method available in the bioinformatics literature and without a more comprehensive comparison, the utility of this approach may not be clear and adoption of the approach may be limited.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "FO-nWN2pCJ0", "reviewer": "Reviewer_YnoU", "summary": "The authors model the time evolution of discrete states (e.g. genetic mutations) using a continuous-time Markov chains. They show that the resulting learning task is generally underspecified given cross-sectional data. The authors suggest including additional independent items that can help determine time order, and resolve this underspecification.  The authors implement an approximate likelihood maximization method for learning continuous-time Markov chains, which shows good performance on high-dimensional data, is faster than competing methods, and shows promising results on synthetic and real cancer data.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "Originality: The work is a novel combination of well-known techniques. Related work needs to be adequetely cited\n\nQuality: The submission is technically sound justified by theoretical analysis and experimental results. This is a complete piece of work and is built heavily using Schill 2019. \n\nClarity: The submission is clearly written. \n\nSignificance: The study and results are important in many fields. There is potential for others to use the ideas. The work provides a unique theoretical approach.\n\n* Comments to authors\n\n   * The problem setup section as well as Section 5.1 heavily lean on Schill 2019: https://academic.oup.com/bioinformatics/article/36/1/241/5524604\nThis fact is not clearly cited in the paper or Problem section.\nIt seems that the current paper builds on Schill 2019, while also uses the same cancer data for real-world experiments, with the value-add of theoretically proving the usage of additional data. \n\n   * Regarding the additional unimportant data, how much and how complete should this be? The Markov chain would require that these additional data are present at (t-1), so that the chain can learn at t. What if additional data is missing at (t-1) or is available only intermittently? Would the equations in Section 4.2 still hold?\n\n  * Is the maximum dim of \\Theta_{full}, sum of dim of V and V_{+}? (Section 3.2). So in cases, when the dim of \\Theta_{full} is < sum of dim of V and V_{+}, shouldn’t the off-diagonal information be taken into account as well?\n  \n  * From Theorem 1, Appendices C and G, it seems that \\Theta_{full} has to maintain a block-diagonal structure. By adding additional unimportant data, across time points, there is a higher likelihood to see dependency across time points (i.e important signals in the off-diagonal elements). This is generally the case when identifying trajectories of cell differentiation in single-cell data across multiple (pseudo) time points and adding incomplete information for certain time points when available. It is a non-standard and underspecified problem, in which case, the dependency of additional terms is what helps connect time points. Based on this, aren’t the authors losing important information by enforcing a block-diagonal structure? How easy is it to incorporate information from off-diagonal elements into the CTMC?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations were not mentioned in the paper. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4 hours"}, "extra_text": {}}, {"review_id": "3xj3ZN8oeaQ", "reviewer": "Reviewer_HeAE", "summary": "The authors tackle the problem of time evolution of discrete sets through the view of continuous-time markov chains. The authors provide a theoretical justification of why these problems in general underspecified in the setting of cross-sectional data. They go on to demonstrate that the underspecification can be overcome by use of additional items / features that can help determine time order by deriving a bound on the mean and variance of the observed time for these additional items. An approximate formulation for likelihood maximization is then presented without any assumptions of prior knowledge of the additional items or the structure of the additional items. The theoretical framework is then tested using simulation studies and demonstrate that accuracy of state sequence recovery can be significantly improved with use of the additional items. The framework is also evaluated using a real world cancer dataset from TCGA by ordering of pairs of mutation and copy number aberration events. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "This is a very interesting submission that demonstrates how underspecification in continuous-time markov chains can be resolved by use of additional independent features/ which might not directly relevant to the problem under consideration. The paper is well written and the use of a consistent example throughout the manuscript conveys the concepts clearly. My comments are largely related in practical applications:\n1. The simulation results presented by the authors are convincing and a clear demonstration of the concept. However it is not clear to me how well the formulation would hold up under noisy data - perhaps larger markov chains where some observations might be missing subset of states would provide a practical view of the limitations \n2. Any individual cancer type has been demonstrated to contain distinct sub types each with associated mutations and/or copy number alterations. For eg, the IDH mutation in GBM highlighted by the authors is observed in only a subset of individuals and is associated with better survival. Therefore a broad span of patient data will most likely include mutation /copy number aberration drivers that are distinct and mutually exclusive to each other. How is the learning process influenced in the presence of such confounding factors and how does it affect the interpretation of the results. \n\t\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors address the limitations of the behavior of the algorithm when assumption of independence of additional features is violated.\n\nThe work does not appear to lead to any potential negative societal impacts in my opinion. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "12"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper describes a method for modeling the time-series changes in a discrete set of items. The particular application the authors have in mind is phylogenetic inference where the items are genetic mutations. The method involves reframing and using continuous-time Markov chains. While moving to the continuous domain helps address the combinatorial nature of the problem, it introduces another problem - underspecification - which results in poor generalization, sensitivity, and false positive correlations. The authors propose a solution based on prior work (Schill et. al 2019) that employs so-called "unimportant items" to help recover the time order.

The reviewers noted two specific aspects of significance of the work. First, the paper provides bounds on the expected amount of information gain from including these "unimportant items". Second, the paper addresses an important problem in bioinformatics. Inferring the tumor evolution process is a well-studies problem in genomics and has been studied extensively. In clinical applications, the tumor evolution may be of limited concern because the physician is primarily interested in treating the whole tumor as it presents (including any heterogeneity). But, in fundamental research, understanding the tumor evolution can inform our basic understanding of the processes that lead to tumor growth and metastatic processes and therefore phylogenetic inference is a significant problem to address. It should be noted that for this method to be used in practice, a distinction must be made between "interesting" and "uninteresting" items. This may be done by the model or with prior knowledge. The authors have constructed their algorithm to allow the model to make this distinction which may improve the generality for problems outside of phylogenetic inference.

I'm in agreement with the reviewers after considering the author response and discussion that this paper should be recommended for acceptance. The theoretical contributions may prove useful for a broad range of similar problems and the practical utility on a data set relevant for practical phylogenetic inference is demonstrated.

Decision: Accept (Poster)