Distance-Driven Compositional Markov Models for Analyzing Time-Evolving Populations with Application to Single-Cell Dynamics

Distance-Driven Compositional Markov Models for Analyzing Time-Evolving Populations with Application to Single-Cell Dynamics Komlan Atitey1*, Benedict Anchang1 1Biostatistics and Computational Biology Branch, National Institute of Environmental Health Sciences, 111 T W Alexander Dr Rall Building, Research Triangle Park, NC 27709 E-mails: [email protected]; [email protected] Abstract Modeling time-evolving, heterogeneous populations as stochastic transitions among states is natural, yet the distance metric that governs neighborhood structure and transition probabilities, remains underexamined. We present a distance-driven compositional Markov framework that treats observations at each time point as probability vectors on the simplex and advances them via distance-parameterized transition operators. Distances (e.g., Euclidean, Manhattan, diffusion, or learned metrics) are elevated to first-class, testable components that shape state-state interactions. Our reference implementation, MarkovCellNet, integrates metric-aware kernels with time-informed embeddings and evaluates model fit using a multinomial likelihood, with a Bayesian extension for uncertainty quantification. This provides a unified objective for comparing metrics, embeddings, and regularization choices. Using single-cell transcriptomics as a motivating example, we benchmark on synthetic datasets that encode divergent, oscillatory, and convergent dynamics. Across scenarios, diffusion distances paired with geometry-preserving embeddings (e.g., PHATE/UMAP) yield higher likelihoods and more coherent longitudinal structures than Euclidean or Manhattan alternatives, demonstrating that metric-embedding synergy materially affects inference. Although instantiated here for single-cell data, the framework is modality-agnostic and applies to diverse high-dimensional, time-resolved systems, including ecological communities, longitudinal patient cohorts, online user populations, financial portfolios, transportation flows, and imaging series. Taken together, our results provide a general foundation for evaluating and selecting distance metrics in probabilistic transition modeling, with MarkovCellNet serving as a practical, extensible reference implementation. Key Words: Markov modeling, High-dimensional data, Dimensionality reduction, Temporal dynamics, Cell interactions, Trajectory inference 1. Introduction Modeling the temporal evolution of complex populations, whether cells, ecological species, patient cohorts, online communities, financial portfolios, or transportation flows requires probabilistic frameworks that capture stochastic transitions, structural heterogeneity, and gradual state changes(Van Kampen 1992; Levin and Peres 2017). To meet this need, we introduce MarkovCellNet, a distance-driven compositional Markov framework that elevates the distance (similarity) metric from a preprocessing convenience to a fundamental or first-class, testable parameter that directly governs transition probabilities. Unlike existing pipelines that fix a kernel, treat distance implicitly, or optimize heuristics without a proper likelihood (Coifman et al. 2005; Maaten and Hinton 2008), MarkovCellNet makes the role of distance explicit and testable. It provides a unified, modality-agnostic recipe that (i) converts metric-defined neighborhoods into rowstochastic transition operators (Coifman et al. 2005), (ii) evolves population states as probability compositions on the simplex (Aitchison 1982), and (iii) scores competing metrics and embeddings with a principled likelihood (Anderson and Goodman 1957; Aggarwal, Hinneburg, and Keim 2001). Together, these steps allow formal model selection and ablation, elevating distance choice from a hidden preprocessing step to a testable modeling hypothesis. This resolves a long-standing gap: to our knowledge, no prior method systematically quantifies how the choice of distance metric shapes probabilistic dynamics and interpretability over time. As illustrated in Figure 1, alternative metrics (Euclidean, Manhattan, diffusion) induce different neighborhood geometries that, when used to construct distance-driven transition kernels, generate different trajectory structures even on the same 2D embedding (Aggarwal, Hinneburg, and Keim 2001; Coifman et al. 2005). By making metric choice explicit and testable, MarkovCellNet transforms it from a hidden preprocessing step into a scientifically actionable parameter, thereby improving both explanatory power and predictive fidelity in modeling dynamic populations. Figure 1. Distance metrics govern neighborhood geometry and, through MarkovCellNet, the inferred population dynamics. Left to right: alternative metrics including Euclidean, Manhattan, and diffusion, define distance-sensitive transition kernels (green) on the same dataset. The data are embedded in 2D for neighborhood construction (middle) and passed to MarkovCellNet, which converts metric-defined neighborhoods into a row-stochastic transition operator and predicted compositions (right). Because each metric emphasizes a different geometry, Euclidean: isotropic proximity; Manhattan: axis-aligned separation; diffusion: manifold connectivity, the resulting trajectories and their interpretability differ. The workflow is modality-agnostic and applies to any high-dimensional, time-resolved population. Despite the growing availability of time-course data sets across different domains, many analytical pipelines sidestep how metric choice governs transitions. Clustering each time point independently discards temporal dependence and cannot quantify uncertainty about movement between states. Pseudotime and trajectory methods including Monocle (Trapnell and Cacchiarelli 2014), Slingshot (Street et al. 2018), and scVelo (Bergen et al. 2020) recover progressions but stop short of formalizing transitions as distance-scaled probabilities that adapt to scale, density, or anisotropy. Markov-based tools such as Palantir(Setty et al. 2019) and CellRank(Lange et al. 2022) advance a probabilistic perspective but typically adopt a fixed kernel structure, leaving the role of distance metrics implicit and largely unevaluated across time. In fact, traditional Markov models are agnostic to geometry: they operate directly on state-state transition probabilities without regard to how distances between states are defined. This gap matters, because alternative metrics (e.g., Euclidean, Manhattan, diffusion) induce different neighborhood geometries that can substantially alter inferred dynamics. Visualization frameworks such as PHATE (Moon et al. 2017) and PAGA (Zhang et al. 2023), preserve geometry but remain primarily descriptive, offering no generative, likelihood-based account of temporal evolution. Related lines of work in optimal transport (Schiebinger et al. 2019) and diffusion maps (Coifman et al. 2005) emphasize geometry, but computational intensity, sensitivity to cost definitions/bandwidths, and the lack of a compositional (simplex) interpretation limit their applicability for comparing distance choices or propagating uncertainty across discrete time points. To overcome these limitations, we present MarkovCellNet, a distance-driven compositional Markov modeling framework. The approach is modality-agnostic and integrates three key components: (i) distance-sensitive transition kernels that convert pairwise distances (Euclidean, Manhattan, diffusion-based, or learned/kernelized metrics) into movement probabilities between states; (ii) compositional observations represented on the probability simplex, which naturally capture dynamic population structures through counts or proportions across time; and (iii) time-aware lowdimensional embeddings that preserve local geometry and neighborhood relationships prior to kernel construction. From the selected distance metric, we derive stochastic transition matrices, introduce perturbations to account for intrinsic variability (Seneta 2006; Von Luxburg 2007), and propagate population compositions across discrete time steps via Markov updates (Von Luxburg 2007; Zhou et al. 2025). Model adequacy is quantified through a likelihood function on the observed compositions, enabling systematic comparison and selection among competing distance measures, embeddings, and transition schemes. Although our framework applies broadly to any time-evolving population with well-defined features, we focus on single-cell data as our motivating example, since its high dimensionality, compositional structure, and temporal variability provide a natural testbed for both real and synthetic case studies (Luecken and Theis 2019). In this setting, cells that are proximate in appropriate latent embeddings tend to share near-term fates or interactions, consistent with gradual biological change. We construct distance-aware kernels from pairwise relationships in widely used embeddings (t-SNE (Van der Maaten and Hinton 2008), UMAP (McInnes, Healy, and Melville 2018), and PHATE (Poličar and Zupan 2025)) and assess model fit using a multinomial-style likelihood on time-indexed compositions (Van der Maaten and Hinton 2008; Moon et al. 2019). To probe behavior under controlled ground truths, we also design synthetic scenarios that capture common dynamical motifs: divergence/branching, oscillations, and convergence. We then applied the same workflow to insilico, time-course, high-dimensional datasets (simulated single-cell profiles), following identical preprocessing, modeling, and evaluation steps. Using these benchmarks, we show that the distance metric is a first-order modeling choice: diffusion-style distances, when paired with geometry-preserving embeddings often yield higher predictive likelihoods and more coherent longitudinal structures, while alternative metrics emphasize different and sometimes less faithful neighborhood geometries. Taken together, these results position MarkovCellNet as a unified, statistically interpretable approach for analyzing timeevolving populations across modalities, with single-cell dynamics serving as a detailed case study and as a reference implementation that highlights how distance choices drive model behavior and interpretability. 2. Method We aim to statistically characterize state-state interactions in complex populations by analyzing high-dimensional matrices whose rows are entities and columns are modality-specific features (e.g., expression profiles, spatial coordinates, interaction scores) (Poličar and Zupan 2025). We treat these data as multivariate observations over time and model dynamics with a Markov-process view. Standard first-order transition models capture only immediate, one-step dependencies between neighboring states. However, many real systems exhibit mediated effects, where influence propagates indirectly through intermediate states, so restricting to first-order transitions risks missing critical structure. To capture mediated effects beyond immediate neighbors, we construct second-order transition maps that summarize two-step movements and reveal indirect dependencies. Figure 2. MarkovCellNet pipeline: from partitioned features to distance-driven transitions and evaluation. A. Operator construction. Starting with an 𝑁×𝑀 observation-by-feature matrix (𝑁≫ 𝑀), the Cyclic Structural Decomposition Principle (CSDP) reorganizes the data and extracts a sequence of overlapping 𝑀×𝑀 square submatrices (blocks). Within blocks, distances are computed under a chosen metric (e.g., Euclidean, Manhattan, diffusion), converted to affinities, and row-normalized to yield local stochastic operators. The Graph-Augmented Stochastic Transition Modeling (GASTM) step reconciles overlaps and leverages a graph built from those affinities to capture context-specific connections, assembling a global 𝑁×𝑁 transition-frequency/transition matrix that encodes feasible state-to-state movements in the Markov process. B. Model evaluation and metric comparison. Given compositions (e.g., cluster sizes) at an earlier time point (left), MarkovCellNet applies the distance-aware transition operator to predict population shifts (center). Arrows labeled Euclidean, Manhattan, and Diffusion illustrate how different distance metrics induce distinct propagation patterns in the embedded space. He predicted compositions are evaluated against the actual subsequent compositions by utilizing the multinomial cross-entropy (negative log- likelihood) between the observed proportions and the probabilities predicted by the model (with scores aggregated over time; multinomial deviance reported in a similar manner), which facilitates a quantitative ranking of distance metrics and operator configurations, while also illustrating how the geometry of the neighborhood influences the inferred dynamics and their interpretability. We introduce two complementary concepts to build these maps: Cyclic Structural Decomposition Principle (CSDP) and Graph-Augmented Stochastic Transition Modeling (GASTM). CSDP decomposes the learned transition operator into multi-step cyclic components, revealing recurrent temporal motifs (e.g., feedback loops, oscillatory modules, and stable attractors that unfold over more than one step) and providing diagnostics such as cycle strength, persistence, and phase assignments. GASTM extends the Markov framework by embedding transitions in a similarity graph induced by domain-appropriate distance metrics (e.g., Euclidean, Manhattan, diffusion), thereby regularizing estimates, respecting topology, and enabling principled propagation of transition mass to sparsely observe or out-of-block states. Together, CSDP captures how dynamics recur over time, while GASTM ensures transitions are where the geometry says they are, yielding topology-aware, multi-step maps that sharpen interpretability and robustness. Here, edge weights encode local neighborhood structure, while higher-order walks propagate context-specific dependencies. This graph-based augmentation allows transitions to incorporate geometry, density, and anisotropy from the underlying feature space, producing richer stochastic operators that adapt to system topology (Figure 2A). Together, CSDP and GASTM provide a modality-agnostic toolkit for analyzing time-evolving populations. They enable researchers to uncover coordinated processes, progression pathways, indirect couplings, and regime synchrony in heterogeneous systems ranging from single-cell populations to ecological communities, patient cohorts, or user behavior networks. Their integration within our framework ensures that both local transitions and emergent higher-order structures are captured in a unified probabilistic language. We detail their implementation and application in the following subsections. 2.1 Scalable Transition Operator Construction via Matrix Partitioning and Factorization for Stable CSDP/GASTM Analysis In this section, we introduce a block-partitioning framework for large 𝑁×𝑀 observation-byfeature data tables (matrix A, Figure 2A). The table is sliced into overlapping 𝑀×𝑀 square submatrices along both the observation and feature axes. This serves two purposes simultaneously. First, it respects that the dataset contains many observations described by distinct features: local blocks preserve feature-specific patterns, accommodate heterogeneous scales and sparsity, and focus distance calculations on the most informative neighborhoods. Second, it enables enforcement of the stochasticity conditions required for Markov transition modeling. We achieve this by performing blockwise normalization and quality checks, then reconciling overlaps via averaging and renormalization so the assembled operator is globally row-stochastic, numerically stable, and free of boundary artifacts. In practice, partitioning speeds computation by producing smaller kernels, improving cache locality, natural parallelism)(Frigo et al. 1999; Goto and Geijn 2008). It also reduces memory pressure during distance and kernel construction (Malkov and Yashunin 2018), and improves robustness by isolating noisy or under-sampled regions (Ester et al. 1996). Conceptually, it aligns the pipeline with the data’s geometry: local feature neighborhoods yield distance-sensitive transition kernels (Coifman et al. 2005; Von Luxburg 2007), and carefully normalized block-level transitions combine into a coherent, population-level Markov model of dynamics(Kemeny and Snell 1969). This approach integrates cleanly with matrix factorization and out-of-core workflows(Halko, Martinsson, and Tropp 2011) and is the strategy we adopt in MarkovCellNet for scalable, domain-agnostic analysis. Figure 3. Blockwise construction of local transition operators and their Markov interpretation. A. Starting from a large observation-by-feature matrix 𝐴 of size 𝑁×𝑀&&(𝑁≫𝑀&) an 𝑀×𝑀 sliding window is moved down the observation axis to extract overlapping square submatrices 𝐵!,⋯,𝐵"#$%!. Each block contains 𝑀 neighboring observations across all 𝑀 features, providing a locality-preserving view that supports scalable computation and later stitching into a global operator. B. Within each block, pairwise affinities (e.g., from distance metrics such as Euclidean/Manhattan/diffusion) are converted row-wise to a stochastic transition matrix via a SoftMax transformation, yielding non-negative entries whose rows sum to one. The resulting 𝑀×𝑀 matrix defines a discrete Markov process over the 𝑀 local states (schematics at right, with directed edges weighted by the corresponding transition probabilities 𝑃&'). Repeating this procedure across all blocks produces a collection of local operators that can be reconciled to form a coherent population-level Markov model of dynamics. 2.1.1 Cyclic Extraction of Overlapping Submatrices for Localized Analysis To establish the mathematical rigor and structural coherence of our submatrix extraction scheme, we formalize it and prove two results in Appendix 1: (i) an existence theorem showing that, for any 𝑁×𝑀 data matrix with 𝑁≥𝑀, a cyclic sliding window yields a well-defined family of 𝑀×𝑀 overlapping submatrices; and (ii) a topology-preservation theorem showing that this cyclic cover preserves the row ordering and adjacency structure of the original matrix, so local neighborhoods in the observation-by-feature matrix remain local in every extracted block. Together, these theorems justify using blockwise operators that can later be stitched into a coherent global model. We consider an observation-by-feature matrix 𝐴∈ℝ"×$ with 𝑁 rows (observations) and M columns (features) where 𝑀≪𝑁. To analyze dynamics locally while preserving global coherence, we partition 𝐴 into a sequence of overlapping square blocks {𝐵)})*! ", each of size 𝑀×𝑀. For each block index 𝑘∈{1,⋯,𝑁}, we take 𝑀 consecutive rows starting at row 𝑘 and include all 𝑀 columns. When the window extends past row 𝑁, indices wrap around cyclically so that continuity is maintained. Using the row-index mapping 𝑟(𝑖)=8(𝑖−1)&mod&𝑁=+1, the 𝑘-th block is 𝐵)=𝐴[𝑟(𝑘),𝑟(𝑘+1),⋯𝑟(𝑘+𝑀−1);1:𝑀]. Under this scheme, the first row of the first block 𝐵! is “1”, the row orders across blocks are explicit: • 𝐵!:1,2,⋯,𝑀 • 𝐵+:2,3,⋯,𝑀+1 • ⋯ • 𝐵):𝑘,𝑘+1,⋯,𝑘+𝑀−1&(𝐢𝐧𝐝𝐢𝐜𝐞𝐬&𝐭𝐚𝐤𝐞𝐧&𝐦𝐨𝐝𝐮𝐥𝐨&𝑁) • ⋯ • 𝐵":𝑁,1,2,⋯,𝑀−1 This stride-1, overlapping construction has three practical benefits. First, cyclic wrap-around ensures continuity at the matrix boundaries, preserving structural or temporal coherence so that no row is excluded from the analysis. Second, because each row is included in exactly 𝑀consecutive blocks, the design provides balanced participation, giving every observation uniform weight in downstream procedures such as kernel construction and normalization to row-stochastic transition matrices. Finally, the blockwise focus on 𝑀×𝑀 submatrices improve computational efficiency: distances can be computed locally, quality checks and stochastic constraints can be enforced within each block, and overlapping contributions are reconciled through averaging, followed by a global renormalization step that restores row-stochasticity. In implementation, we avoid materializing all 𝐵) in memory; instead, we use views/slices into 𝐴 and iterate over 𝑘, which yields linear memory cost in 𝑀+. While we default to stride 1 (maximal overlap), the procedure generalizes to a stride 𝑠≥1 when faster, coarser scanning is desired; all results reported here use the stride-1 specification above (see Figure 2A). 2.1.2. Constructing Probability Transition Matrix for Stochastic Population Dynamics Using our block-partitioning and factorization procedure, we construct a probability transition matrix that encodes stochastic movements among states using explicitly chosen distance metrics (e.g., Euclidean, Manhattan, diffusion) computed within each block. This probability transition matrix characterizes the movement behavior of each entity (Moon et al. 2019), with probabilities assigned to two possible outcomes: (1) remaining in its current location, or (2) transitioning to the position of a neighboring target entity (Haghverdi et al. 2016). Blockwise normalization enforces valid probabilities (nonnegative probabilities that sum to one), while overlaps are reconciled by averaging and renormalizing to yield a stable, global operator. The result is a compact, modalityagnostic Markov model of population dynamics that scales efficiently and treats the choice of distance metric as a core, testable modeling decision. 2.1.2.1. Assumptions Underlying the Probability Transition Matrix Before specifying the operator, we make explicit the modeling assumptions that render the construction identifiable, interpretable, and computationally tractable: 1. Discrete State Space: Entities move among a predefined set of discrete states or locations (e.g., grid cells, graph nodes, or feature-defined clusters). This reflects structured environments, improves interpretability, and keeps computation tractable. 2. Distance-Dependent Transitions: Movement propensities depend on proximity under a chosen distance metric suited to the modality (e.g., Euclidean, Manhattan, diffusion, or learned/graph distances). The metric choice is a first-class modeling decision and may include self-transitions to capture persistence. 3. Softmax-Based Probabilities: Distances are converted into probabilities via a Softmax function (Bridle 1989). This ensures nonnegativity and row-wise normalization, while smoothly biasing toward closer states. A temperature-like parameter controls how sharply the model prefers nearest neighbors. These assumptions yield a coherent, scalable population-level Markov model of dynamics that generalizes across domains 2.1.2.2. Compute the Distance Matrix To construct transition probability matrices that capture system dynamics, we begin by computing pairwise distances between entities in an 𝑀-dimensional feature space, derived from each squared matrix 𝐵), where rows correspond to the coordinates of individual entities (Schiebinger et al. 2019). We consider three commonly used distinct distance metrics including Euclidean, Manhattan, and diffusion distance which are widely applied in stochastic modeling and trajectory inference and serve as the foundation for defining similarity-based transition kernels. For each block 𝐵) (row = entities, columns = 𝑀 features), we compute a finite family of distance metrics indexed by 𝑟∈{1,⋯,𝑞}. In this study 𝑞=3 where 𝑟=1 Euclidean, 𝑟=2 Manhattan, and 𝑟=3 diffusion. We keep a single, consistent notation 𝑑),&' (.) for the distance between entities 𝑖 and 𝑗 within block 𝐵) under metric 𝑟. - Euclidean (𝑟=1). 𝑑),&' (!)=WX8𝐵),&,0 −𝐵),',0=+ $ 0*! , the straight-line proximity in feature space (intuitive, fast; can be noise-sensitive in high𝑀). - Euclidean (𝑟=2). 𝑑),&' (+) =XY𝐵),&,0 −𝐵),',0Y $ 0*! , aggregating absolute deviations and often more robust to outliers. - Diffusion (𝑟=3). We first build a local random-walk operator 𝑃) (1) on the block graph (e.g., by converting a base distance to similarities and row-normalizing). The diffusion distance at walk length 𝑡 compares probability flows: 𝑑),&' (1) =[X\𝑝2 (1)(𝑖,𝑙)−𝑝2 (1)(𝑗,𝑙)_+ 𝜋(1)(𝑙) 3a! + ⁄, where 𝑝2 (1)(𝑖,𝑙) is the 𝑡-step transition probability and 𝜋(1) the stationary distribution of 𝑃) (1). A common, efficient approximation uses the spectral decomposition of 𝑃) (1) : retaining the top 𝐿 eigenpairs (𝜆3,𝜓3), 𝑑),&' (1) ≈fX𝜆3 +2 5 3*! (𝜓3(𝑖)−𝜓3(𝑗))+g! + ⁄ which preserves long-range, manifold-aware similarities. For each metric 𝑟 and block 𝑘, we assemble the pairwise distance matrix 𝐷) (.) =i𝑑),&' (.) j, convert distances to similarities (e.g., a monotone kernel 𝜅(𝑑), and row-normalize to obtain a distancedependent Markov transition matrix 𝑃) (.) (self-transition mass encodes persistence; light regularization mitigates sparsity). This yields a family of operators l𝑃) (.)m)*! " for each candidate metric 𝑟. Finally, although we analyze each metric individually, we select the best-performing metric by conditional maximization of the data likelihood under the corresponding distance-dependent Markov process. Concretely, letting ℒ(𝑟) denote the aggregated multinomial (cross-entropy) likelihood of the observed next-time compositions given the transitions l𝑃) (.)m, we choose 𝑟∗=arg max .∈{!,⋯,:}ℒ(𝑟), and carry 𝑟∗forward for reporting and downstream interpretation. This formulation makes metric choice a first-class, testable modeling decision, while keeping notation and estimation uniform across metrics. 2.1.2.3. Compute the Transition Probability Matrix To statistically characterize the temporal evolution of population states, we define a family of blockwise, finite-state Markov chains indexed by block 𝑘. For each overlapping block 𝐵), the local system is represented by a right-stochastic transition matrix 𝑃)∈ℝ$×$, where 𝑀 is the number of local states within that block’s feature or embedded space. Each entry 𝑃)(𝑖→𝑗) is the conditional probability of moving from state 𝑖 to state 𝑗 within the block. We construct 𝑃) by computing pairwise distances among states inside the block under a chosen distance metric and neighborhood rule, transforming those distances into nonnegative affinities, and then normalizing rows to enforce stochasticity; a self-transition term captures persistence, and light regularization mitigates sparsity and noise. The collection {𝑃)} is subsequently reconciled across overlapping blocks to yield a stable global operator. This formulation is modality-agnostic and applies to any setting where entities evolve through locally plausible transitions in feature space. Let 𝑑),&' (.) denote the distance between states 𝑖 and 𝑗 at step 𝑘, computed using one of several candidate distance metrics: Euclidean (E), Manhattan (M), or Diffusion (D). To model the dependency of transition probabilities on state proximity, we employ a distance-aware probabilistic kernel defined via a scaled softmax transformation: 𝑃)\𝑖→𝑗t𝑑),&' (.),𝐷) (2<=>)_=𝑒𝑥𝑝8−𝛽.𝑑),&' (.)= ∑𝑒𝑥𝑝\−𝛽.𝑑),&' (.)_ $ 3, where 𝑃)(𝑖→𝑗|𝐷) (.)) is the conditional probability of transitioning from state 𝑖 to state 𝑗 at the k-th step, given the distance matrix 𝐷) (.) (Figure 3B). 𝑑),&' (.) is the distance between states 𝑖 and 𝑗 based 𝑝8Γ&YΨ&,𝐷(.)== 𝑛&! ∏𝛾&'! " '*! ¡𝜓&' H"# " '*! . Step 3: Posterior Inference via Dirichlet-Multinomial Conjugacy Due to conjugacy, the posterior distribution for Ψ& is also Dirichlet: Ψ&|Γ&,𝐷(2<=>)~Dirichlet(𝛼&! +𝛾&!,𝛼&+ +𝛾&+,⋯,𝛼&" +𝛾&") yielding the joint posterior distribution: Ψ|Γ,𝐷(.)~¡Dirichlet(𝛼&+Γ&) " &*! This Bayesian framework allows us to capture the full posterior uncertainty over transition probabilities while remaining agnostic to the specific choice of distance metric. Additionally, normalized or scaled versions of Γ can be accommodated by interpreting them as effective counts, thus enabling robust inference across varying data granularities and similarity measures. 2.4. Modeling Population Composition Dynamics We now move from local, metric-aware state transitions to global composition dynamics. After partitioning the observation-by-feature table for scalability, computing distance metrics within blocks, and assembling a reconciled, row-stochastic transition operator that encodes feasible stateto-state movements, we link this operator to the observed population counts over time. In this subsection, we represent these counts as compositions on the probability simplex and specify a probabilistic evolution that advances compositions across time using the metric-informed transition operator. The likelihood formulation naturally accounts for sampling variability and enables principled comparisons of alternative metrics and embeddings. This closes the pipeline from distances → kernels → transition operator → composition updates, yielding a coherent, scalable population-level Markov model of dynamics that we instantiate in MarkovCellNet. Let Π∈ℝEF "×I denote the observed count matrix, where column Π(2) =(𝜋!2,𝜋+2,⋯,𝜋"2)I contains counts of 𝑁 subpopulations at time 𝑡. We normalize each column to obtain empirical compositions Π¨(2) =8𝜋!2 𝑁(2) ⁄,⋯,𝜋"2 𝑁(2) ⁄=I∈∆"#!,&&&&𝑁(2)=X𝜋&2 " &*! , where ∆"#! is the probability simplex. Observed counts at 𝑡+1 follow 𝑋(2%!)~Multinomial&8𝑁(2%!),𝜇(2%!)=. with mean composition 𝜇(2%!)∈∆"#! generated from the previous composition Π¨(2) by a distanceaware transition operator plus noise that is guaranteed to remain on the simplex. Building on this operator, we structure the update in three layers that separate signal from uncertainty while guaranteeing compositions stay on the simplex. First, a deterministic propagation step applies the distance-derived transition kernel (with optional constraints), followed by elementwise clipping and renormalization, yielding the baseline mean-field trajectory. Second, to capture stochastic variability while preserving positivity and unit-sum constraints, we add noise on the simplex by perturbing the baseline in a compositionally coherent space (e.g., log-ratio) and then map back to the simplex, with a few hyperparameters controlling perturbation size. Third, as a robustness check, we replace this noise model with a simplex-native distribution (e.g., Dirichlet or logistic-normal), ensuring conclusions are not tied to a single parameterization. Together, these three layers: deterministic propagation, simplex-respecting stochasticity, and a distributional robustness via alternative perturbations yield an update rule that is both principled and practical for evaluating metrics, tuning locality, and scoring out-of-sample fit. 2.4.1. Baseline (deterministic) update From the previous empirical composition Π¨(2) ∈∆"#!, we compute the mean-field (noise-free) next-step composition by propagating through the distance-derived transition kernel and the optional constraint operator, then stabilizing and renormalizing: 𝜇¬(2%!)=𝐹8Π¨(2);𝒦8𝐷(.)=,Ψ2,𝛿==ℛ\clipJ8Ψ2 I𝒦8𝐷(.)=Π¨(2)=_. Here 𝒦8𝐷(.)=∈ℝ"×" is the transition/kernel operator built from metric 𝑟 (e.g., Euclidean, Manhattan, diffusion) on the chosen embedding/feature space; Ψ2∈𝑅"∗" is a (possibly timevarying) stochastic operator that encodes modeling constraints or sampling adjustments (e.g., conservation constraints, exposure offsets, or controlled re-weighting); clipJ(𝑥) applies an elementwise floor 𝛿>0 to prevent zeros/negatives with clipJ(𝑥)K=max8𝑥K,𝛿= with 𝛿 small (e.g., min&(10#L, 0.01/𝑁)) to avoid numerical issues (e.g., log0) while having negligible effect on the fitted dynamics; and ℛ(𝑥)=𝑥𝟏I𝑥 ⁄ renormalizes to the simplex. Clipping is immediately followed by renormalization so 𝜇¬(2%!)∈∆"#!. 2.4.2. Noise on the simplex (ALR-Gaussian perturbation) To introduce stochastic variability while preserving positivity and unit-sum constraints, we perturb the baseline in an additive log-ratio (alr) coordinate system, add Gaussian noise and map back to the simplex: 𝑧(2%!) =alr8𝜇¬(2%!)==&alr\𝐹8Π¨(2);𝒦8𝐷(.)=,Ψ2,𝛿=_&∈𝑅"#!, 𝑧MN&O> (2%!) =𝑧(2%!) +𝜀,&&&&&&&&&𝜀~𝒩(𝑚2,Σ2), 𝜇(2%!)=alr#!8𝑧MN&O> (2%!)==softmax8i𝑧MN&O> (2%!),0j=∈∆"#!. We make the dependence of the noise explicit by allowing 𝑚2=𝑚8Π¨(2),𝒦8𝐷(.)=,ΨI,𝛿= (typically 𝑚2=0) and Σ2=∑8𝜇¬(2%!),𝑁(2%!),𝒦8𝐷(.)=,ΨI,𝛿= (e.g., 𝜎2 +𝐼 or diagonal with 𝜎2 + learned by marginal likelihood or cross-validation, or scaled as 𝑐𝜇¬(2%!) ⁄). Here, the perturbation is applied to 𝜇¯(2%!) produced by 𝐹(⋅), and the noise hyperparameters can adapt to the kernel, constraint, sample size, and stabilization choices used to generate that baseline. 2.4.3. Robustness alternative (Dirichlet perturbation) As a simplex-native sensitivity check, we also consider a Dirichlet perturbation centered at 𝜇¬(2%!) 𝜇(2%!)~𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡8𝜅2𝜇¬(2%!)=, with 𝜅2=𝜅8𝜇¬(2%!),𝑁(2%!),𝒦8𝐷(.)=,Ψ2,𝛿= The concentration 𝜅2>0 is a scalar, with small 𝜅2 producing wide variability around 𝜇¯(2%!), while large 𝜅2 concentrating mass near the baseline. 𝜅2 can be chosen to reflect sampling depth 𝑁(2%!), a target coefficient of variation, or selected by out-of-sample likelihood. For scoring and model selection (e.g., comparing metrics 𝑟), we compute the multinomial negative log-likelihood of 𝑋(2%!) under 𝜇(2%!)aggregated over time, and we renormalize any intermediate vector so that every update remains in Δ"#!. 2.5. Evaluation of Model Fit via Multinomial Likelihood We assess model fit by scoring time-ahead predictions against observed counts using a multinomial (or Dirichlet-multinomial) likelihood, consistent with the transition operator used in construction. Recall that the deterministic mean update 𝜇¯(2%!) is produced by 𝜇¯(2%!) =𝐹8Π¨(2);𝒦8𝐷(.)=,Ψ2,𝛿=, and the stochastic prediction incorporates simplex-respecting noise: 𝜇(2%!) =𝐺8𝜇¯(2%!);𝜀=. Specifically 𝜇(2%!) =𝐺8Π¨(2);𝒦8𝐷(.)=,Ψ2,𝛿,𝜀=. Here, 𝐷(.) encodes the distance metric and its neighborhood/kernel hyperparameters (e.g., 𝑘, bandwidth), 𝒦(⋅) is the corresponding transition kernel assembled block-wise and reconciled across overlaps, 𝛹2 denotes any time-varying constraint/sampling operator (e.g., exposure offsets, conservation constraints, re-weighting), and 𝛿 is the stabilization floor followed by renormalization. Let 𝑋(2%!) =8𝑥! (2%!),…,𝑥" (2%!)&= denote observed counts over 𝑁 states, with total 𝑁(2%!) = X𝑥& (2%!) & . Conditional on 𝑁(2%!), we score the prediction 𝜇(2%!) with the multinomial likelihood ℒMult (2%!) =𝑃8𝑋(2%!) ∣𝑁(2%!),𝜇(2%!)== 𝑁(2%!)! ¡𝑥& (2%!) &!¡8𝜇& (2%!)=P" (%&') " &*! . In the presence of extra-multinomial variation (over-dispersion), we use the Dirichlet-multinomial marginal likelihood instead with 𝜇(2%!) ∼Dirichlet8𝜅2 𝜇¯(2%!)=,& 𝑋(2%!)&~Dirichlet−Multinomial&(𝑁(2%!),𝜅2 𝜇¯(2%!)). The model objective is the aggregated negative log-likelihood over time 𝒥(𝜃)=−Xlog&ℒ(2%!)(𝜃) I#! 2*F where 𝐿(2%!)is chosen as multinomial or Dirichlet-multinomial depending on dispersion diagnostics and 𝜃={ 𝑟, 𝑘, bandwidth, Ψ⋅, 𝛿, 𝜺, block size/overlap, regularization }. We evaluate 𝒥(𝜃) in held-out time folds (or rolling-origin splits) to obtain an out-of-sample score, and we report complementary summaries such as multinomial deviance, calibration checks (probability integral transform on held-out counts, coverage of predictive intervals), and residual diagnostics on per-state log-odds. Because the likelihood depends directly on 𝒦(𝐷(.)), constraints Ψ2, stabilization 𝛿, and the noise layer, improvements or degradations in 𝒥(𝜃) can be attributed to specific design choices. In practice, higher out-of-sample likelihoods indicate that the distance-aware operator and its hyperparameters offer a coherent explanation for the observed state shifts, whereas lower likelihoods flag misspecification (e.g., overly global kernels creating short-circuits, insufficient locality, unmodeled transitions, sampling bias, or under/over-dispersion). 3. Results In this work, we apply MarkovCellNet to model the temporal evolution of complex single-cell populations, where each cell is represented by a high-dimensional gene expression vector(Van Kampen 1992). As outlined in Figure 6, the workflow proceeds from a high-dimensional gene expression matrix to a geometry-preserving 2D embedding (Becht et al. 2019; Luecken and Theis 2019) and then to a neighborhood map defined by distances between cells. Pairwise proximities are computed using alternative distance metrics (Euclidean, Manhattan, diffusion), which define local neighborhoods and are converted into distance-driven transition kernels(Coifman et al. 2005). These kernels allocate probability mass to either remaining in place or moving to nearby states(Metzner, Schütte, and Vanden-Eijnden 2009). We then apply MarkovCellNet to propagate population compositions through time, assessing model fit with a multinomial likelihood(Grün and van Oudenaarden 2015). Figure 6. MarkovCellNet workflow for geometry-aware cell-cell interaction modeling. Starting from a high-dimensional gene-by-cell expression matrix (left), cells are embedded into a lowdimensional manifold (middle; example 2-D layout). On this embedding, MarkovCellNet computes pairwise distances using alternative metrics, (1) Euclidean, (2) Manhattan, and (3) diffusion, and converts these distances into transition probabilities to define a stochastic (Markov) operator over cells (right; arrow thickness reflects transition weight). Each metric is evaluated for preservation of biological programs and for out-of-sample predictive fit; the best-scoring metric is selected and used to model cell-cell interactions for downstream analyses. 3.1. Synthetic dynamics for benchmarking (linear, oscillatory, and convergent processes) To illustrate performance under distinct geometries, we stress-test distance choices using three archetypal dynamic regimes that recur in biology and population systems. First, divergence/branching captures progressive separation of states, as in lineage specification or treatment-induced phenotypic drift. Second, periodic/oscillatory behavior reflects rhythmic programs (e.g., cell-cycle waves, circadian regulation) in which states recur with phase. Third, convergence/relaxation models stabilization toward a common attractor, consistent with recovery after perturbation or homeostatic re-equilibration. These well-studied motifs allow us to tie simulation parameters to interpretable processes while retaining full control over ground truth for benchmarking. To instantiate these motifs, we generated in-silico, time-course, high-dimensional datasets in which cluster centroids follow explicit, mathematically specified trajectories while covariances are held fixed to control noise. In the Linearly Divergent States (LD-5) dataset, five centroids move at constant velocity over time (a standard linear transport), producing progressively increasing separation. In the Periodically Oscillating States (PO-4) dataset, four centroids follow sinusoidal paths with prescribed amplitudes, periods, and phases, mimicking phase-advancing programs. In the Logarithmically (exponentially) Convergent States (LC-3) dataset, three centroids relax toward a shared attractor with first-order decay, yielding coordinated convergence. We generated 5,000, 4,000, and 3,000 samples with 50 features for LD-5, PO-4, and LC-3, respectively, and retained the time index to support time-aware evaluation. This design captures linear drift, sinusoidal oscillation, and exponential relaxation, canonical forms with well-defined dynamics while remaining simple enough to isolate the impact of metric choice and transition construction. To make our simulations fully reproducible, we formalize the data-generating mechanism used to create the in-silico time-course datasets. We generate a gene-by-cell matrix 𝑋∈ℝ=×" with 𝑝=50 genes (indexed by 𝑔) and 𝑁cells/samples (indexed by 𝑐). Each cell is assigned a latent state (cluster) 𝑘∈{1,…,𝐾} and a continuous time 𝑡@∈[0,1]. Conditional on (𝑘,𝑡@), expression is drawn from a Gaussian with stateand time-dependent mean and fixed noise covariance, Here 𝑚)(𝑡)∈ℝ=is the centroid trajectory for state 𝑘. Unless stated otherwise, we sample 𝑘 uniformly and 𝑡@∼Uniform(0,1). We use the standard basis vectors 𝑒' in ℝ= to define directions of motion; coordinates not referenced in the formulas remain mean-zero noise. - Linearly Divergent States (LD-5) To model linear drift with increasing separation, we set 𝐾=5 states and 𝑁=5,000&cells. Centroids start at a common origin and move at constant velocity along orthogonal axes: 𝑚)(𝑡) = 𝑣 𝑡 𝑒),&&&𝑣=3.0,&&&&&&𝑘=1,⋯,5,&&&&&&&&&𝑡∈[0,1]. Thus, the distance between any two state centroids grows linearly with 𝑡; at 𝑡=1 the pairwise separation is 𝑣√2=3√2. - Periodically Oscillating States (PO-4) To capture cyclic programs with phase offsets, we set 𝐾=4 and 𝑁=4,000. Each centroid traces a circle of radius 𝑅 in a distinct two-dimensional subspace with a fixed phase shift: 𝑚)(𝑡) = 𝑅[cos(2𝜋𝑡+𝜙))+sin(2𝜋𝑡+𝜙))],&&&&&&𝑅=2.5, 𝜙)∈&Ê0,𝜋2,𝜋,3𝜋 2Ë All states complete one full oscillation over 𝑡∈[0,1], producing phase-advanced trajectories. - Logarithmically (Exponential) Convergent States (LC-3) To model relaxation toward a common attractor, we set 𝐾=3 and 𝑁=3,000. Centroids begin at distance 𝑠 along orthogonal axes and decay exponentially toward the origin: 𝑚)(𝑡) = 𝑠 𝑒#R2 𝑒),&&&&&&&&&&𝑠=3.0, 𝜆=2.5,&&&&&&&&&𝑘=1,2,3. These yields coordinated convergence with a first-order decay rate 𝜆; by 𝑡=1 the amplitude is reduced by 𝑒#R. Across all motifs we retain the true time 𝑡@ for evaluation, keep Σ fixed to control noise, and balance state proportions to avoid confounding difficulty with prevalence. These explicit centroid functions, parameter values (𝑣,𝑅,𝜙),𝑠,𝜆,𝜎), and sample sizes (𝑁) enable exact reproduction of LD-5, PO-4, and LC-3. Crucially, MarkovCellNet does not require discrete labels or low-dimensional embeddings. The transition operator is trained directly on continuous features or embeddings using neighborhood structure alone. Labels, when present, are used only for benchmarking (e.g., reporting agreement with known regimes) or post-hoc interpretation. In real data, where labels are rarely available, approximate groupings can be obtained via a standard unsupervised pipeline: dimensionality reduction for denoising, k-nearest-neighbor graph construction, community detection (e.g., Leiden) with resolution sweeps and bootstrap stability checks, optional density-based refinement (e.g., HDBSCAN), and conservative label propagation across time to reduce fragmentation. This preserves the unsupervised nature of the framework while providing interpretable summaries for evaluation. To enhance interpretability, we refer to the simulated datasets by their characteristic temporal behaviors: Linearly Divergent States capture stepwise differentiation, Periodically Oscillating States reflect rhythmic dynamics, and Logarithmically Convergent States represent gradual unification in synchronized systems (Bishop and Nasrabadi 2006). This nomenclature links simulated cluster trajectories to well-recognized dynamical patterns, providing an intuitive and biologically interpretable framework for evaluating how different distance metrics and transition operators capture temporal population structure. 3.2. Dimensionality Reduction to Preserve Temporal and Structural Features We used dimensionality reduction methods (DRMs) to project high-dimensional, time-stamped measurements into lower-dimensional latent spaces while retaining neighborhood structure relevant for transition modeling. We selected three complementary DRMs including t-SNE, UMAP, and PHATE because they emphasize different aspects of geometry: t-SNE prioritizes local neighborhoods (useful for resolving fine-grained, transient states), UMAP balances local and global relationships (often preserving inter-cluster topology), and PHATE preserves diffusion geometry, which can reveal smooth progressions along putative trajectories(Van der Maaten and Hinton 2008; McInnes, Healy, and Melville 2018; Poličar and Zupan 2025). Applied independently to our synthetic time courses, all three methods recovered the intended structures: t-SNE produced compact neighborhoods but occasionally fragmented long paths; UMAP yielded smoother arrangements capturing both cluster separation and between-state bridges; PHATE most clearly delineated continuous progressions and intermediate states. We acknowledge that our three canonical motifs (linear drift, sinusoidal oscillation, exponential relaxation) emphasize continuous evolution. To mitigate bias, we include additional stress tests and real-data ablations: piecewise-constant dynamics with abrupt jumps, disconnected mixtures with intermittent bridging, and switching regimes with sparse observation windows. In these settings, likelihood-based selection typically favors Euclidean/Manhattan over diffusion, and our diagnostics reflect the lack of a reliable manifold. We discuss these outcomes and practical guidance in the Discussion, including when diffusion geometry is beneficial versus when ambient metrics are preferable(Coifman and Lafon 2006; McInnes, Healy, and Melville 2018). DRMs provide useful summaries, but we do not assume a manifold by fiat; instead, we test the assumption, compete metrics under a common likelihood, and adapt kernel construction accordingly. This safeguards MarkovCellNet against violations of continuity while retaining the advantages of geometry-aware modeling when the assumption is supported. To characterize temporal dynamics in high-dimensional single-cell time-series data, we first applied dimensionality reduction methods (DRMs) to project the data into lower-dimensional latent spaces while preserving key structural features. We selected three widely used DRMs: t-SNE, UMAP, and PHATE based on their complementary strengths in capturing local versus global structures. t-SNE emphasizes local neighborhood preservation, enhancing detection of subtle or transient cell-state changes. UMAP balances local and global structure, capturing both discrete clusters and their intercluster relationships. PHATE preserves diffusion-based distances, effectively revealing continuous trajectories in biological processes (Van der Maaten and Hinton 2008; McInnes, Healy, and Melville 2018; Poličar and Zupan 2025). We applied each DRM independently to our synthetic datasets, designed with five, four, and three clusters respectively (Figure 7). All methods recovered the expected cluster structure, though interpretability varied. t-SNE generated compact, locally coherent clusters but sometimes fragmented trajectories. UMAP produced smoother, globally coherent embeddings, preserving cluster separation while reflecting temporal transitions. PHATE excelled at continuous trajectories, revealing intermediate states along smooth developmental paths (Moon et al. 2017; Zhang, Song, and Song 2023). By combining these complementary DRMs, we captured local, global, and continuous temporal structures, providing a robust foundation for downstream Markov-based transition modeling and ensuring latent spaces retained sufficient biological and statistical signal. Figure 7. Temporal embeddings of a synthetic evolving population: PHATE vs. t-SNE vs. UMAP. Rows correspond to three time points (Time 1-3); columns show 2D embeddings produced by PHATE, t-SNE, and UMAP from the same high-dimensional data at each time. Colors indicate ground-truth subpopulations/states (labeled “Cell 1–5” in the panels). PHATE (left) preserves a smooth, nearly circular manifold across times, with colors changing gradually along the continuum, evidence of strong global continuity and a clear temporal progression. t-SNE (middle) emphasizes local neighborhoods, fragmenting the manifold into separated islands and distorting between-time relationships, which can yield diffuse or unstable transition neighborhoods. UMAP (right) recovers a dominant one-dimensional trajectory with better global coherence than t-SNE but still compresses/warps the ends relative to PHATE. These geometric differences directly affect distancesensitive transition kernels in MarkovCellNet, and thus the inferred dynamics and likelihood scores: embeddings that preserve global structure support smoother, more interpretable state-to-state transitions over time. 3.3. Distance-Aware Markov Modeling of Reduced Cellular Dynamics We generated synthetic single-cell datasets with known dynamics (linear, branching, and cyclic), added noise, nuisance genes, sparsity, disconnected populations, and batch effects, then normalized and defined a ground-truth transition reference for scoring. We embedded each dataset with UMAP, t-SNE, and PHATE, computed Euclidean, Manhattan, and diffusion distances, built local transition kernels in MarkovCellNet, and evaluated all embedding, metric pairs by held-out transition likelihood, edge precision, recall, pseudotime concordance, and topology diagnostics. In well- sampled, continuous settings, diffusion distances on PHATE or UMAP performed best; under sparsity or disconnected mixtures, ambient metrics (Euclidean or Manhattan on UMAP) outperformed diffusion, while t-SNE preserved very local neighborhoods but weakened long-range ordering. These results support a simple rule aligned with our objectives: let held-out likelihood choose the operator, use diffusion when continuity is strong, and switch to ambient distances with tighter locality when diagnostics flag violations to more reliably reconstruct cellular dynamics. 3.3.1. Evaluating Cell Interactions Across Distance Metrics via Log-Likelihood We characterized stochastic cell behavior within reduced-dimensional manifolds by constructing Markov transition probability matrices that capture transitions between discrete cell states. These matrices were derived in latent spaces using three distance metrics: Euclidean, Manhattan, and diffusion distances, the latter incorporating structural continuity via graph-based diffusion. We assessed model quality using log-likelihood, quantifying how well each transition matrix reflected observed cell-state transitions. Higher log-likelihoods indicate stronger agreement between modeled and empirical dynamics, allowing direct comparison of distance metrics in combination with each DRM. We observed that diffusion-based matrices achieved the highest log-likelihoods when applied to PHATE embeddings, consistent with PHATE’s design for capturing continuous and branched trajectories. Across all datasets and time points, diffusion distance consistently outperformed Euclidean and Manhattan metrics, highlighting its ability to preserve both local and global relationships in the reduced space. Combinations like Euclidean + t-SNE and Manhattan + t-SNE performed well, but diffusion + PHATE emerged as the most robust pairing (Figure 8A-D). These results demonstrate that the choice of distance metric interacts with embedding geometry to shape transition modeling. Our framework provides a quantitative basis for selecting compatible DRM-distance metric pairs, enabling accurate inference of dynamic cell-state transitions and guiding the analysis of temporal single-cell datasets. Figure 8. Likelihood-based comparison of distance metrics and embeddings across time. Panels show multinomial log-likelihood scores (higher/less negative = better fit) for MarkovCellNet’s distance-aware transition operators evaluated on three time points. Colors within panels A-C indicate the embedding used (PHATE, t-SNE, UMAP). A. Euclidean distance: performance varies with embedding; PHATE/UMAP generally exceed t-SNE. B. Manhattan distance: similar trend, with t-SNE typically the weakest and PHATE/UMAP yielding higher likelihoods. C. Diffusion distance: consistently the best-performing metric across time, especially when paired with geometry-preserving embeddings (PHATE/UMAP). D. Aggregate comparison (mean ± SD across embeddings) ranks metrics Diffusion > Manhattan > Euclidean, highlighting that the choice of distance metric, and its interaction with the embedding materially shapes model fit. Error bars denote variability over repeated runs; all panels use the same vertical scale for comparability. 3.3.2. Heatmap-Based Evaluation of Distance-Aware Cell Interaction Models We further examined how dimensionality reduction and distance metrics shape modeled cell interactions by visualizing transition probability matrices as heatmaps. Each heatmap depicts the probability of a cell of type i transitioning to type j in the reduced space, providing an intuitive view of directional interaction strengths across cell populations. Figure 9. Transition structure depends on the pairing of embedding and distance metric. Heatmaps show one-step transition probability matrices estimated by MarkovCellNet for a five-state system at a representative time point. Rows of each heatmap are source states and columns are target states (labels “Cell 1–5” denote state IDs); each row sums to one, so a uniform baseline is ≈0.20. The rows of the panel correspond to the embedding used to construct neighborhoods (PHATE, t-SNE, UMAP); the columns correspond to the distance metric used to convert pairwise distances into transition probabilities (Euclidean, Manhattan, Diffusion). Dendrograms indicate hierarchical similarity among rows/columns. With Euclidean or Manhattan distances on PHATE/UMAP, matrices are nearly uniform, indicating weak directional information. t-SNE often yields scattered, less stable patterns due to its strong local emphasis. In contrast, Diffusion distance produces pronounced block-diagonal structure and selective off-diagonal flows, most evident with PHATE/UMAP reflecting manifold connectivity and coherent progression. The figure highlights that metric-embedding synergy controls the sparsity, anisotropy, and interpretability of the transition operator, which ultimately governs dynamic predictions in MarkovCellNet. These visualizations allowed us to assess specificity and clarity: red regions indicate highprobability transitions, while blue regions reflect weak or infrequent interactions. By comparing heatmaps across all nine model configurations (three DRMs × three distance metrics) at Time 1, we found that the best-performing combinations particularly PHATE or UMAP paired with diffusion distance produced sharp, localized high-probability regions, such as distinct diagonals or focused off-diagonal clusters, reflecting coherent and biologically meaningful interaction dynamics (Figure 9). In contrast, t-SNE combined with Euclidean or Manhattan distances often yielded diffuse, less interpretable patterns. Notably, PHATE paired with diffusion distance produced the clearest, most distinct transitions. Across all DRMs, diffusion distance consistently enhanced interpretability, aligning with our loglikelihood analyses. These results reinforce diffusion distance, particularly when combined with PHATE or UMAP, as a robust metric for modeling temporally coherent and biologically meaningful cell interactions. 4. Discussion We set out to address the broader challenge of characterizing temporal evolution in complex populations, where state changes can reflect development, adaptation, rhythmic cycles, synchronization, regime shifts, or recovery. To this end, we developed MarkovCellNet, a modalityagnostic framework that models state-state interactions via distance-aware Markov transition operators defined on low-dimensional embeddings and evaluated on population compositions over time. While we focused on single-cell transcriptomics as a concrete example, the modeling pipeline, distance metrics → kernels → transition operator → composition updates, extends naturally to other domains (e.g., ecological communities, patient cohorts, online user populations, financial portfolios, transportation flows). Many DRMs implicitly assume that the high-dimensional data concentrate on a low-dimensional, continuous manifold. Real datasets (including single-cell experiments) can deviate from this ideal: mixtures of disconnected populations, abrupt state switches, or sparse sampling can break continuity. Our distance-driven framework is designed to be robust to such violations in three ways. First, we evaluate multiple metric families in parallel (Euclidean, Manhattan, diffusion) and select the metric that maximizes out-of-sample likelihood of observed compositions; if continuity is weak, ambient-space metrics (Euclidean/Manhattan) often win over diffusion. Second, we run diagnostics that do not assume a manifold: trustworthiness/continuity scores for neighborhood preservation, kNN graph connectivity (components, short-circuit edges), and local intrinsic dimensionality estimates to flag regions where a low-dimensional embedding is implausible(Venna and Kaski 2006; Amsaleg et al. 2019). Third, our block-wise construction limits spurious long-range links: distances and kernels are computed locally and reconciled across overlaps, which reduces the chance that a global embedding artifact induces erroneous transitions. When diagnostics indicate manifold violations, we (i) down-weight diffusion-based kernels, (ii) increase locality (smaller neighborhood sizes, adaptive bandwidths), or (iii) operate without a nonlinear embedding (e.g., on PCA-whitened features), and let the likelihood-based model selection pick the safer operator. A central contribution is to elevate distance metrics to explicit modeling choices rather than neutral preprocessing. Different metrics emphasize different structures, short-range isotropic proximity (Euclidean), axis-aligned separations (Manhattan), or manifold connectivity and multi-step coherence (diffusion), and therefore induce different transition neighborhoods, interaction strengths, and inferred trajectories. Across our synthetic benchmarks, diffusion distances paired with