summaries: 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.  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. 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. 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.
abstractText: 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.
authors.name: Alkis Gotovos
authors.name: Rebekka Burkholz
references.authors: L.B. Alexandrov P.H. Jones D.C. Wedge J.E. Sale P.J. Campbell S. Nik-Zainal M.R. Stratton
references.title: Clock-like mutational processes in human somatic cells
references.venue: Nature Genetics,
references.authors: L.B. Alexandrov J. Kim N.J. Haradhvala M.N. Huang A.W. Tian Ng Y. Wu A. Boot K.R. Covington D.A. Gordenin E.N. Bergstrom S.M.A. Islam N. Lopez-Bigas L.J. Klimczak J.R. McPherson S. Morganella R. Sabarinathan D.A. Wheeler V. Mustonen G. Getz S.G. Rozen M.R. Stratton
references.title: The repertoire of mutational signatures in human cancer
references.venue: Nature,
references.authors: C. Archambeau M. Opper
references.title: Approximate inference for continuous-time markov processes
references.venue: Bayesian Time Series Models
references.authors: N. Beerenwinkel R.F. Schwarz M. Gerstung F. Markowetz
references.title: Cancer evolution: Mathematical models and computational inference
references.venue: Systematic Biology,
references.authors: M. Bibinger
references.title: Notes on the sum and maximum of independent exponentially distributed random variables with different scale parameters
references.venue: arXiv, art
references.authors: E. Cerami J. Gao U. Dogrusoz B.E. Gross S.O. Sumer B.A. Aksoy A. Jacobsen C.J. Byrne M.L. Heuer E. Larsson Y. Antipin B. Reva A.P. Goldberg C. Sander N. Schultz
references.title: The cBio cancer genomics portal: An open platform for exploring multidimensional cancer genomics data
references.venue: Cancer Discovery,
references.authors: S. Cristea J. Kuipers N. Beerenwinkel
references.title: pathTiMEx: Joint inference of mutually exclusive cancer pathways and their progression dynamics
references.venue: Journal of Computational Biology,
references.authors: D.T. Crommelin E. Vanden-Eijnden
references.title: Fitting timeseries by continuous-time markov chains: A quadratic programming approach
references.venue: Journal of Computational Physics,
references.authors: S. Yadlowsky T. Yun X. Zhai D. Sculley
references.title: Underspecification presents challenges for credibility in modern machine learning
references.authors: J. Duchi E. Hazan Y. Singer
references.title: Adaptive subgradient methods for online learning and stochastic optimization
references.venue: Journal of Machine Learning Research (JMLR),
references.authors: M. Gerstung M. Baudis H. Moch N. Beerenwinkel
references.title: Quantifying cancer progression with conjunctive bayesian networks
references.authors: S.C. Sahinalp Y. Ji M. Peifer F. Markowetz V. Mustonen K. Yuan W. Wang Q.D. Morris P.T. Spellman D.C. Wedge P.V. Loo
references.title: The evolutionary history of 2,658 cancers
references.venue: Nature,
references.authors: G. Grimmett D. Stirzaker
references.title: Probability and random processes
references.authors: C. Jolly P. Van Loo
references.title: Timing somatic events in the evolution of cancer
references.venue: Genome Biology,
references.authors: M.D.M. Leiserson D. Blokh R. Sharan B.J. Raphael
references.title: Simultaneous identification of multiple driver pathways in cancer
references.venue: PLoS Computational Biology,
references.authors: Liu Y.-Y S. Li F. Li L. Song J.M. Rehg
references.title: Efficient learning of continuous-time hidden markov models for disease progression
references.venue: In Neural Information Processing Systems (NeurIPS),
references.authors: K.L. Mengersen R.L. Tweedie
references.title: Rates of convergence of the Hastings and Metropolis algorithms
references.venue: The Annals of Statistics,
references.authors: P. Miron A. Fatima A. Langerød Aparicio S.A.J. R A. Tutt A.M. Sieuwerts Borg G. Thomas A.V. Salomon A.L. Richardson Børresen-Dale A.-L P.A. Futreal M.R. Stratton P.J. Campbell
references.title: The life history of 21 breast cancers
references.venue: Cell,
references.authors: M. Opper G. Sanguinetti
references.title: Variational inference for markov jump processes
references.venue: In Neural Information Processing Systems (NeurIPS),
references.authors: T. Perkins
references.title: Maximum likelihood trajectories for continuous-time markov chains
references.venue: In Neural Information Processing Systems (NeurIPS),
references.authors: V. Rao Y. Teh
references.title: Mcmc for continuous-time discrete-state systems
references.venue: In Neural Information Processing Systems (NeurIPS),
references.authors: B.J. Raphael F. Vandin
references.title: Simultaneous inference of cancer pathways and tumor progression from cross-sectional mutation data
references.venue: Journal of Computational Biology,
references.authors: F. Sanchez-Vega M. Mina more
references.title: Oncogenic signaling pathways in the cancer genome
references.venue: atlas. Cell,
references.authors: R. Schill S. Solbrig T. Wettig R. Spang
references.title: Modelling cancer progression using mutual hazard
references.authors: Y. Song D.P. Kingma
references.title: How to train your energy-based models
references.venue: arXiv, art
references.authors: M.A. Suchard R.E. Weiss J.S. Sinsheimer
references.title: Bayesian selection of continuous-time markov chain evolutionary models
references.venue: Molecular Biology and Evolution,
references.authors: S. Turajlic N. McGranahan C. Swanton
references.title: Inferring mutational timing and reconstructing tumour evolutionary histories
references.venue: Biochimica et Biophysica Acta (BBA) - Reviews on Cancer,
references.authors: T. Watanabe S. Nobusawa P. Kleihues H. Ohgaki
references.title: IDH1 mutations are early events in the development of astrocytomas and oligodendrogliomas
references.venue: The American Journal of Pathology,
references.authors: J. P P.A. Alldredge T.N. Wylie J.R. Walker M.A. Watson S.E. Heath W.D. Shannon N. Varghese R. Nagarajan J.E. Payton J.D. Baty S. Kulkarni J.M. Klco M.H. Tomasson P. Westervelt M.J. Walter T.A. Graubert J.F. DiPersio L. Ding E.R. Mardis R.K. Wilson
references.title: The origin and evolution of mutations in acute myeloid leukemia
references.venue: Cell,
sections.heading: 1 Introduction
sections.text: Modeling the time evolution of physical processes with discrete states is an important machine learning problem that spans a wide range of application domains, including molecular dynamics (Crommelin & 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, supported by extensive research into learning and inference in such models (Opper & Sanguinetti, 2008; Perkins, 2009; Archambeau & Opper, 2011; Rao & Teh, 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 to consider a 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 a 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 is, they consist of (unordered) sets of items observed at unknown time points without any further information about the history of the underlying process.
  In this paper, we focus on the problem of learning a particular parametric family of continuous-time Markov chains from such cross-sectional data, and show that the resulting problem is in general underspecified. The issue of underspecification is many-faceted and has been shown to permeate a
  35th Conference on Neural Information Processing Systems (NeurIPS 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 a perhaps surprising remedy: including a number of additional (approximately) independent items can help determine the time order of process events, and hence resolve underspecification. We theoretically show that these extra items act as a “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 a small subset of relevant, highly-interacting items, our insight suggests that scaling up the learning procedure can be crucial for the robustness of the inferred models. We thus show that items deemed a priori unimportant 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 state-of-the-art method by Schill et al. (2019) for learning continuous-time 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 likelihood 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 problems 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.
sections.heading: 2 Problem setup
sections.text: Given a ground set V = {1, . . . , n}, we consider a continuous-time Markov chain (CTMC) {X(t)}t≥0 on state space 2V (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 commonly represented by a generator matrix Q ∈ R2n×2n . Since the rows and columns of this matrix correspond to states in 2V , 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 qS R represents the infinitesimal rate of transitioning from state S to state R, that is,
  qS R = lim δt→0+ P(Xt+δt = R |Xt = S) δt , for S ̸= R,
  with qS R ≥ 0. The quantity q̃S := −qS S represents the total infinitesimal rate of leaving state S, with q̃S ≥ 0. Since the total rate of leaving S is equal to the sum of the rates of transitioning to any other state R, we have q̃S = ∑ R∈2V qS R, for all S ∈ 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(0) = ∅. Second, we only allow transitions that add a single item to the current set, that is, for S ̸= R, qS R ̸= 0 only if R = S ∪ {i} for some i ∈ V . Since items are never removed, this implies that limt→∞Xt = V .
  In this paper, we focus on 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 Θ ∈ Rn×n (also denoted by θ as a vector), and parameterize the off-diagonal entries of the generator matrix as
  qS S∪{j}(θ) = exp ( θjj +
  ∑ i∈S θij
  ) ,
  for all S ⊆ V , and j ∈ V \S. The off-diagonal elements of Θ encode positive (attractive) or negative (repulsive) pairwise interactions between items in V . The presence of item i increases (θij > 0) or decreases (θij < 0) the rate of adding item j by a multiplicative factor of wij := eθij . The diagonal elements θjj 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 following sequential procedure. At each step, given the current state S, we draw a “holding time” h
  Algorithm 1: Sampling a set from the marginal CTMC Input :Parameters θ Draw tobs ∼ Exp(1) t← 0, S ← ∅ while t < tobs do
  Draw h ∼ Exp(q̃S(θ)) Compute pi = qS S∪{i}(θ)/q̃S(θ), for i ∈ V \ S Draw x ∼ Cat(V \ S, (pi)i∈V \S) t← t+ h, S ← S ∪ {x}
  return S
  from an exponential distribution with parameter q̃S . This represents how much time will pass until we add the next item. Then, we draw the next item x ∈ V \S according to a categorical distribution with probabilities ∝ qS S∪{x}. The result is a sequence of items (σ1, . . . , σn), σi ∈ V , together with a sequence of holding times (h1, . . . , hn), hi ∈ R. Ideally, our data set would consist of a number of such item sequences and holding times; we would then proceed to learn the parameters θ 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 to 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, it 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 a data setD = {S(1), . . . , S(N)}, S(i) ⊆ V , with each S(i) representing a state observed at a potentially different unknown time Tobs. More specifically, we assume that Tobs is a random variable following an exponential distribution with density p(t) = e−t (Gerstung et al., 2009; Schill et al., 2019), and investigate the problem of maximizing the marginal likelihood,
  p(D;θ) = N∏ i=1 p(S(i);θ) = N∏ i=1 ∫ ∞ 0 p(S(i) | t;θ)p(t)dt. (1)
  Algorithm 1 shows how to sample a set from p(S;θ). The difference from our previous description of sampling from a CTMC lies in the fact that the procedure is only run until we reach the observation time tobs, and the result is an unordered set without any information about the holding times.
sections.heading: 3 The value of “unimportant” items
sections.text: Our ultimate goal is to retrieve the ordered interaction structure present in the data by learning the parameter vector θ. Before attempting to maximize the marginal likelihood of eq. (1), though, it is natural to first ask whether the limited information available in the data is sufficient to infer such ordered interactions in the first place. We start with a simple example to showcase this issue.
sections.heading: 3.1 Warmup example
sections.text: Consider a ground set V = {1, 2} and true model parameters θ∗11 = θ∗21 = 0, and θ∗12 = −θ∗22 = α, for some α > 0. When α is sufficiently large, the probability of item 2 occurring before item 1 goes to zero, therefore the model encodes the fact that item 1 is a prerequisite (and therefore, always appears before) item 2. As a consequence, we will also have p({2};θ∗) ≈ 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 1. There is a one-dimensional family of parameters θ = θ(s), and s1, s2 ∈ R, such that p(· ;θ(s)) = p(· ;θ∗), for all s1 < s < s2.
  The proof follows by simple algebra, and can be found in Appendix A together with an illustration of the parameter family. Interestingly, for s → s1, the model θ(s) encodes the fact that, when we
  observe both items, item 2 always appears before item 1, which is the opposite of what is encoded in our true model.
sections.heading: 3.2 Independent items as a background clock
sections.text: 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 are given another ground set V+ whose items have no interaction with the items in V , that is, the parameter matrix Θfull of the resulting model over Vfull = V ∪ V+ has a block diagonal structure with blocks Θ ∈ R|V |×|V | and Θ+ ∈ R|V+|×|V+|. Then, it is easy to see that the distributions over V and V+ are conditionally independent given time t. That is, for any S ⊆ V , S+ ⊆ V+, we can write p(S ∪S+ | t;θfull) = p(S | t;θ) p(S+ | t;θ+). Note that the same statement does not hold when considering the marginal distributions, i.e., p(S∪S+;θfull) ̸= p(S;θ) p(S+;θ+). Using this conditional independence property, we can rewrite the marginal probability of S ∪ S+ as follows:
  p(S ∪ S+;θfull) = ∫ ∞ 0 p(S ∪ S+ | t;θfull)p(t)dt
  = ∫ ∞ 0 p(S | t;θ) p(S+ | t;θ+)p(t)dt (by cond. ind.)
  = ∫ ∞ 0 p(S | t;θ) p(t |S+;θ+)p(S+;θ+)dt (by Bayes’ rule)
  ⇒ p(S |S+;θfull) = ∫ ∞ 0 p(S | t;θ) p(t |S+;θ+)dt (dividing with p(S+;θ+)). (2)
  Comparing equations (1) and (2), we see that the information about S gained by observing S+ can be explained via a posterior observation time distribution p(t |S+;θ+) that refines the prior p(t) by taking into account S+.
  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 θ+, that is, Θ+ = θ+Im. In this case, |S+| | t follows a binomial distribution over m variables with success parameter 1−e−w+t, wherew+ := eθ+ . Intuitively, we expect these independent items to act as a “background clock”: counting the number of items |S+| 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+ ⊆ V be randomly drawn according to the CTMC with parameter matrix Θ+, and let t∗ be the true observation time of S+. Then, for any δ ∈ (0, 1), there exists an m0, such that for all m ≥ m0, the mean and variance of the posterior time distribution p(t |S+;θ) can be bounded as follows with probability at least 1− δ,
  |Mpost − t∗| ≤ C1(θ+, t∗) √ logm
  m +O
  ( logm
  m ) Vpost ≤ C2(θ+, t∗) 1
  m +O
  ( 1
  m2
  ) .
  The quantities C1 and C2 are constant w.r.t. m and encode how suitable a particular rate (quantified through θ+) of an item is for estimating the specific observation time t∗. Intuitively, larger rates are better suited for estimating smaller times and vice versa. We further illustrate this effect in Appendix B together with the detailed proof of the theorem. Note that for m → ∞ we get p(t |S+;θ+) → δ(t− t∗) in distribution, and therefore p(S |S+;θfull)→ p(S | t∗;θ). This means that observing an arbitrarily large amount of independent items is equivalent to knowing the true observation time. In Appendix C we show that a similar asymptotic convergence behavior can be observed even when the items in V+ are not identically distributed or exhibit some dependencies.
sections.heading: 3.3 Discussion
sections.text: To illustrate how this result applies to our two-item example discussed before, we define Θ+ = Im[θ (1) + , . . . , θ (m) + ] T , where each θ(i)+ is drawn uniformly in [−4,−2]. This corresponds to individual
  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 θfull. Figure 1 shows the posterior time density corresponding to each sample with S = {1} (red) and S = {1, 2} (blue), for different numbers m of independent items. To avoid clutter, we do not show the densities for samples with S = ∅. As expected, whenm = 0 all samples follow the prior time density, and we have no way to distinguish their time order. However, as 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 1 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.
  Due to computational considerations, it is common practice to constrain the ground set to a small subset of items that interact strongly and are deemed interesting to the application at hand; for example, genetic alterations that are known to have a role in cancer development (Raphael & Vandin, 2015; Schill et al., 2019). Our results in this section challenge this practice. We show that items believed 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, we will not assume any prior knowledge about the block-diagonal structure of Θfull or the values of the parameters θ+. Furthermore, we will make no distinction between “items of interest” and independent items. If there is indeed a block-diagonal parameter structure in the data, the addition of a regularization term that promotes sparsity can help recover this structure in 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 L1-regularized marginal log-likelihood,
  F (D;θ) = 1 N N∑ d=1 log p(S(d);θ)− λ ∑ i̸=j |θij |. (3)
sections.heading: 4 Efficient approximate likelihood maximization
sections.text: To put the previously discussed theoretical insight into practice, we need to be able to efficiently maximize 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 log p(S;θ) with respect to the parameters θ. Computing the exact gradient requires O(2n) 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.
sections.heading: 4.1 Full sequences
sections.text: First, we consider a setup in which we observe the arrival order of the items in the CTMC, which corresponds to building up an ordered sequence σ instead of a set S in Algorithm 1. Furthermore, we assume that tobs = +∞, which means that the returned sequence will be a permutation of the ground set V . By definition of the CTMC, the probability of observing a sequence σ = (σ1, . . . , σn) ∈ SV , where SV denotes the permutation group over V , can be written as
  p(σ;θ) = n∏
  i=1
  P ( Xi = σ[i] | Xi−1 = σ[i−1] ) = n∏ i=1 qσ[i−1] σ[i](θ) q̃σ[i−1](θ) . (4)
  We denote subsequences using the bracket notation, σ[i] := (σ1, . . . , σi), with σ[0] = (). Also, with a slight abuse of notation we directly use sequences in some places that require sets; therefore, the notation q̃σ is to be understood as q̃set(σ), where set(σ) := {σi | i ∈ {1, . . . , |σ|}}.
sections.heading: 4.2 Partial sequences given time
sections.text: As a next step, instead of tobs = +∞, we assume that we are given the observation time tobs, and would like to derive the probability of observing a partial sequence σ = (σ1, . . . , σk), with k ≤ n, at 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 k ordered elements agree with our partial sequence σ. Event B identifies the outcomes for which the observation time tobs falls between the addition of the k-th and k + 1-th elements; that is, Tk < tobs < Tk+1. We can then express the desired probability as p(σ | tobs;θ) = P (A)P (B | A). The probability of event A can be directly written using (4), except that, in this case the product is taken from i = 1 to k, instead of n. Regarding event B, first note that the sequence of events( {Ti < tobs | A} )n i=1 is decreasing. As a result, we can write
  P (B | A) = P (Tk < tobs < Tk+1 | A) = P (Tk < tobs | A)− P (Tk+1 < tobs | A) . (5)
  By definition of the CTMC, each holding time Hi is an exponentially distributed random variable whose parameter depends on the i − 1 items that have been added up to that point. We can also define the jump time Tk := ∑k i=1Hi, which represents the time at which the k-th change of state occurs. When we condition on knowing the first k items of the sequence, i.e., event A, the jump times Tk and Tk+1 are distributed as sums of k and k+1 independent exponential random variables respectively. The CDF of the sum of r independent exponential variables with rates λ1, . . . , λr, also known as the hypoexponential distribution (Bibinger, 2013), is given by
  FHEXP(y;λ1, . . . , λr) =
  ( r∏
  i=1
  λi
  ) r∑
  i=1
  1− e−λiy λi ∏ j ̸=i(λj − λi) .
  In our case, the rates λi of each exponential distribution are q̃[σi](θ), thus we can compute the terms of (5) as P (B | A) = FHEXP(tobs; q̃[σ1], . . . , q̃[σk])− FHEXP(tobs; q̃[σ1], . . . , q̃[σk+1]).
sections.heading: 4.3 Marginal partial sequences
sections.text: Next we consider the probability of observing a partial sequence σ = (σ1, . . . , σk), with k ≤ n, without knowledge of the observation time tobs. We define events A and B as before, except that in
  this case the observation time Tobs is a random variable, and event B is defined by Tk < Tobs < Tk+1. In particular, we assume that Tobs is an exponential random variable with rate λobs = 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 partial sequence σ = (σ1, . . . , σk) can be written as
  p(σ;θ) =
  ( k∏
  i=1
  qσ[i−1] σ[i](θ) 1 + q̃σ[i−1](θ)
  ) 1
  1 + q̃σ[k](θ) . (6)
  The proof can be found in Appendix D, but the form of this expression allows for an intuitive interpretation. We can modify the original continuous-time Markov chain by adding an extra dummy state Xobs, such that for all X ⊆ V , qX Xobs(θ) = 1, and qXobs X(θ) = 0. The first property implies that there is a fixed rate 1 to transition from any state to Xobs, while the second property implies that Xobs is a terminal state. It is easy to see then, that eq. (6) expresses exactly the probability of observing a “full” sequence σ in this modified chain. Note that a “full” sequence in this chain is no longer a permutation of V , but rather any sequence that ends with Xobs.
sections.heading: 4.4 Marginal sets
sections.text: Finally, we consider the probability of observing a set S ⊆ V , without knowledge of the order in which the individual items arrived, and still under the assumption that the observation time Tobs is an exponential random variable with rate λobs = 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;θ) = ∑ σ∈SS p(σ;θ). 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;θ). This gradient can be written as
  ∇θ log p(S;θ) = 1
  p(S;θ) ∇θp(S;θ) =
  1
  p(S;θ) ∑ σ∈SS ∇θp(σ;θ) = ∑ σ∈SS p(σ;θ) p(S;θ) ∇θ log p(σ;θ).
  In the last expression, we observe that the gradient corresponding to each permutation σ is weighed by the probability of observing that permutation of the given set S. This suggests a stochastic approximation of the desired gradient by first sampling m permutations σ(1), . . . , σ(M) according to p(· |S;θ) := p(· ;θ)/p(S;θ), and then computing the average of the resulting gradients,
  ∇θ log p(S;θ) ≈ 1
  M M∑ i=1 ∇θ log p(σ(i);θ). (7)
  Conceptually similar approaches have been used to approximate the gradient of the log-normalizer when maximizing the likelihood of energy-based probabilistic models (Song & Kingma, 2021).
  It remains to show how to obtain samples from p(· |S;θ), which we do by using a Markov chain Monte Carlo method on SS . More concretely, we employ a Metropolis-Hastings chain (Levin et al., 2009); at each time step, given the current permutation σ, the chain proposes a new permutation σnew according to proposal distribution Q(σnew |σ), and transitions to σnew with probability
  paccept = min ( 1, p(σnew |S;θ)Q(σ |σnew;θ) p(σ |S;θ)Q(σnew |σ;θ) ) .
  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 is based on the form of the marginal probability (6), and leads to faster and more stable learning in practice. While a detailed discussion and theoretical analysis of mixing is beyond the scope of this paper, in Appendix E we present our proposal in detail, and compare it to the uniform.
sections.heading: 5 Further related work
sections.text: There has been a long line of research focused on understanding and reconstructing the evolutionary history of tumors (Nik-Zainal et al., 2012; Welch et al., 2012; Turajlic et al., 2015; Mitchell
  et al., 2018; Jolly & Van Loo, 2018). However, these papers do not attempt to systematically model the timed interaction of genetic alterations across a data set. More recently, Gerstung et al. (2020) showed that modern data resulting from deep sequencing procedures may contain partial order information about specific genetic events. Incorporating this additional information into an analysis like ours, could potentially help with further reducing underspecification. Finally, there is evidence that specific classes of mutational processes generate mutations at a constant rate, and are largely unaffected by other factors (Alexandrov et al., 2015; Campbell et al., 2020; Alexandrov et al., 2020). The resulting so-called mutational signatures have been likened to a biological clock, and could be candidates for explicitly defining the set V+ discussed earlier (cf. Section 3.2).
sections.heading: 6 Experiments
sections.text: The following experimental setup is common to all our experiments. To optimize the objective (3), we use a proximal AdaGrad method (Duchi et al., 2011). We fix the initial step size to η = 1, and the regularization weight to λ = 0.01. To initialize the parameter matrix Θ, 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 in addition to 10 burn-in samples that are discarded. The code used to run our experiments can be found at https://github.com/3lectrologos/time/ tree/clean.
sections.heading: 6.1 Synthetic data
sections.text: We start by evaluating our approach on two synthetic data sets drawn from known true models. The first data set is identical to the one discussed in Figure 1, and contains observations of two items that tend to follow a specific time order (item 1 comes before item 2). The second data set contains observations of five items following the graphical structure shown on the upper right of Figure 2. Edges with arrows ← represent positive interactions (θij > 0), while edges with vertical lines ⊢ represent negative interactions (θij < 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 to examine the effect of adding a number m of independent items to the data. To quantify the quality of recovery, we compare the true distribution ptrue of marginal sequences, computed using eq. (6), to the distribution of marginal sequences p∗ of the learned model, approximated by drawing 106 marginal sequences from that model. In Figure 2 we depict the KL divergence dKL(p∗ ∥ ptrue) for increasing numbers m 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 a dashed line in the figure. We repeated each experiment 96 times; sources of randomness include the choice of independent items, the 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 a large enough number of additional items is practically equivalent to knowing the true observation times, and thus validates Theorem 1.
sections.heading: 6.2 Real cancer data
sections.text: We next evaluate our approach on a data set that contains N = 378 tumor samples of glioblastoma multiforme, an aggressive type of brain cancer. The data is part of the TCGA PanCancer Atlas project (https://www.cancer.gov/tcga), 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 a ground set of n = 410 items, each of which represents a point mutation, amplification, or deletion of a specific gene. Analyzing the interaction structure between such genetic alterations is of fundamental importance to cancer research, as it can help illuminate the processes that are involved in cancer initiation and progression.
  For the remainder of this section, we assume that the items in V are ordered by decreasing frequency, and when selecting a subset we keep the most frequent items. In Table 1, we show that our learning approach runs about 1000 times faster on a 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 to evaluate how some of the learned parameters behave as we increase the size of the ground set. Given a pair of items a, b ∈ V , and n−2 other (most frequent) items, we first learn a CTMC model on the n total items using our approach. Next we approximate the sequence distribution of just a and b (marginally over time and over all other items) by sampling 106 marginal sequences from the learned model. Finally, we compute the proportion of sequences in which a occurred before b, 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 3 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 in 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 cooccurrence 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 of the PDGFRA alterations can be robustly inferred only after including more than 70 additional items in the analysis.
  Second, we confirm some previously reported interactions, for example, the tendency of IDH1 mutations to occur before TP53 ones, as shown in the bottom left of the figure. This was observed before by Schill et al. (2019), and is also supported by some biological evidence (Watanabe et 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 to a fortuitous choice of the optimization setup, e.g., fixed initialization (see Appendix G). The final result in the bottom right of the figure seems to support this claim: the amplifications of MDM2 and CDK4 do not seem to have a clear time order, even though previous work has inferred both that CDK4(A) occurs before MDM2(A) (Cristea et al., 2017), as well as the opposite order (Schill et al., 2019).
  Finally, the parameter matrix Θ learned from this data set (see Appendix G) has indeed an approximately block-diagonal structure, which confirms the practical relevance of Theorem 1.
sections.heading: 7 Conclusion
sections.text: Our experiments demonstrated that, in the presence of underspecification, special care is required when interpreting results that may be significantly affected by particular choices in the optimization setup. In this paper, we proposed an 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 underspecification and learning robust time evolution models.
  Acknowledgments and Disclosure of Funding
  This work was partially supported by a Swiss National Science Foundation Early 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 (1R35CA220523).
title: Scaling up Continuous-Time Markov Chains Helps Resolve Underspecification