Reconstructing disease dynamics for mechanistic insights and clinical benefit

Abstract

Diseases change over time, both phenotypically and in their underlying molecular processes. Though understanding disease progression dynamics is critical for diagnostics and treatment, capturing these dynamics is difficult due to their complexity and the high heterogeneity in disease development between individuals. We present TimeAx, an algorithm which builds a comparative framework for capturing disease dynamics using high-dimensional, short time-series data. We demonstrate the utility of TimeAx by studying disease progression dynamics for multiple diseases and data types. Notably, for urothelial bladder cancer tumorigenesis, we identify a stromal pro-invasion point on the disease progression axis, characterized by massive immune cell infiltration to the tumor microenvironment and increased mortality. Moreover, the continuous TimeAx model differentiates between early and late tumors within the same tumor subtype, uncovering molecular transitions and potential targetable pathways. Overall, we present a powerful approach for studying disease progression dynamics-providing improved molecular interpretability and clinical benefits for patient stratification and outcome prediction.

Fig. 1. TimeAx discovers the shared disease dynamics across multiple patients.

A Disease progression is not captured by patient subtyping. While current patient stratification and disease state comparison requires clustering of patients into coarse subgroups (top), understanding disease dynamics as a common ground between patients’ disease trajectories enables patient stratification in a higher resolution (bottom). B TimeAx seed features selection process. TimeAx selects a “conserved-dynamics-seed” with shared dynamics across all subjects. In the heatmap (left), each row is a patient and each column is a sample collected at a different time. Each feature spans multiple patients (multiple rows) and is also presented as a line plot (right). C An illustration describing TimeAx’s utility for revealing the shared disease progression dynamics and the projection of patients-specific disease pseudotime positions, integrating the differences between patients’ disease trajectories over time. D Disease pseudotime can be utilized to discover novel disease mechanisms as well as to support new clinical frameworks for patient stratification and outcome prediction. We thank Yuval Abraham for his contribution in the design and creation of (A, C and D).

Fig. 2. Disease pseudotime captures disease progression dynamics better than chronological time.

A An illustration of the influenza infection dynamics TimeAx modeling and disease pseudotime inference, based on the longitudinal influenza cohort. B TimaAx disease pseudotime (y-axis) is different from chronological time (time from infection, x-axis). Shown are the differences between symptomatic (green lines) and asymptomatic (red lines) patients from the longitudinal influenza cohort. P value was calculated by comparing the prediction of disease pseudotime (by ANOVA), using only time or the interaction between time and symptoms as predictors. Trend lines represent the average levels in each of the groups ± standard error. C Gene associations (−log10 transformed, FDR-corrected, Q values based on linear regression), across all symptomatic and asymptomatic patients, using either sampling time (x-axis) or disease pseudotime (y-axis), and applying a Q value threshold of 10⁻⁵ (Dashed lines) (See Methods). Genes are colored based on their association with the two time axes. D An illustration of the UBC dynamics TimeAx modeling and disease pseudotime inference, based on the UBC longitudinal cohort. E TimeAx disease pseudotime (y-axis) is different from chronological time (time from primary tumor, x-axis), exemplified in six patients in the UBC longitudinal cohort. F Gene associations (−log10 transformed, P values based on linear regression) with sampling time (x-axis) and disease pseudotime (y-axis), using a P value threshold of 10⁻² (Dashed lines) (See Methods). Genes are colored based on their association with the two time axes, displaying significant associations almost entirely only with the disease pseudotime. G An illustration of the AMD dynamics TimeAx modeling and disease pseudotime inference, based on segmented features extracted from OCT scans of the patients’ retina. Optical coherence tomography (OCT). H TimeAx disease pseudotime (y-axis) is different from sampling time (time from first encounter, x-axis), exemplified in five patients in the AMD train cohort. I The distribution of AMD test cohort patients’ disease pseudotime positions (top; n = 29205 biologically independent samples) and times from diagnosis (bottom; n = 11075 biologically independent samples) (y-axis), across different visual severity states, determined according to the patients’ visual acuity levels (logMAR; x-axis). Boxes represent the 25th, 50th, and 75th percentiles; whiskers show maxima and minima. P values were calculated based on linear regression. We thank Yuval Abraham for his contribution in the design and creation of (A, D and G).

Fig. 3. TimeAx uncovers an advanced tumor state with unfavorable clinical outcomes.

A Disease pseudotime clinical applicability. Shown are disease pseudotime (y-axis) relation with tumor stage (x-axis) for samples within the UBC longitudinal cohort. P value was calculated based on linear regression. B Tumor purity scores (y-axis) along the disease pseudotime (left) and the time from primary tumor (right) (x-axis), displaying a sharp decrease in high disease pseudotime positions (SPIP; dashed line). Disease pseudotime ranges pre and post the SPIP are marked by colored bars. C Cell type deconvolved cell contributions (y-axis), displaying major differences between pre- and post- SPIP samples. * p < 0.05, ** p < 0.01 based on a two-sided t-test (p < 10⁻⁴, 0.003, 0.05, 0.05, 10⁻⁴; Urothelial cells, Macrophages, naive T cells, memory T cells and Fibroblasts, respectively). D Survival plot for UBC patients within the TCGA test cohort, comparing patients’ tumors with disease pseudotime positions lower and higher the SPIP (pre versus post, respectively; color coded). P value was calculated based on a log-rank test. E Survived (black) versus deceased (white) percentages of basal/squamous patients within the TCGA test cohort pre and post SPIP (x-axis). p value was calculated using Fisher’s exact test. In (A and C), boxes represent the 25th, 50th, and 75th percentiles; whiskers show maxima and minima and n = 84 biologically independent samples. Stromal pro-invasion point (SPIP).

Fig. 4. Disease pseudotime captures variation undetectable by current stratification frameworks.

A Disease pseudotime distribution across tumor molecular classifications in patients within the UBC longitudinal cohort based on the ‘LundTax’ molecular subtyping framework. B Distribution of ‘LundTax’ tumor molecular classifications for pre (red) and post (blue) stromal pro-invasion samples within the UBC longitudinal cohort. C Comparison of disease pseudotime (y-axis) between primary and recurrent tumors in pre-SPIP samples within the UBC longitudinal cohort. D Disease pseudotime distribution across tumor molecular classifications in patients pre-SPIP within the UBC longitudinal cohort based on the ‘LundTax’ molecular subtyping framework. E Percent of surviving patients in the TCGA test cohort with UroA tumors (y-axis) across disease pseudotime bins (bin size = 0.1; x-axis) pre-SPIP. p value was calculated based on linear regression. In (A, C and D), boxes represent the 25th, 50th, and 75th percentiles; whiskers show maxima and minima and n = 84 biologically independent samples.

Fig. 5. Disease pseudotime uncovers molecular mechanisms promoting UBC progression.

A Co-expression matrix for genes (rows, columns) differentially expressed (q < 0.05) between early and late UroA tumors. Two gene sets were discovered and are highlighted in green and purple. B High enrichment of different pseudogene families within the downregulated module, compared to the upregulated module set and a random gene set (color-coded). C Pathway enrichment for the downregulated module, including pathway enrichment scores (left) and a heatmap of the expression levels of these pathways genes (columns), in early and late UroA tumors (rows) within the UBC longitudinal cohort. D Pathway enrichment scores for the upregulated module. E An illustration of the molecular model of UBC tumorigenesis, discovered by the TimeAx-based analysis, based on the increasing gene module. The illustration includes the association of known hallmarks of cancer with disease pseudotime (top) and the suggested mitotic kinetochore- spindle microtubules (MT) interaction (middle), which is also presented as boxplots of differential expression of its subunits (bottom; n = 40 biologically independent samples, boxes represent the 25th, 50th, and 75th percentiles; whiskers show maxima and minima). G protein-coupled receptors (GPCRs). We thank Yuval Abraham for his contribution in the design and creation of (E).