Dynamic Distribution Decomposition for Single-Cell Snapshot Time Series Identifies Subpopulations and Trajectories during iPSC Differentiation

Recent high-dimensional single-cell technologies such as mass cytometry are enabling time series experiments to monitor the temporal evolution of cell state distributions and to identify dynamically important cell states, such as fate decision states in differentiation. However, these technologies are destructive, and require analysis approaches that temporally map between cell state distributions across time points. Current approaches to approximate the single-cell time series as a dynamical system suffer from too restrictive assumptions about the type of kinetics, or link together pairs of sequential measurements in a discontinuous fashion. We propose Dynamic Distribution Decomposition (DDD), an operator approximation approach to infer a continuous distribution map between time points. On the basis of single-cell snapshot time series data, DDD approximates the continuous time Perron-Frobenius operator by means of a finite set of basis functions. This procedure can be interpreted as a continuous time Markov chain over a continuum of states. By only assuming a memoryless Markov (autonomous) process, the types of dynamics represented are more general than those represented by other common models, e.g., chemical reaction networks, stochastic differential equations. Additionally, the continuity assumption ensures that the same dynamical system maps between all time points, not arbitrarily changing at each time point. We demonstrate the ability of DDD to reconstruct dynamically important cell states and their transitions both on synthetic data, as well as on mass cytometry time series of iPSC reprogramming of a fibroblast system. We use DDD to find previously identified subpopulations of cells and to visualize differentiation trajectories. Dynamic Distribution Decomposition allows interpreting high-dimensional snapshot time series data as a low-dimensional Markov process, thereby enabling an interpretable dynamics analysis for a variety of biological processes by means of identifying their dynamically important cell states. Author summary High-dimensional single-cell snapshot measurements are now increasingly utilized to study dynamic processes. Such measurements enable us to evaluate cell population distributions and their evolution over time. However, it is not trivial to map these distribution across time and to identify dynamically important cell states, i.e. bottleneck regions of state space exhibiting a high degree of change. We present Dynamic Distribution Decomposition (DDD) achieving this task by encoding single-cell measurements as linear combination of basis function distributions and evolving these as a linear system. We demonstrate reconstruction of dynamically important states for synthetic data of a bifurcated diffusion process and mass cytometry data for iPSC reprogramming.

Data-driven reconstruction of dynamic processes constitutes a central aim of systems biology. High-dimensional 2 single-cell molecularly resolved time series data is becoming a key data source for this task [1,2]. However, 3 these technologies are destructive, and consequently result in snapshot time series data originating from 4 batches of cells collected at time points of interest. A longstanding and still challenging problem is to 5 reconstruct dynamic biological processes from this data, to the end of identifying dynamically important 6 states, i.e. regions of state space that cells preferentially pass through, these include transitionary states 7 (bottlenecks) when differentiation decisions are made, and terminal states. With snapshot time series, it is 8 challenging to identify these states as we cannot track the state of an individual cell at one time point to a 9 new state at a later time point, one has to temporally map between state distributions. 10 Chemical reaction networks are a popular class of parametric models assuming that the temporal state 11 evolution is well described by chemical kinetics. Ordinary differential equations (ODEs) are used to describe 12 smooth deterministic dynamics, and stochastic differential equations (SDEs) for dynamics in the low copy 13 number/concentration regimes affected by stochastic fluctuations. Chemical reaction network models require 14 explicit definition of the model structure, i.e. set of reactions or interactions among the system components. 15 This task is manageable for small, well defined systems, such as small signaling systems [3]. However, by 16 means of high-dimensional measurements, we typically observe larger systems comprising at least dozens of 17 components with largely a priori undefined interactions. This situation results in a combinatorial explosion 18 of model variants that cannot be exhaustively evaluated [4,5]. Alternative approaches are agnostic with 19 regards to parametric form and model structure and use a probabilistically-motivated rule to map between 20 distributions, e.g., one optimal transport method maps neighbours at one time point to the nearest neighbour 21 at the next time point [6]. However, such generic approach are rather extreme in their agnosticism and 22 abandon reasonable assumptions on the dynamics of cellular systems, e.g., that cells can be modelled as an 23 autonomous dynamical systems in continuous time, such as a Markov chain, where the cell's current state 24 infers its likely future state independent of the current time within the experiment. 25 Operator approximation methods constitute an alternative class of models that are agnostic to model 26 structure and yet allow for encoding of general system properties such as autonomy, conservation of mass, 27 and boundary conditions. These methods approximate both the Perron-Frobenius operator [7] and the 28 Koopman operator [8]. These operators describe the evolution of distributions and other functions of a 29 dynamical system's state. The early theory on these operators was developed to describe systems in classical, 30 statistical, and quantum mechanics [9][10][11][12][13], and in probability theory [14,15]. The operators fully describe a 31 nonlinear dynamical system as a linear system in higher, possibly infinite, dimensions. Hence, techniques 32 from linear analysis can be utilised to gain insight from the systems; in particular, calculation of eigenvalues 33 and eigenfunctions allow for timescale separation. Eigenfunctions of linear operators show the fundamental 34 building blocks of possible behaviours available to a dynamical system, e.g., exponential growth/decay, 35 oscillations, a steady state. Data-driven approximations to the operators have been investigated in the past 36 years, originating in the computational fluid dynamics community [16][17][18][19][20]. Their focus is to approximate 37 finite-dimensional projections of the Koopman operator with a family of algorithms known as Dynamic 38 Mode Decomposition (DMD). The algorithm has further been applied to other areas such as neuroscience, 39 infectious disease epidemiology, and control theory [21][22][23] and parameter estimation [24,25]. When carrying 40 out approximations of these operators, eigenvalues can be ordered in terms of magnitude to extract slow 41 behaviour (approximated well) to fast behaviour (representative of noise) [25]. Dynamic Mode Decomposition 42 assumes that the data is recorded at equally spaced time points, whilst our work extends the technique to 43 support data recorded at arbitrary time points. 44 We adapt Dynamic Mode Decomposition to identify dynamically important states from single-cell snapshot 45 time series. Our method is based on representing the distribution at each time point via basis functions 46 and calculating an approximation to the Perron-Frobenius operator by minimising an error term -akin 47 to least squares when fitting ODEs. The error terms then give an indication of how well the data fits 48 into the model assumptions: primarily that the data is generated by an autonomous dynamical system. 49

2/22
Because the Perron-Frobenius operator describes the evolution of distributions, we name our approach 50 Dynamic Distribution Decomposition (DDD) in keeping with the DMD naming convention. DDD leads to 51 the calculation of a Markov rate matrix but over a continuum of states -as opposed to discrete states. 52 As previously mentioned, one can then use standard methods of analysis for linear operators based on 53 eigen-decompositions. As Markov processes can be represented as directed weighted graphs, our graph can 54 be evaluated in two dimensions and then the high-dimensional operator and its corresponding eigenfunctions 55 have a natural low dimensional representation. Our approach also allows for visualisation of inferred state 56 trajectories as a branching structure when cell fates are stochastic, and approximation of fitting error when 57 matching model prediction to sample data. By using a Markov rate approach over distributions, we overcome 58 the difficulties listed above. We demonstrate DDD on a synthetic stochastic dynamical system representing 59 cells making a cell differentiation decision as well as for a mass cytometry time series taken from an iPSC 60 reprogramming of a fibroblast cell line taken from Zunder et. al. [26] to re-identify subpopulations of cells as 61 first elucidated in the original manuscript.  We developed a method herein referred to as Dynamic Distribution Decomposition to analyse snapshot 66 time series, consisting of the following stages (illustrated in Figure 1 Markov rate matrix; the structure of this matrix is often dense but its dominating structure can be elucidated 73 via Lasso regularisation, see Figure 1(f). We applied our method to two systems: first, simulated particles in 74 a potential well; and second, experimental data of iPSC reprogramming of a fibroblast system.

75
Particles in Potential Well with Fluctuations 76 The first numerical example is for illustrative purposes whereby we know the stochastic process generating the sample points. We consider simulated particles in a bistable potential well undergoing fluctuations. After initialisation around point (1, 1) /2, particles stochastically switch between one of two paths: y = 2x or y = x/2 to finally settle in one of the two final state (2,4) or (4,2) . We model this process by the two-dimensional SDE {X t = (X 1t , X 2t ) ∈ R 2 where W t is a two-dimensional Wiener process. The potential well is of the form As an initial condition at t = t 1 , the sample is placed with a multivariate normal distribution with mean 77 µ = (1, 1) /2 and covariance matrix Σ = I 2 /2 where I 2 is the identity matrix. The diffusion constant is 78 chosen to be D = 1/4. Along the lines x = 0 and y = 0, the system has reflecting boundaries imposed. 79 For simulations, the Euler-Maruyama numerical scheme is used 1 with time step δt = 2 −9 . Three sample 80 trajectories are visualised, see Figure 2(a); and the potential well is plotted, see Figure 2(b).

4/22 a.)
Eigenfunction plots Graph representation d.) x 1 x 2 e.) x 1 Dynamic distribution decomposition applied to data generated by stochastic dynamical system described by equation

86
The eigenfunctions of the approximated Perron-Frobenius operator allow us to identify the steady state and 87 bistable paths in the above system. In low dimensions we visualise eigenfunctions of P as a continuous function, 88 see Figure 3(a,b,c); or a graph, see Figure 3(d,e,f) and Methods section. Eigenfunctions corresponding to 89 eigenvalues with large absolute value are approximated with larger error than in cases with small eigenvalues, 90 a point also noted in Ref. [25]; notice in Figure 3(c,f) there are a few small fluctuations around (1, 1) /2. 91 Also, eigenvalues and eigenfunctions are basis function dependent, so changes in basis functions change the 92 eigen-decomposition. However, regardless of changes to the basis functions, the key dynamic states (as visible 93 in the eigenfunctions) remain the same provided the changes to the basis functions are not drastic. Since 94 we use 30 basis functions, hypothetically we can find 30 eigenfunctions. However, we just plot the first 95 three eigenfunctions; these are real with no imaginary component. From these figures, it is clear that three 96 basins around (2,4) and (4,2) and at the initial condition (1, 1) /2 are dynamically important. Therefore, 97 examination of the first few eigenfunctions allows for detection of dynamically important states.

99
We evaluated the robustness of our inference procedure to measurement noise. Specifically, three further 100 modified data sets are also considered: (i.) considering a single time point perturbed by additive random 101 noise; (ii.) removing this perturbed time point and fitting the model; and (iii.) randomly perturbing all time 102 points by additive random noise. For all perturbations, sample points are modified by an additive error term 103 drawn from a zero-mean multivariate Gaussian with covariance matrix Σ = I 2 /4.

5/22
a.) decomposition applied to data generated by stochastic dynamical system described by equation (1)- (2). The representation of P as a Markov rate matrix, colours of edges denote magnitude of rate change from one node to another, and node colours denote rate at which node decays. Colours scaled to unit interval. for r = 1, . . . , R, see Figure 2(c). We find that all data sets have consistently low error, but with an 106 increase in error at the beginning of the realisation; this is due to the boundary conditions which were 107 not incorporated into the choice in basis functions, which one would typically do when solving PDEs via a 108 Galerkin approximation. The data set perturbed at time t = 8 (dashed red line) leads to increased error 109 immediately before and after this time point (circled in black); after removing this erroneous data (fitting 110 without t = 8) one obtains reduced errors comparable to the original data set (dotted red line). When adding 111 systematic error to all time points (dashed blue line), one observes similar errors to the original data set; the 112 reason for this is that the Gaussian basis functions now have a covariance matrix with larger entries (whilst 113 the means remain similar). The eigenfunction plots are also similar to those generated by the original data set 114 but more spread out (not plotted). Therefore we see that DDD is robust with regards to noisy observations. 115 Lasso Regularisation Reveals Sparse Topology 116 We utilize Lasso regularization to identify key transition states, see Methods section. Specifically, P as a 117 Markov rate matrix with the nodes located at the mean of the components of the Gaussian mixture model. 118 The resulting network is cluttered and is hard to identify meaningful states or transition, see Figure 4(a). 119 Lasso regularisation encourages sparsity and reveals the simple underlying structure, see Figure 4(b). The 120 skeletal structure shows that around the initial condition the particle becomes strongly committed to one 121 branch over the other -an accurate reflection of the dynamical system.

122
Mass Cytometry Data: iPSC Fibroblast Reprogramming 123 We studied the process of iPSC reprogramming using Dynamic Distribution Decomposition. We considered 124 data from a study established by Zunder et. al. [26]. Specifically, the reprogramming of a fibroblast cell 125 line differentiating into an induced pluripotent stem cell state was studied using mass cytometry. Cells 126 are labelled using mass-tag cell barcoding, stained with antibodies before being measured via CyTOF. We 127 focus our study to the cell line with the largest amount of cell events, i.e. on a Nanog-Neo secondary mouse 128 embroyic fibroblasts (MEF) that expresses neomycin resistance gene from the endogenous Nanog locus. 129 Reprogramming was monitored by Dox induction for 16 days followed by subsequent addition of LIF; the 130 experiment was carried out over 30 days. Experiments were initialised together and cells harvested every 2 131 days until 24 days with a final measurement taken at the final time point; 18 protein markers were used as 132 proxies to measure pluripotency, differentiation, cell cycle status, and cellular signalling. 133 We now briefly state DDD for choosing basis functions. In the synthetic data example, we used prior 134 information that there were 3 clusters, so used a 3 component Gaussian mixture model for each time (therefore 135 N = 30 for 10 time points). For non-synthetic data we do not necessarily have this information, therefore we 136 developed an approach to choose an expressive set of basis functions without letting their number grow too 137 large and thereby ensure efficient solving of the minimisation problem later presented in equation (17). We 138 fit multiple Gaussian mixture models to each time point, varying the number of components until the AIC 139 curve flattened out [27]; in our case this happens at approximately 8 basis functions per time point. To avoid 140 overfitting, we use regularisation to specify minimum diagonal entries of the covariance matrix. As our data 141 has been scaled via the commonly used transformation function f (x) = arcsin(x/5) and then standardised by 142 z-scoring, we choose a regularisation value of 1/2; smaller values can be used should one wish to capture 143 sharp peaks, but at the cost of additional basis functions. Finally, for each time point we cluster the data 144 into these Gaussian mixture models and remove poorly populated components. Here, after clustering, we 145 remove any basis functions which represent less than α × 100% of the data. Therefore, it should be noted 146 that obtaining a good fit to the matrix P is a payoff between: (i.  We now decrease α and evaluate whether we have sufficient basis functions. We plot the percentage 151 fitting error at each time point and the mean percentage error as a function of α, see Figure 5(a,b). These 152 figures show that as α decreases, the error only minimally decreases for large increases in the total number of 153 basis functions N . We can also view the eigenvalues plotted in the complex plane for various values of α, see 154 Figure 5(c). We rescaled time to the unit interval, therefore one will not be able observe eigenfunctions with a 155 corresponding non-zero eigenvalue (λ) > −1, i.e., we cannot observe timescales slower than the observation 156 window. We notice for α = 0.005 that (λ 1 ) = −1.06 so we are confident decreasing α will not offer much 157 benefit. Additionally in the cases where α = 0.01 and α = 0.005, the extrema of the first few eigenfunctions 158 correspond to the same basis functions, not plotted. 159

7/22
Loss of Dynamic Autonomy after Stimulus Removal 160 When a stochastic dynamical system is autonomous, the current state of the system determines the likely 161 future states; here we show that after and including t = 16 days the system becomes less autonomous, once 162 the Dox induction had ended. We see the error plotted again at each time point for α = 0.005, see Figure 163 5(d). We notice that time points after and including t = 16 days contain the vast majority of fitting error. 164 To rule out the possibility that the dynamical system instantaneously changed at t = 16 days, we fit two 165 Perron-Frobenius matrices, one using the first 8 time points and a second using the last 6 time points with 166 all fits using the same basis functions. We find that there is still much more error contained in the final 6 167 time points compared to the first 8.

168
The autonomous dynamical system assumption means that using the data presented, the future states 169 of the system depend on the current state. While this is likely true within a cell culture system, we only 170 observe a tiny fraction of the state space of the dynamical system as we do not measure the transcriptome 171 and the vast majority of the proteome. Therefore, it seems reasonable to assume that from t = 16 days, we 172 are not observing enough of the dynamical system to obtain a linear map between distributions. This insight 173 suggests further single-cell experiments at these later time points using technologies allowing greater 'omic' 174 profiling, e.g., single-cell RNA-Seq.

Inferred Dynamically Important States Agree with Previously Described Cell Subpopulations 176
We evaluated the extreme of the eigenfunctions of the approximated Perron-Frobenius operator to re-identify 177 cell subpopulations found in Zunder et. al. [26]. We first plot the first 6 eigenfunctions, see Figure 6. Nodes 178 that are close together to each other in 18 dimensions (using Euclidean distance) as plotted as close to each 179 other in 2 dimensions. Protein expression of the basis functionsare also plotted using the same coordinates as 180 the graph, see extra figure in Appendix.

181
When examining the extrema of the eigenfunctions, basis functions seem to cluster in 3 groups: group A 182 centred around basis function 32; group B with members 56, 61, and 65; and group C with one member, 183 basis function 66. Our algorithm recovers the same populations as stated in Zunder et. al. [26]: cells with 184 low Ki-67 expression do not successfully reprogram and remain MEF-like (group A); cells with high Ki-67 185 expression then subdivide into two populations, an embryonic stem cell-like (ESC-like) population with 186 Nanog + , Sox2 + , and CD54 + (group C) and a mesendoderm like population with Nanog − , Sox2 − , Lin-28 + , 187 CD24 + expression (group B). As our basis functions were added sequentially per time point, the MEF-like 188 population appeared first.

189
DDD suggests a few new insights previously not elucidated in Zunder et. al. [26]. We find according to 190 the fitted Perron-Frobenius operator, MEF-like cells form the steady state (when λ = 0). Therefore, the 191 model predicts all cells would revert to fibroblasts if enough time passes -although one has to be careful 192 over interpreting predictions due to the higher error after t = 16 days.

193
Lasso Regularisation Reveals Sparse Topology of iPSC Dynamics 194 Reminiscent of the SDE example, the graph as induced by the transition matrix P is cluttered due to an 195 abundance of low weighted edges, see Figure 7(a); we apply the Lasso modification to reveal a two branching 196 points, see Figure 7(b) and Methods section. Finally, to focus on the 3 groups previously identified, we prune 197 edges leading to unannotated nodes to obtain an easily interpretable branching structure, see Figure Figure 198    described by basis function 29 was previously described in Zunder et. al. [26], but we are able to include an 205 additional transitionary state by means of basis function 16. We can conclude that the cell decision towards 206 becoming remaining MEF-like is made early during the course of the experiment: basis function 16 was 207 placed with the data recorded at 6 days; basis function 29 was placed with the data recorded at 12 days.

209
We now give the mathematical set-up to our problem, additional technical details are given in the appendix. 210 The method follows the following steps: (i.) the statement is posed that the temporal evolution of cell states 211 follows a linear partial differential equation; (ii.) the distribution of the sample points at each time point 212 can be encoded into a sequence of basis functions; (iii.) the weights of these basis functions can change 213 dynamically interpolating between sample points; and (iv.) we fit the form of the matrix approximation of 214 differential operator around these changing basis functions; and (v.) study the eigenfunctions. The workflow 215 of the method is also as an illustration in Figure 1.

217
Assume we have a sequence of R experimental readings at time points t 1 < · · · < t R ; without loss of 218 generality we choose t 1 = 0. At each time point, n r cells are harvested with states X r = {x r,1 , . . . , x r,nr } for 219 r = 1, . . . , R. The state of each cell is located in a (measurable) space, x ∈ M. We note that this space may 220 not be the full dimension of the data set, but after a dimensionality reduction technique has been applied, 221 e.g., PCA, diffusion maps etc. For example, in the case of RNA-Seq data, the state of a cell would consist of 222 thousands of genes which would be too high to apply kernels to. We wish to find a probability distribution 223 = (t, x) such that when t = t r , the probability of observing cells in states X r would be highly probable 224 for r = 1, . . . , R. 225

10/22
Immediately necessary to ensure conservation of mass, we require for all t ∈ (t 1 , t R ). We now make the crucial assumption that our method relies on: each cell follows a (well behaved) autonomous dynamical system -implicit in this assumption is that cells do not interact, alternatively cell interactions can be accounted for via stochastic noise terms. Under these assumptions, we can interpret (t, x)∆x as the probability a randomly selected cell has state in the interval [x, x + ∆x) at time t. We write down the (continuous-time) Perron-Frobenius equation for the dynamics of the density profile as for initial condition (t = 0, x) = 0 (x). The P term is the continuous-time Perron-Frobenius operator 226 [19,20]. This operator is known by many names depending on the underlying dynamical system for the state 227 evolution x(t) ∈ M for t ∈ (t 1 , t R ). For example, within SDEs equation (4) is a second order parabolic PDE 228 known as the Fokker-Planck equation [28]; and for chemical reaction networks equation (4) is a system of 229 coupled ODEs known as the chemical master equation [29].

231
We would like to find a finite dimensional approximation of of P; we can do this with non-negative basis functions ψ(x) = [ψ 1 (x), . . . , ψ N (x)] . We take the ansatz that for all t ∈ (t 1 , t R ), we can expand as the linear combination If the basis functions are themselves probability density functions, then Λ is the probability simplex. In 232 Figure 1(a), we show how a distribution can be represented as a sum of normal distributions, that is, the 233 density is given by a Gaussian mixture model. 234 We can derive a linear system of ODEs for the coefficients c(t). We do this by noting in weak form [30,31] equation (4) is Choosing g = ψ i and expanding as in equation (5), then where M ij = ψ i , ψ j and Q ij = ψ i , Pψ j . Assuming M is invertible, define which is the projection of P onto the basis functions {ψ 1 , . . . , ψ N }. That is, for g(x) = c ψ(x), we have the equality Pg(x) = (P c) ψ(x). In order to preserve probability density and positivity, we require P ∈ P for P = P ∈ R N ×N : N k=1 ω k P kj = 0 and P ij ≥ 0 for i, j = 1, . . . , N and i = j .
The explanation behind equation (10) is contained in the Appendix. 235 We can solve the dynamics for equation (4) using the approximation in equation (5) by using the matrix 236 exponential operation, specifically 237 c(t) = e tP c * , where c * are the coefficients at corresponding to the chosen initial condition 0 = c * ψ(x).

238
Selection of P matrix 239 We now need to address the issue of how to determine P from data. Consider a linear operator P with matrix 240 representation P on the space spanned by ψ. The L 2 norm gives a measure of how well P represents the 241 evolution of the densities.

242
For an initial condition of the squared relative prediction error at time t = t r for r = 1, . . . , R is Notice that we specify that the error is a function of both the Perron-Frobenius operator and the initial condition, which is then treated as a parameter of the model. We then define the mean squared relative prediction error as It would be ideal now to find {P, 0 } := arg min P, 0≥0, 0∈L 1 ε 2 (P, 0 ) = arg min Of course, we do not know what this error is without using our finite dimensional approximation; therefore, by using equation (5) we calculate the time t = t r error as and therefore, our objective function is modified by using this finite dimensional approximation to {P, c * } := arg min P ∈P, c * ∈Λ ε 2 (P, c * ) = arg min Unmentioned at this point is that: for a large quantity of basis functions the size over which this optimisation problem occurs is challenging. That is: one has N 2 free parameters in the matrix P , with zero
2: return P column sums one has N (N − 1) degrees of freedom; but one also has the initial condition to choose adding N parameters, so N − 1 degrees of freedom with unit column sum -in total (N − 1)(N + 1) degrees of freedom. Therefore, the problem is unapproachable without gradient calculations to speed up the optimization algorithm. Using the exponential matrix derivative (see Appendix), one can calculate the t = t r relative error with respect to P as and the derivative with respect to c * as The terms can be calculated using the recursion relation Lasso Regularisation

243
For display purposes, we can promote sparsity in P by using a Lasso regularisation. We modify the error term in equation (17) to where • denotes the Hadamard product (or entrywise product) and vec(·) denotes the vectorisation of a matrix. In the case where the basis functions are probability density functions, we can calculate the derivative of this expression as

13/22
By using the mass matrix M as a weighting in front of the PF matrix P , we are promoting edges between 244 basis functions located apart from each other. One can interpret the matrix P as a Markov rate matrix, in which case the entry P i,j shows the rate at which 248 state j transitions into state i. Therefore, a cell in cluster i will switch to cluster j in time interval [t, t + ∆t) 249 for ∆t > 0 with probability P i,j ∆t + O(∆t 2 ). To reiterate, cells exist in states between the basis functions so 250 instead of being at a single state, they are in a state which is a weighted combination of the basis functions. 251 However, this interpretation allows us to plot a directed network with weighted adjacency matrix P i,j . Nodes 252 can then be placed using one of a multitude of algorithms, in our case we use force-directed node placement 253 with weights inversely proportional to the mass matrix M . We also plot the size of node i proportional to 254 P i,i (as this is the rate at which state i remains in state i).

Eigenfunction visualisation 256
To investigate key dynamical behaviours of a linear operator, a common theme is the study of the corresponding 257 eigenproblem. By solving the eigenproblem, one can decompose the solution of the operator into components 258 (known as eigenfunctions or eigenvectors) that will dynamically change with respect to the eigenvalue. 259 By studying the eigenproblem, one can break down the solution into key behaviours and find important 260 transitionary states.

261
For an eigenfunction satisfying using the finite dimensional Galerkin approximation, there is the corresponding eigenvector In low dimensions, one can simply plot this function as a linear combination To ensure consistent scales when plotting, we demand or in our finite dimensional representation In high dimensions, our ability to visualise functions is limited. However, we have represented the function as a linear combination of basis functions and so we only need to present the coefficients in the eigenvector.
To visualise if the eigenvector values have similar or dissimilar values, we can consider representing the eigenfunction as a graph. We specify the adjacency matrix for an undirected weighted graph as the outer product

14/22
For Perron-Frobenius eigenfunctions with eigenvalue λ = 0, the function will have both a positive part and a 262 negative part (indicating where probability mass is flowing from and to). Therefore, by examining G λ , a 263 positive value in entry (i, j) in G λ indicates basis functions i and j both have the same sign (positive or 264 negative) and a negative value indicates they have opposite signs. By using the weighting of M 1 2 in front of 265 the eigenvector v λ we ensure that the sizes of the eigenvectors are bounded. Occasionally it is the case that 266 we get complex eigenvalues, in which case they appear as complex conjugates and one can plot the real and 267 imaginary parts separately.

269
We presented Dynamic Distribution Decomposition for identification of dynamically important states of 270 biological processes. This method operates on snapshot time series data and infers dynamically important 271 states by mapping between distributions.. We applied our approach to synthetic data generated for a simple 272 test system, and then further to a mass cytometry time series data set for iPSC reprogramming. Our 273 approach performed well for both systems showing key dynamical states. For the experimental system of 274 iPSC reprogramming of fibroblasts, we could also identify key time points where the current experimental 275 set-up is insufficient elucidate the reprogramming process and where further investigation is warranted.

276
DDD can be computed efficiently, e.g., via the recursion relation given by equation (55), we have not 277 optimised the implementation. In the case where there is more than 100 basis functions, our minimisation 278 procedure using the inbuilt MATLAB multivariate minimisation algorithm can be unreasonably slow; therefore, 279 we will investigate more efficient implementations.

280
DDD depends on a few design decisions, such as the choice of basis functions. In this manuscript, we 281 used basis functions as components of a Gaussian mixture model and gave parameters that needed tuning to 282 alter the fit. Our method for choosing basis functions does not have an optimal configuration with regards 283 to minimising error. This is because by sending the regularisation value to infinity one obtains perfect fits, 284 and sending the regularisation value to zero can lead to ill fitting solution (basis functions with same mean 285 etc). However, one can use the α parameter as a rule of thumb to ensure enough basis functions are included. 286 For other applications other choices in basis functions are conceivable, for example, radial basis functions 287 [32]; piecewise linear basis functions [33]; and global basis functions [19,25] to name a few. When the basis 288 functions have finite support, the mass matrix M will be sparse, in which case the Lasso step will not be 289 necessary. It would also make sense to use basis functions built around the data type, for example negative 290 binomial distributions are often used to model UMI counts from single-cell RNA-Seq data. When one uses a 291 single basis function centred around each data point, one refers to this as kernel density estimation, of which 292 there are optimal methods to choose the basis function [34]; when using a small number of basis functions for 293 a large number of sample points, there are likely optimal ways of choosing them which we will investigate in 294 future work.

295
DDD could be applied to investigate pseudo-time ordered single-cell data of single-time-point experiments. 296 Here, one uses single-cell data measured at only a single time point to carry out trajectory inference and 297 subpopulation identification to infer biological processes, e.g., the cell cycle; a review of such methods can 298 be found in Ref. [35]. It may be possible to improve our fits by combining approaches: while cells are 299 monitored according experimental time coordinate, individual cell time coordinates might might deviate due 300 to asychronity of process initiation; this could be incorporated to get smoother Perron-Frobenius operators 301 between time points, see Ref. [36]. This would then be a biologically motivated method to account for delays 302 in the system.

303
The work presents Dynamic Distribution Decomposition, linking operator theory to the practical world of 304 high-dimensional data analysis. While we focus on application of DDD to mass cytometry measurements, 305 it is conceivable to expand to applications to single-cell RNA sequencing time series as well as biological 306 processes other than an iPSC reprogramming. We expect DDD and method variations will be instrumental 307 in providing intuitive understanding of such biological processes.

Additional Mathematical Details Choice in basis functions
One option for choosing the basis functions is to use probability density functions, so M ψ j (x)dx = 1 , for j = 1, . . . , N . We can find the values of c at the observed time points by noting that the value of the coefficient at time t = t r must be proportional to the probability that basis function j created the data at that time point, so One can then normalise N j=1 c j = 1 for j = 1, . . . , N to find the coeffient vector c r . Assuming multivariate Normal distributions as basis functions, we have where N is the probability density function for a multivariate Normal distribution with mean µ and covariance matrix Σ. To determine the mass matrix, adapting results from Ref. [37] we can analytically calculate the inner product between two multivariate normal distributions as for η 1 = Σ −1 1 µ 1 , η 2 = Σ −1 2 µ 2 , and η * = η 1 + η 2 .

Requirements on Perron-Frobenius matrix approximation
To preserve probability, it must be the case that N i=1 ω i P ij = 0 ∀j = 1, . . . , N , = ω P c(t) Therefore, for mass to be conserved, equation (10) holds.
We also require the operator e tP preserves positivity, that is, e tP u 0 ≥ 0 for functions u 0 ≥ 0. This implies that the off-diagonal entries of P must be non-negative. The reason is as follows. As ψ ≥ 0, the function u 0 = c u ψ must have c u ≥ 0 due to linear independence of the normalised basis functions. It therefore follows that e tP ≥ 0 element wise, because the positivity condition and linear independence of the normalised basis functions require that c u e tP ≥ 0 element wise. Thus, 0 ≤ e tP = I + tP + O(t 2 ) for all t ≥ 0, which in particular means that tP i,j + O(t 2 ) ≥ 0 for i = j. The condition P i,j ≥ 0 follows by considering arbitrarily small t > 0. This is sufficient for e tP to be element wise non-negative, as the matrix exponential of a matrix with non-negative off-diagonal elements is always non-negative element wise. Note that this, together with the mass conservation condition, implies that the diagonal entries of P are all negative.

Gradient Calculation
For relevant background reading on functions of matrices, see Ref. [38]. Let H denote the Hilbert space R N with inner product p, q H = p M q, P ∈ R N ×N . For t > 0, and p, q ∈ H, we derive the gradient of f (P ) = e tP q − p 2 H that is used for the gradient calculation in Equation (18). Let D f (P ; dP ) be the directional derivative of f at P , in the direction dP . Let ·, · F denote the Frobenius inner product on R N ×N .
The gradient of f is the matrix G ∈ R N ×N such that D f (P ; dP ) = G, dP F for any dP ∈ R N ×N .
We therefore need an expression for the directional derivative D e t· (P ; dP ) of P → e tP . By definition of the matrix exponential, we have that e t(P +dP ) = e tP + ∞ k=0 t k+1 (k + 1)! k j=0 P k−j (dP )P j + o( dP H→H ).