Pith. sign in

REVIEW 4 major objections 5 minor 13 references

Graph Pseudotime Analysis and Neural Stochastic Differential Equations for Analyzing Retinal Degeneration Dynamics and Beyond

T0 review · 4 major / 5 minor · reviewed 2026-08-08 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read The paper claims that graph-level pseudotime analysis combined with neural stochastic differential equations places an irreversible bifurcation at step 4 of the estimated JR5558 retinal-degeneration trajectory, with pathway stability…

desk verdict The paper's headline bifurcation at Step 4 is a displacement-threshold artifact, and the pseudotime ordering is partly built from the severity label—but the new dataset and the graph-pseudotime idea are worth a serious look. read the letter →

arxiv 2502.06126 v1 pith:RG7TJU7G submitted 2025-02-10 cs.LG

Dai Shi , Kuan Yan , Lequan Lin , Yue Zeng , Ting Zhang , Dmytro Matsypura , Mark C. Gillies , Ling Zhu
show 1 more author
Junbin Gao
This is my paper · ORCID
classification cs.LG
keywords pseudotimeanalysisgraphneuralnetworksstochasticdifferentialequationsretinaldegenerationdiseasebifurcationpathwaystabilitytranscriptomicspointofnoreturn
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper tries to turn static snapshots of retinal degeneration into a time-ordered molecular trajectory. It builds one pathway graph per mouse, orders the graphs with graph-level pseudotime analysis, and fits a neural stochastic differential equation to the resulting trajectory. From that fit it defines pathway stability and formal bifurcation conditions, and it reports a system-wide point of no return at step 4 in the JR5558 dataset. If this is right, the disease has an early therapeutic window, and the stable versus non-stable pathway rankings offer concrete intervention candidates.

What carries the argument

The machinery is a staged pipeline. Each mouse becomes a weighted graph whose nodes are canonical molecular pathways defined by gene-set annotations and whose edges connect pathways sharing genes, with edge weights computed from feature distances. Graph-level Pseudotime Analysis embeds these graphs through spectral positional encoding, pooling, dimensionality reduction, and concatenation of each graph's severity label, then builds a K-nearest-neighbour graph, a minimum spanning tree, and a shortest-path trajectory between the lowest- and highest-severity endpoints. A temporal graph convolutional network is trained on this trajectory to predict severity and to produce node-sensitivity scores. On the same trajectory, a neural stochastic differential equation with learned drift $\psi_\theta$ and diffusion $\xi_\phi$ is fitted, and stability is defined by the time-averaged squared change of the diffusion term, while bifurcation conditions are evaluated through a stability potential $J(x,t)=-\int \psi_\theta(x,t)\,dx$ and the diffusion variance. The load-bearing object is the learned diffusion term, because it carries both the stability ranking and the step-4 bifurcation claim.

What would settle it

Randomly permute the 23 severity labels and rerun GPA and the neural SDE; if a step-4 bifurcation appears under permuted labels as often as under real labels, the point of no return is an artifact of the forced low-to-high severity ordering, whereas its disappearance under permutation would support the claim.

Watch

Extended reading notes

Core claim

The central claim is that disease progression can be reconstructed at the pathway level from cross-sectional graphs and then analysed for dynamic structure. On the estimated pseudotime trajectory, the paper defines pathway stability as the time-averaged squared change in the diffusion term of a learned Itô SDE, and defines a bifurcation point as the first time one of three conditions holds: a vanishing stability potential with deteriorating trend, an irreversible transition to a new steady state, or exploding diffusion variance. Applying this to the JR5558 mouse dataset, all major pathways satisfy the transition-to-a-new-steady-state condition at step 4, which the paper interprets as a system-wide irreversible transcriptomic shift and therefore the latest effective time for intervention. The paper also reports stable pathways such as FoxO signaling and non-stable pathways such as caffeine metabolism and phototransduction, and it extends the stability and bifurcation definitions to interacting pathways.

Load-bearing premise

The load-bearing premise is that the ordering of mice produced by the embedding-and-trajectory algorithm is real disease time and not just a re-encoding of the severity score, since every stability score and the step-4 bifurcation is indexed by that ordering.

Editorial extensions

If this is right

  • If the step-4 bifurcation is real, intervention in JR5558 retinal degeneration must occur before that transition; after it, the transcriptomic shift is treated as irreversible by the model.
  • The stable versus non-stable pathway ranking gives a short list of candidate intervention targets, with non-stable pathways such as caffeine metabolism and phototransduction being the most responsive to perturbation.
  • The formal conditions in the paper's definition of a point of no return turn that concept from a metaphor into a computable predicate on an estimated trajectory.
  • The interaction extension implies that positive pathway-to-pathway coupling tends to stabilise the system by homogenising feature variation, while negative couplings can drive divergence.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The authors do not validate the pseudotime axis against an independent temporal marker; one direct test would be to compare the inferred ordering with serial fundus images or longitudinal transcriptomes in a second cohort, which would turn the step-4 threshold into a clinically actionable timing rule.
  • The paper explicitly flags that the pathway-annotation database is too broad for retinal tissue; re-running the pipeline under a retina-specific ontology would show how much of the stability ranking depends on the node definitions.
  • Because the neural SDE is trained on only 23 ordered graphs, the diffusion variance used for stability is sensitive to the ordering; bootstrapping permutations of the trajectory would give confidence intervals for the step-4 bifurcation.
  • The same formal definitions could be exported to other progressive diseases with cross-sectional molecular data, turning 'point of no return' into a general quantity to estimate for early intervention windows.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper introduces a new JR5558 mouse retinal-degeneration transcriptomics dataset (23 bulk RNA-seq samples with lesion-severity scores) and proposes a pipeline that builds pathway graphs from KEGG annotations, identifies disease-sensitive pathways via graph regression, orders the graphs with a new Graph-level Pseudotime Analysis (GPA), fits temporal GCNs to the resulting sequence, and models pathway-feature dynamics with neural stochastic differential equations. From the fitted SDEs it derives pathway stability scores and defines a disease bifurcation point via Definition 1, reporting a system-wide point of no return at Step 4 of the estimated trajectory and ranking stable versus non-stable pathways. The paper also contains an extension to graph-interacted pathway dynamics in Appendix C.2.

Significance. If the framework were valid, it would provide a genuinely useful way to turn static bulk transcriptomic profiles into pathway-level dynamical statements, and the Step-4 bifurcation claim would be a concrete, testable prediction about the therapeutic window in JR5558 retinal degeneration. The curated JR5558 dataset itself is a potentially valuable resource. However, the central results are not currently supported: the GPA time axis is constructed with the severity label inside the embedding, the bifurcation criterion in Definition 1 reduces to an arbitrary displacement threshold, and the stability lemma in Appendix B is not proved under verifiable conditions. These are load-bearing issues for the headline empirical findings rather than presentation details.

major comments (4)
  1. [Section 3.3 (GPA construction)] The GPA time axis is constructed with the target label inside it: z(i) = [\hat z(i) || y_i], and the MST trajectory is anchored at the objects with lowest and highest disease severity. Consequently, the ordering that defines the pseudotime is at least partly a re-parameterization of the severity score y_i, and the temporal GCN of Section 3.4 then predicts y(t) from a sequence sorted by y. This circularity is load-bearing: the neural SDE trajectories of Section 4, the stability scores of Eq. (8), and the Step-4 bifurcation reported in Section 5.1 are all indexed by this ordering, so they do not provide independent evidence about disease dynamics. Section 6 does not acknowledge this issue.
  2. [Appendix B, Lemma 2] The proof of Lemma 2 derives only the one-step Itô increment identity Eq. (14) and then invokes Khasminskii Theorem 6.13 to conclude mean-square divergence or stability. The theorem's hypotheses, such as Lipschitz and linear-growth conditions on the coefficients and the existence of a suitable Lyapunov function, are not verified for the neural-network drift ψ_θ and diffusion ξ_ϕ. The ratio bound in Eq. (10) is an assumption, not a property shown to hold for the learned networks. As written, the lemma is not proved, and the pathway stability metric PS(x) in Eq. (8) lacks the claimed theoretical grounding.
  3. [Section 4, Definition 1] Condition (2) of Definition 1, lim_{t→∞} E||x(t) − x(t*)|| > C with C = 1 set by hand, is not a marker of irreversibility. For any trajectory that drifts monotonically toward a different state, this limit is positive for essentially every t*, not only at a bifurcation. The paper gives no procedure for estimating the t→∞ limit from the 23-step fitted SDE, so in practice the condition reduces to checking whether cumulative displacement exceeds C = 1 in raw feature units. This also explains the Section 5.1 report that all major pathways exhibited their own bifurcation at Step 4: all pathways move along the same pseudotime and cross the same absolute threshold at nearly the same index. The Step-4 bifurcation claim is therefore not supported by the stated criterion.
  4. [Section 4, Eq. (7)] The Itô SDE model is introduced as an inductive assumption, and no model-checking is reported: there is no goodness-of-fit test, no comparison against a non-SDE baseline for the 23-point trajectories, and no uncertainty quantification on the fitted drift and diffusion. Since the trajectories contain only 23 steps total, the neural SDE parameters are severely underdetermined, and the downstream stability rankings and bifurcation times inherit this unvalidated model choice. This is load-bearing for the paper's formal claims.
minor comments (5)
  1. [Section 6] The limitations section acknowledges dataset and annotation limitations but does not mention the severity-label leakage in the GPA construction or the sensitivity of the bifurcation step to the arbitrary constants C and δ; these should be addressed or explicitly discussed.
  2. [Section 3.3 vs. Appendix C.1] Section 3.3 says the trajectory is obtained by fixing pairs of starting and ending objects, but Appendix C.1 says two starting points with zero severity and ending points were fixed and reports four trajectories; the number of endpoints is inconsistent.
  3. [Section 3.2, Eq. (2)] The sensitivity score is a sum of Jacobian norms with no normalization over feature dimensions or graphs; the definition should be made precise and its scale dependence should be discussed.
  4. [Appendix C.2, Corollary 1] The statement that nonnegative A_ij will keep homogenizing every node feature so that x_j = x_i as t→∞ is not a general consequence of the graph-convolved SDE in Eq. (23) unless additional assumptions on the graph Laplacian are imposed, and the proof is only a sketch.
  5. [Appendix B] There is a typo in the last sentence of the proof: 'this completes the graph' should read 'this completes the proof'; also, 'Theorem 6.13' from Khasminskii [2012] is cited without a theorem statement or page reference.

Circularity Check

2 steps flagged · score 6.0 of 10

The Step-4 bifurcation and the GPA disease-time axis reduce to construction choices: a hand-set displacement threshold and a severity-label-informed embedding.

  1. fitted input called prediction [Section 3.3, GPA trajectory construction (Eq. 4 and MST paragraph)]
    "to sufficiently leverage the disease severity information [Shi et al., 2024], we further concatenate each bz(i) with its yi, i.e., z(i) = [ bz(i)||yi]. ... by fixing pairs of starting and ending objects, e.g., objects with the lowest/highest disease severity, the disease evolution trajectory can be obtained by leveraging the shortest path distance algorithm"

    The GPA ordering is used as the time index t for the TGCN forecast (Eq. 5), for identifying stage-transition sensitive pathways, and for fitting the neural SDE that underlies stability and bifurcation. This ordering is constructed by concatenating the disease-severity label y_i into the UMAP embedding and by fixing the MST endpoints as the least and most severe objects. The time axis is therefore a function of the target variable y that the models are later said to predict. Any trajectory, stage transition, or dynamical quantity computed along this axis is partly determined in advance by the severity labels, so the inferred disease dynamics are not independent evidence of progression.

  2. self definitional [Section 4, Definition 1 condition (2); Section 5.1]
    "Specifically, it is the first time at which at least one of the following conditions holds: ... (2): Transition to a New Steady State: limt→∞ E∥x(t) − x(t∗)∥ > C, C >0; ... Although the quantity of C and δ vary between datasets, in our implementation, we simply set C = 1 and δ = 0 .1 ... a key outcome of our study was the identification of a bifurcation point at Step 4 ... all major pathways exhibited their own bifurcation at this step by satisfying our proposed condition (2)"

    Condition (2) is satisfied whenever the trajectory eventually ends more than C away from x(t*). Any non-constant trajectory that converges to a different state passes this check for all sufficiently early t*, whether or not a qualitative transition occurs there. The paper gives no procedure for evaluating the t→∞ limit from the fitted 23-step SDE, and C = 1 is simply imposed. The reported 'bifurcation at Step 4' is therefore a distance-threshold property of the fitted curve, not an independent marker of irreversibility. Since every pathway is driven along the same monotone GPA time axis, all pathways pass the same threshold near the same index, which is exactly the reported 'all major pathways' phenomenon.

full rationale

The paper is not a self-citation chain: references such as [Shi et al., 2024] support a design choice rather than a uniqueness theorem, and no external benchmark is claimed. The circularity is internal to the construction. The GPA pseudotime is built from the severity label (y_i is concatenated into the embedding and endpoints are chosen by severity), yet it is then treated as ground-truth disease time for predicting severity and for defining all dynamics, making the trajectory and stage-transition findings partially self-definitional. In addition, the headline bifurcation at Step 4 is forced by Definition 1, condition (2): with C = 1 chosen by hand and no evaluation of the t→∞ limit, the condition reduces to a distance check that any monotone trajectory passes at early times, so 'all major pathways bifurcate at Step 4' is a threshold artifact rather than a detected critical transition. The dataset curation and the graph-regression sensitivity analysis retain independent content, so the circularity is partial rather than total; a score of 6 reflects that the paper's two most prominent dynamical predictions reduce to construction choices.

Assumptions & free parameters 6 free parameters · 6 assumptions · 0 invented entities

The central claims rest on hand-set thresholds (C=1, delta=0.1 in Definition 1), the Gaussian kernel variance (Appendix C.1), the KNN parameter K=10 (Appendix C.1), the GMM stage count (Figure 3a), and the neural SDE weights fitted to the 23-step pseudotime sequence (Section 4). The axioms include the validity of the label-influenced pseudotime, the Ito-SDE model class for pathway features, the unverified ratio bound (eq. 10), the applicability of Khasminskii's theorem to learned neural networks, and the meaningfulness of Euclidean distances in 9703-dimensional FPKM space with KEGG overlap edges. No new physical or biological entities are postulated: PS, GPS, and the point of no return are composite metrics defined from fitted quantities, so invented_entities is empty.

free parameters (6)
  • Bifurcation threshold C = 1
    Definition 1 condition (2) sets C=1 by hand; combined with the fitted SDE trajectory it determines the reported step-4 bifurcation.
  • Diffusion variance threshold delta = 0.1
    Definition 1 condition (3) sets delta=0.1 by hand; the bifurcation time is sensitive to this choice.
  • Gaussian kernel variance for edge weights = 1
    Appendix C.1 defines W_ij = exp(-||x_i - x_j||^2), 'intrinsically let the variance of the Gaussian kernel as 1'; this determines all graph edge weights used downstream.
  • KNN neighborhood size K in GPA = 10
    Appendix C.1 sets K=10 for the KNN graph that feeds the MST; this shapes the estimated trajectory.
  • Number of GMM stages = 4 (from Figure 3a)
    The GMM clustering segments disease into stages and defines stage-transition pairs; the text does not state how the component count was chosen.
  • Neural SDE drift and diffusion weights = Trained on the 23-step pseudotime sequence
    Section 4: stability rankings and bifurcation conditions are functions of the fitted diffusion network xi_phi.
assumptions (6)
  • domain assumption The GPA-estimated pseudotime orders the graphs along a true disease progression axis.
    Section 3.3 treats the MST/shortest-path ordering as time t for the TGCN and the neural SDE; if the ordering is an artifact of the label-containing embedding, the temporal analysis is circular.
  • ad hoc to paper Pathway feature evolution follows an Ito stochastic differential equation with learnable drift and diffusion.
    Section 4: 'Based on our inductive assumption that the transaction between graph pathway features is governed by stochastic differential equations'; this model class is asserted, not tested against alternatives.
  • ad hoc to paper The drift-to-diffusion ratio bound (eq. 10) holds for the learned drift and diffusion networks.
    Appendix B Lemma 2 assumes C1 <= ||psi||^2 / Tr(xi xi^T) <= C2 for all x and t, but the proof never uses or verifies it; it exists to enable the Khasminskii citation.
  • standard math Khasminskii's stochastic stability theorem (Theorem 6.13) applies to the learned neural SDE.
    Appendix B imports mean-square stability conclusions from Khasminskii [2012]; the theorem's regularity conditions are not checked for the neural network coefficients.
  • domain assumption Euclidean distance between 9703-dimensional FPKM feature vectors is a meaningful measure of pathway similarity.
    Section 3.1 and Appendix C.1 define edge weights from feature distances in a very high-dimensional space, where distances concentrate; the resulting graph structure is not validated.
  • domain assumption Bulk RNA-seq FPKM values from 23 mice support population-level trajectory and SDE inference.
    Section 2: n=23, only FPKM normalization, no batch correction, no differential expression or QC; all downstream dynamics are fit on these values.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Graph Pseudotime Analysis and Neural Stochastic Differential Equations for Analyzing Retinal Degeneration Dynamics and Beyond." pith.science (2026). https://pith.science/paper/RG7TJU7G

@misc{pith2026250206126,
  author       = {Pith},
  title        = {Pith review of: Graph Pseudotime Analysis and Neural Stochastic Differential Equations for Analyzing Retinal Degeneration Dynamics and Beyond},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RG7TJU7G}},
  note         = {Machine review of arXiv:2502.06126}
}
read the original abstract

Understanding disease progression at the molecular pathway level usually requires capturing both structural dependencies between pathways and the temporal dynamics of disease evolution. In this work, we solve the former challenge by developing a biologically informed graph-forming method to efficiently construct pathway graphs for subjects from our newly curated JR5558 mouse transcriptomics dataset. We then develop Graph-level Pseudotime Analysis (GPA) to infer graph-level trajectories that reveal how disease progresses at the population level, rather than in individual subjects. Based on the trajectories estimated by GPA, we identify the most sensitive pathways that drive disease stage transitions. In addition, we measure changes in pathway features using neural stochastic differential equations (SDEs), which enables us to formally define and compute pathway stability and disease bifurcation points (points of no return), two fundamental problems in disease progression research. We further extend our theory to the case when pathways can interact with each other, enabling a more comprehensive and multi-faceted characterization of disease phenotypes. The comprehensive experimental results demonstrate the effectiveness of our framework in reconstructing the dynamics of the pathway, identifying critical transitions, and providing novel insights into the mechanistic understanding of disease evolution.

Figures

Figures reproduced from arXiv: 2502.06126 by the authors.

Figure 1
Figure 1. Illustration of the key steps of the quantification method. (A) A circular region of interest [PITH_FULL_IMAGE:figures/full_fig_p002_1.png] view at source ↗
Figure 2
Figure 2. Illustration on how the pathway graphs are formed and how SPs are identified through two [PITH_FULL_IMAGE:figures/full_fig_p003_2.png] view at source ↗
Figure 3
Figure 3. Process of generating graph trajectory via GPA. (a): Graph dimensionality reduction and [PITH_FULL_IMAGE:figures/full_fig_p005_3.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

13 extracted references · 5 canonical work pages

  1. [1]

    Additionally, prior works in neural stochastic differential equations [Kidger et al., 2021] enforce similar regularization strategies to maintain numerical sta- bility

    The additional assumption is com- monly observed in physical and biological systems, such as Langevin dynamics and gene regulatory networks [Elowitz and Leibler, 2000]. Additionally, prior works in neural stochastic differential equations [Kidger et al., 2021] enforce similar regularization strategies to maintain numerical sta- bility. C Experiment Detail...

    work page 2000

  2. [2]

    in which graphs are formed by the genetic pathways (as nodes) and pathway similarities (as edges), serving as the profiles of the experimental objects (e.g., mice). In addition, we also label these graphs with mice’s lesion severity scores (as a continuous variable) obtained from fundus photographs of these mice, where severity was quantified by measuring...

    work page 2014

  3. [6]

    Inductive represen- tation learning on temporal graphs

    Da Xu, Chuanwei Ruan, Evren Korpeoglu, Sushant Kumar, and Kannan Achan. Inductive represen- tation learning on temporal graphs. arXiv preprint arXiv:2002.07962,

    arXiv 2002

  4. [7]

    Equivariant graph neural operator for modeling 3d dynamics

    Minkai Xu, Jiaqi Han, Aaron Lou, Jean Kossaifi, Arvind Ramanathan, Kamyar Azizzadenesheli, Jure Leskovec, Stefano Ermon, and Anima Anandkumar. Equivariant graph neural operator for modeling 3d dynamics. arXiv preprint arXiv:2401.11037,

    arXiv

  5. [13]

    For the hyperparameters, we let the hidden dimension of the GCN model as 64 and the dropout ratio as 0.5 with the learning rate as 1e−3 and weight decay as 1e−4

    = GCN(X(ℓ), bA), by = MLP(MeanPool(X(ℓ + 1)), (16) and we let ℓ be the number of layers, which is set as 2 in our modeling. For the hyperparameters, we let the hidden dimension of the GCN model as 64 and the dropout ratio as 0.5 with the learning rate as 1e−3 and weight decay as 1e−4. The model is trained with 10 runs in which every run owns 200 epochs. A...

    work page 2019

  6. [2012]

    We start by discretizing the SDE with a small step ∆t, we obtain: x(t + ∆t) = x(t) + ψθ(x(t), t)∆t + ξϕ(x(t), t)∆B(t)

    in which the stability of SDE is analyzed in a more solid manner. We start by discretizing the SDE with a small step ∆t, we obtain: x(t + ∆t) = x(t) + ψθ(x(t), t)∆t + ξϕ(x(t), t)∆B(t). (11) Taking the expectation of the squared difference, we get: E∥x(t + ∆t) − x(t)∥2 = E∥ψθ(x(t), t)∆t + ξϕ(xt, t)∆B(t)∥2. (12) By independence of the drift and diffusion te...

    work page 2012

  7. [2013]

    Umap: Uniform manifold approximation and projection for dimension reduction

    Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426,

    arXiv

  8. [2016]

    Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting

    Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. arXiv preprint arXiv:1709.04875,

    arXiv

Show all 13 references
  1. [2019]

    By considering the so-called adjacency information stored in the graph, GNNs propagate graph node features by aggregating its neighboring information [Wu et al., 2020]

    A Related Works Graph Neural Networks and Different Types of Graphs GNNs were originally proposed to resolve the challenge of data point dependencies via the traditional convolution neural networks, which, in general, treat every input data point independently of each other [K...

    2017

  2. [2020]

    Diffusion convolutional recurrent neural net- work: Data-driven traffic forecasting

    9 Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural net- work: Data-driven traffic forecasting. arXiv preprint arXiv:1707.01926,

    arXiv

  3. [2021]

    Variational graph auto-encoders

    Thomas N Kipf and Max Welling. Variational graph auto-encoders. arXiv preprint arXiv:1611.07308,

    arXiv

  4. [2023]

    A graph autoencoder approach to causal structure learning

    Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang. A graph autoencoder approach to causal structure learning. arXiv preprint arXiv:1911.07420,

    1911 arXiv

  5. [2024]

    Exposition on over-squashing problem on GNNs: Current methods, benchmarks and challenges

    Dai Shi, Andi Han, Lequan Lin, Yi Guo, and Junbin Gao. Exposition on over-squashing problem on GNNs: Current methods, benchmarks and challenges. arXiv:2311.07073,

    arXiv

Pith tools

Reviewed August 8, 2026 · model on record in the stance chip above.