WO2023225618A2

WO2023225618A2 - Method for estimating a dynamic molecular program of a cell

Info

Publication number: WO2023225618A2
Application number: PCT/US2023/067200
Authority: WO
Inventors: Manik KUCHROO; Alex TONG; Smita Krishnaswamy; Christine CHAFFER
Original assignee: Garvan Institute of Medical Research; Yale University
Current assignee: Garvan Institute of Medical Research; Yale University
Priority date: 2022-05-18
Filing date: 2023-05-18
Publication date: 2023-11-23
Anticipated expiration: 2024-11-18
Also published as: WO2023225618A3; US20250308625A1

Abstract

Aspects of the present invention relate to a method of estimating a dynamic molecular program of a population of cells including the steps of providing a set of at least two static snapshots of a population of cells undergoing a state transition at a corresponding set of at least two time indices to a neural network, calculating a set of possible population flows between the at least two time indices based on the at least two static snapshots, negatively weighting any of the set of population flows which are unrealistic, and inferring an estimated population flow of the cells between the set of static snapshot data by selecting a population flow from the set of possible population flows with the neural network.

Description

METHOD FOR ESTIMATING A DYNAMIC MOLECULAR PROGRAM OF A CELL

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to U.S. Provisional Application No. 63/343,142, filed on May 18, 2022, incorporated herein by reference in its entirety.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

[0002] This invention was made with government support under 1F30-AI157270 and 1R01GM135929 awarded by the National Institutes of Health and 2047856 awarded by the National Science Foundation. The government has certain rights in the invention.

BACKGROUND OF THE INVENTION

[0003] Cell state plasticity provides a mechanism for cancer cells to rapidly and dynamically evolve in a manner that facilitates primary tumor growth, metastasis and the development of therapy -resistant disease. While the advent of single-cell technologies have allowed detailed characterization of static cell states within tumors, further technological development is required to elucidate the mechanisms governing dynamic cell state transitions that may not only span days, months or years, but also create a variety of cell states shaping the tumor heterogeneity that drives disease progression. Furthermore, the transcriptional networks used by individual cells to undergo dynamic functional changes have been difficult to dissect due to the high dimensional nature of the data, as well as the computational challenge of resolving cellular trajectories over extended periods of time from static snapshot single-cell data. If these issues were addressed, it would be possible to use longitudinal patient samples to gain an unprecedented insight into the mechanisms governing metastasis and therapy-resistant disease, both of which have eluded scientists for decades.

[0004] Thus, there is the need in the art for a method of elucidating the mechanisms governing dynamic cell state transitions while addressing the computational challenge of resolving cellular trajectories over extended periods of time from a static snapshot of single-cell data. The present invention meets this need.

SUMMARY OF THE INVENTION

[0005] Aspects of the present invention relate to a method of estimating a dynamic molecular program of a population of cells including the steps of providing a set of at least two static snapshots of a population of cells undergoing a state transition at a corresponding set of at least two time indices to a neural network, calculating a set of possible population flows between the at least two time indices based on the at least two static snapshots, negatively weighting any of the set of population flows which are unrealistic, and inferring an estimated population flow of the cells between the set of static snapshot data by selecting a population flow from the set of possible population flows with the neural network.

[0006] In some embodiments, the method further includes an Ordinary Differential Equation (ODE) solver to form a neural ODE system. In some embodiments, the state transition is a mesenchymal-to-epithelial transition (MET), or epithelial-to-mesenchymal transition (EMT). In some embodiments, the calculating step further includes the steps of approximating a timevarying derivative parametrized by a set of network weights and biases and integrating the timevarying derivative across the at least two time indices to calculate a possible population flow of the set of possible population flows.

[0007] In some embodiments, the method further includes the step of limiting a magnitude of the time-varying derivative across the at least two time indices to a predetermined threshold. In some embodiments, the set of at least two static snapshots comprise three-dimensional tumorsphere data. In some embodiments, the method further includes the step of identifying shared and trajectory-specific molecular programs using a time-lapsed causality analysis. In some embodiments, the at least two time indices of the at least two static snapshots are separated by at least 24 hours. In some embodiments, the dynamic molecular program is a gene regulatory network. In some embodiments, the method further includes the step of calculating an interpolated snapshot of the cell population at a third time index based on the inferred population flow. In some embodiments, the set of population flows which are unrealistic comprise energy inefficient biological pathways. [0008] Aspects of the present invention relate to a method of preventing a mesenchymal cell from differentiating into an epithelial cell having the step of contacting the cell with one or more modulator of one or more molecule involved in the mesenchymal-to-epithelial transition (MET) transcriptional network. In some embodiments, the one or more molecules in the MET transcriptional network are one or more selected from the group consisting of: estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAI1 (SNAI1), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1).

[0009] Aspects of the present invention relate to a method of directing one or more cells in a population of cells to a transition state, having the steps of obtaining a population of cells from the subject, identifying a target state transition for the population of cells to undergo, administering at least one modulator of at least one molecule within a MET transcriptional network of the population of cells thereby directing the population of cells to the targeted state transition.

[0010] In some embodiments, the at least one molecule in the MET transcriptional network is selected from the group consisting of: estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAI1 (SNAI1), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1).

[0011] In some embodiments, the targeted state transition is energy optimal.

BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The foregoing purposes and features, as well as other purposes and features, will become apparent with reference to the description and accompanying figures below, which are included to provide an understanding of the invention and constitute a part of the specification, in which like numerals represent like elements, and in which: Figs. 1 A through 1C shows an overview of an algorithm for estimating a dynamic molecular program of a cell according to aspects of the present invention (i.e., TrajectoryNet Algorithm). Fig. 1A shows data for TrajectoryNet is generated by capturing time lapsed single-cell data (RNA-seq, ATAC-seq, CITE-seq, etc.) on an evolving system. Fig. IB is an illustration of the TrajectoryNet model. TrajectoryNet uses an adaptive ODE solver to integrate cell state over time in order to calculate trajectories and relative proliferation rates. Time series single cell data produces disconnected distributions over a developmental time. TrajectoryNet interpolates disconnected distributions to continuously infer transcriptional dynamics as well as cell growth and death via unbalanced dynamic optimal transport. Fig. 1C shows an overview of TrajectoryNet analysis framework: i) identify continuous trajectories from disconnected time- lapsed data; ii) interpolate continuous gene expression dynamics based on terminal cell population; iii) Build transcriptional networks using causality analysis and public gene regulatory network databases.

Figs. 2A through 2G shows an overview of Tumorsphere dataset. Fig. 2A is a schematic illustrating dynamic transitions via the EMT and MET between highly tumorigenic and metastatic CD44^hiZEB 1 ^hiCDH 1 ^lo CSCs and poorly tumorigenic CD44^loZEB 1 ^loCDH 1 ^hi epithelial cells. Fig. 2B is an illustration of the tumorsphere protocol and single cell RNAseq experiment. CD44¹¹' CSCs are seeded in single cell suspension at day 1. By day 30, ten percent of single CD44¹¹¹ cells seeded produce three-dimensional heterogeneous tumorspheres. Fig. 2C shows PHATE [Moon, K. R. et al. Visualizing structure and transitions in high-dimensional biological data. Nature Biotechnology 37, 1482-1492 (2019)] embedding of time-lapsed scRNAseq data generated by the tumorsphere assay described in Figure 2B. Samples are coloured by timepoint of data acquisition. Trend lines are created by TrajectoryNet. Fig. 2D shows a visualization of TrajectoryNet inferred proliferation rate. Fig. 2E shows a visualization of EMT gene expression score [Chakraborty, P., George, J. T., Tripathi, S., Levine, H. & Jolly, M. K. Comparative study of transcriptomics-based scoring metrics for the epithelial hybrid-mesenchymal spectrum. Frontiers in Bioengineering and Biotechnology 8 (2020)] on PHATE embedding. Fig. 2F shows a visualization of inferred continuous gene expression dynamics over time. Genes have been clustered into 5 groups based on their gene expression dynamics . Key EMT genes (EPCAM, TWIST1/2, SNAI1/2 and ZEB 1/2) have been indicated on heatmap. Fig. 2G shows a visualization of inferred continuous gene expression dynamics in each gene cluster shown in Figure 2C. Each line represents a single gene trend overtime.

Figs. 3 A through 3E shows refining identification of the tumorsphere cell-of-origin. Fig. 3A shows PHATE visualization of day 0 scRNAseq timepoint with CD44 expression highlighted. Fig. 3B shows a zoom in on day 0 CD44¹¹' population with TrajectoryNet proliferation rate (top) and inferred cell cycle state (bottom) [Tirosh, I. et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189-196 (2016)] visualized. Fig. 3C shows flow cytometry-based sorting of HCC38 CD44¹¹¹ cells infected with the FUCCI cell cycle sensor system to isolate cells at different stages of the cell cycle. Each cell cycle isolate is measured for in vitro tumorsphere-initiating potential. Fig. 3D shows a visualization in CD44hi population of key differentially expressed cell surface markers and their DREMI [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in single-cell data. Science 346, 1250689-1250689 (2014)] association scores with TrajectoryNet-inf erred proliferation rate. Fig. 3E shows flow cytometry isolated EPCAM^+/- and CAV1^+/- populations are measured for tumorsphere-initiating potential.

Figs. 4A through 4H shows comparing gene regulation in the EMT and MET trajectories. Fig. 4A shows PHATE visualization of day 30 scRNAseq timepoint with three populations highlighted corresponding to epithelial (orange), mesenchymal (green) and apoptotic (blue) populations (left). Populations were computed with Louvain clustering [Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment 2008, Pl 0008 (2008)] and were identified using mitochondrial (MT) and EMT gene signatures (right). Fig. 4B shows microscopy visualization of CDH1 and ZEB1 at day 7 and day 28 of the tumorsphere assay. Dashed white lines highlight distinct populations with high ZEB1 (mesenchymal) or CDH1 (epithelial) expression. Fig. 4C shows tracing individual trajectories of cells that undergo EMT (green) and MET (orange). Fig. 4D shows 3D PHATE visualization of day 12, 18, and 30 showing the divergence of the EMT (green) and MET (orange) trajectories highlighting their mapping to discrete cell populations at day 30. Fig. 4E shows a heatmap of signed Granger values between 5,273 genes and 461 transcription factors organized by gene clusters from Figure 2C for combined trajectories. Red denotes a strong enhancing relationship between a transcription factor and a target gene, while blue denotes a strong repressive relationship. Fig. 4F shows the top regulatory transcription factors from EMT, MET and combined trajectories are separated by gene cluster (as computed in Figure 2C). The overlap between regulatory transcription factors from each of the EMT, MET and combined trajectories. Fig. 4G shows visualizing inferred continuous gene expression dynamics of representative TFs regulating the EMT and MET. Fig. 4H shows visualisation of gene regulatory networks of core EMT transcription factors (green) and MET transcription factors (orange) highlighting cross-talk through common gene interactors (grey nodes). The core MET transcription factors share regulatory relationships with a hub of MET-specific interactors (yellow nodes). Similarly, the EMT core transcription factors share regulatory relationships with a hub of EMT-specific interactors (purple).

Figs. 5A through 5G shows MET temporal gene network identification and validation. Fig. 5A shows an exemplary schematic of the workflow and filtering strategy used to curate the temporal MET gene regulatory network using TRRUST v2 database (Transcriptional Regulatory Relationships Unravelled by Sentence-based Text mining [Han, H. et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions.

Nucleic Acids Research 46, D380-D386 (2017)]), overlay ed with genes associated with EMT (Epithelial Mesenchymal Transition Gene Database, dbEMT 2.0 [Zhao, M., Liu, Y., Zheng, C. & Qu, H. dbEMT 2.0: An updated database for epithelial-mesenchymal transition genes with experimentally verified information and precalculated regulation information for cancer metastasis. Journal of Genetics and Genomics 46, 595-597 (2019)]). Fig. 5B shows the resultant MET network comprises transcription factors identified by TrajectoryNet (rectangles) across Gene Clusters 1-5. Known E/M plasticity genes are highlighted by ellipses, and novel E/M plasticity genes are marked by diamonds. Fig. 5C shows a visualization of ESRRA gene regulatory module within panel B limited to known gene interactions with ESRRA and the known EMT genes (CDH1, ZEB1, SNAI1/2). Genes upregulated in the MET trajectory are highlighted in pink and genes downregulated are highlighted in grey. Fig. 5D shows a visualization of ESRRA, ZEB1 and CDH1 by immunofluorescence staining at 4 timepoints in 3D tumorspheres. Fig. 5E shows a visualization of TrajectoryNet interpolated gene trends as well as ground truth protein trends for ESRRA, ZEB1 and CDH1. Fig. 5F shows the effect of ESRRA knockdown (siRNA) on CDH1 expression in HCC38 CD44*¹¹ cells by western blot. Fig. 5G §hows the effect of ESRRA inhibition using C14 on CDH1 expression in HCC38 CD44¹¹¹ cells by western blot.

Figs. 6A through 6H shows validating cancer cell plasticity trajectories in vivo. Fig. 6A shows an exemplary schematic of in vivo scRNAseq experimental workflow on primary tumors implanted in the mammary fat pad with matched spontaneous lung metastases from the xenograft MDA- MB-231 triple-negative breast cancer model. Fig. 5B shows an exemplary schematic of learned trajectories from scRNAseq data developing 1) within a primary tumor and 2) from a primary tumor to lung metastasis. Fig. 6C shows expression of current and new CSC markers shown on scRNAseq data from primary tumor and matched lung metastasis. Fig. 6D shows expression of current and new CSC markers shown on scRNAseq data on matched spatial transcriptomic profiling of the primary tumor. Fig. 6E shows identification of primary tumor and lung metastasis CSCs (red) based on high expression of CD4-L EPCAM, CAV1, and ZEB1. Fig. 6F shows EMT score of primary tumors and lung metastases [Chakraborty, P., George, J. T., Tripathi, S., Levine, H. & Jolly, M. K. Comparative study of transcriptomics-based scoring metrics for the epithelial-hybrid-mesenchymal spectrum. Frontiers in Bioengineering and Biotechnology 8 (2020)]. Fig. 6G shows clustering of primary tumors and lung metastases into six subpopulations. Gray cluster represents CSCs. Green cluster represents emergent mesenchymal subpopulation. Fig. 6H shows transcriptional dynamics of core EMT transcription factors in primary tumor trajectory (top row) and primary tumor to lung metastasis trajectory (bottom row).

Figs. 7A through 7D shows TrajectoryNet comparisons on synthetic data of different structures. Fig. 7A shows a depiction of ID manifolds with non-branching and branching. Models are trained on data from timepoints to (blue) and ti (orange) with the goal of accurately interpolating ground truth data at t_1/2. Fig. 7B shows comparing TrajectoryNet with other methods at predicting t_1/2. Earth Mover’s Distance between the predicted data at t_1/2 and the ground truth distribution at t_1/2 shown across models, with lower distances indicating a more accurate prediction of t_1/2. Fig. 7C shows visualizing individual paths projected forward from to across methods. Fig. 7D shows comparing interpolated and ground truth t_1/2 for the three other methods (OT, random, RNA velocity) that do not create paths to visualize. Figs. 8A through 8C shows ordering of cells by TrajectoryNet, scVelo, and Diffusion Pseudotime. Fig. 8A shows cells colored by (from left to right) TrajectoryNet, scVelo, and Diffusion inferred time between zero (purple) and one (yellow). scVelo represents RNA- Velocity based methods and Diffusion Pseudotime represents graph-based pseudotime inference methods. Fig. 8B shows plots of inferred cell time vs. ground truth observation time. It can be seen that TrajectoryNet is the only method where the inferred time correlates with the ground truth observation time. Fig. 8C shows a visualization of scVelo inferred velocity streams across timepoints and embeddings. Each embedding is colored by the inferred cell state (S: Green), (Gl: Blue) and (G2M: Orange). It can be seen that the scVelo inferred velocity streams are inconsistent between embeddings.

Figs. 9A through 9D shows dynamic versus static inference of gene regulatory interactions. Fig. 9A shows an exemplary schematic overview of Granger analysis. Signed Granger values are computed using Granger causality inference [Granger, C. W. J. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424 (1969)] between a potentially regulatory transcription factor and a potentially regulated target gene based on r time lag in regulatory effect. The Granger values are signed based on the direction of effect with an enhancing relationship having a + sign and a repressive relationship having a - sign. Fig. 9B shows visualizations of synthetic datasets, in PC dimensions, and gene networks used to evaluate the disclosed method, total Granger causal score (TGCS) in Fig. 9C. Fig. 9 C shows a visual comparison of the performance of TGCS, against three methods of gene network inference, DREMI [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in singlecell data. Science 346, 1250689 (2014)], Pearson correlation and Spearman correlation, on four synthetic datasets with different known ground truth regulatory interaction structures as specified in Bool ODE [Pratapa, A., Jalihal, A. P., Law, J. N., Bharadwaj, A. & Murali, T. M. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature Methods 17, 147-154 (2020)]. Fig. 9D shows a numerical comparison of TGCS against DREMI, Pearson correlation and Spearman correlation at identifying ground truth gene network structure using area under the receiver operator characteristic curve (AUC ROC). Here higher scores indicate that a method is more accurately able to identify ground truth gene network structure. Results were averaged across 100 runs with one standard deviation error bars visualized. Figs. 1 OA and 1 OB shows gene dynamics and gene network calculations for EMT and MET dynamics. Fig. 10A shows recomputed gene expression dynamics based on trajectories that terminate in mesenchymal cellular cluster identified in Figure 4A (left). Signed Granger analysis between 5,273 genes and 461 transcription factors gene trends of cells that terminate in mesenchymal cluster (right). Red denotes a strong enhancing relationship between a transcription factor and a target gene, while blue denotes a strong repressive relationship. Fig. 10B shows recomputed gene expression dynamics based on trajectories that terminate in epithelial cellular cluster identified in Figure 4A (left). Signed Granger analysis between 5,273 genes and 461 transcription factors gene trends of cells that terminate in epithelial cluster (right). Red denotes a strong enhancing relationship between a transcription factor and a target gene, while blue denotes a strong repressive relationship.

Figs. 11A and 1 IB show visualizing expression dynamics of key mesenchymal and epithelial regulatory genes. Fig. 11A shows the top MET-specific regulatory transcription factors separated by gene clusters: i) cluster 2, ii) cluster 3, iii) cluster 4, iv) cluster 5. Fig. 1 IB shows the top EMT-specific regulatory transcription factors separated by gene clusters: i) cluster 1, ii) cluster 2, iii) cluster 3, iv) cluster 4, v) cluster 5.

Figs. 12A through 12J shows visualizing the extended EMT and MET subnetwork along with the key pathways they regulate. Fig. 12A shows MET subnetwork comprising the core MET transcription factors (orange), MET specific interactors (yellow and interactors common with EMT subnetwork (grey). Figs. 12B through 12D show exemplary gene regulatory relationships of core MET transcription factors ESRRA (Fig. 12B), ARNT (Fig. 12C), ZEB1 (Fig. 12D). Fig. 12E shows EMT subnetwork comprising the core EMT transcription factors (green), EMT specific interactors (purple) and interactors common with EMT core transcription factors (grey). Figs. 12F through 12H show exemplary gene regulatory relationships of core EMT transcription factors HES1 (Fig. 12F), SNAI1 (Fig. 12G), F0X03 (Fig. 12H). Fig. 121 shows pathways enriched in the MET subnetwork. Fig. 12J shows pathways enriched in the EMT subnetwork.

Figs. 13 A and 13B show validating predicted CDH1 expression dependence on ESRRA using siRNA knock-down and C14 antagonist. Fig. 13 A shows Western-blot scans of ESRRA with GAPDH as reference using siRNA knockdown of ESRRA. Quantification can be seen in Figure 5F. Fig. 13B shows Western-blot scans of ESRRA with (I-actin as reference using Cl 4 antagonist of ESRRA. Quantification can be seen in Figure 5G.

Fig. 14 shows spatial expression and dynamics of EMT network genes within in vivo datasets. Visualization of EMT-related genes in scRNAseq data of primary tumor and lung metastasis, and their spatial localization in spatial transcriptomic (Visium lOx) data of a primary tumor.

Fig. 15 is an exemplary computing device.

DETAILED DESCRIPTION

[0013] It is to be understood that the figures and descriptions of the present invention have been simplified to illustrate elements that are relevant for a clear understanding of the present invention, while eliminating, for the purpose of clarity, many other elements found in related systems and methods. Those of ordinary skill in the art may recognize that other elements and/or steps are desirable and/or required in implementing the present invention. However, because such elements and steps are well known in the art, and because they do not facilitate a better understanding of the present invention, a discussion of such elements and steps is not provided herein. The disclosure herein is directed to all such variations and modifications to such elements and methods known to those skilled in the art.

[0014] Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. Although any methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present invention, exemplary methods and materials are described.

[0015] As used herein, each of the following terms has the meaning associated with it in this section.

[0016] The articles “a” and “an” are used herein to refer to one or to more than one (i.e., to at least one) of the grammatical object of the article. By way of example, “an element” means one element or more than one element. [0017] “About” as used herein when referring to a measurable value such as an amount, a temporal duration, and the like, is meant to encompass variations of ±20%, ±10%, ±5%, ±1%, and ±0.1% from the specified value, as such variations are appropriate.

[0018] Throughout this disclosure, various aspects of the invention can be presented in a range format. It should be understood that the description in range format is merely for convenience and brevity and should not be construed as an inflexible limitation on the scope of the invention. Accordingly, the description of a range should be considered to have specifically disclosed all the possible subranges as well as individual numerical values within that range. For example, description of a range such as from 1 to 6 should be considered to have specifically disclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numbers within that range, for example, 1, 2, 2.7, 3, 4, 5, 5.3, 6 and any whole and partial increments therebetween. This applies regardless of the breadth of the range.

[0019] In some aspects of the present invention, software executing the instructions provided herein may be stored on a non-transitory computer-readable medium, wherein the software performs some or all of the steps of the present invention when executed on a processor.

[0020] Aspects of the invention relate to algorithms executed in computer software. Though certain embodiments may be described as written in particular programming languages, or executed on particular operating systems or computing platforms, it is understood that the system and method of the present invention is not limited to any particular computing language, platform, or combination thereof. Software executing the algorithms described herein may be written in any programming language known in the art, compiled or interpreted, including but not limited to C, C++, C#, Objective-C, lava, JavaScript, MATLAB, Python, PHP, Perl, Ruby, or Visual Basic. It is further understood that elements of the present invention may be executed on any acceptable computing platform, including but not limited to a server, a cloud instance, a workstation, a thin client, a mobile device, an embedded microcontroller, a television, or any other suitable computing device known in the art.

[0021] Parts of this invention are described as software running on a computing device. Though software described herein may be disclosed as operating on one particular computing device (e g. a dedicated server or a workstation), it is understood in the art that software is intrinsically portable and that most software running on a dedicated server may also be run, for the purposes of the present invention, on any of a wide range of devices including desktop or mobile devices, laptops, tablets, smartphones, watches, wearable electronics or other wireless digital/cellular phones, televisions, cloud instances, embedded microcontrollers, thin client devices, or any other suitable computing device known in the art.

[0022] Similarly, parts of this invention are described as communicating over a variety of wireless or wired computer networks. For the purposes of this invention, the words “network”, “networked”, and “networking” are understood to encompass wired Ethernet, fiber optic connections, wireless connections including any of the various 802.11 standards, cellular WAN infrastructures such as 3G, 4G/LTE, or 5G networks, Bluetooth®, Bluetooth® Low Energy (BLE) or Zigbee® communication links, or any other method by which one electronic device is capable of communicating with another. In some embodiments, elements of the networked portion of the invention may be implemented over a Virtual Private Network (VPN).

[0023] Fig. 15 and the following discussion are intended to provide a brief, general description of a suitable computing environment in which the invention may be implemented. While the invention is described above in the general context of program modules that execute in conjunction with an application program that runs on an operating system on a computer, those skilled in the art will recognize that the invention may also be implemented in combination with other program modules.

[0024] Generally, program modules include routines, programs, components, data structures, and other types of structures that perform particular tasks or implement particular abstract data types. Moreover, those skilled in the art will appreciate that the invention may be practiced with other computer system configurations, including hand-held devices, multiprocessor systems, microprocessor-based or programmable consumer electronics, minicomputers, mainframe computers, and the like. The invention may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote memory storage devices. [0025] Fig. 15 depicts an illustrative computer architecture for a computer 1000 for practicing the various embodiments of the invention. The computer architecture shown in Fig. 15 illustrates a conventional personal computer, including a central processing unit 1050 (“CPU”), a system memory 1005, including a random access memory 1010 (“RAM”) and a read-only memory (“ROM”) 1015, and a system bus 1035 that couples the system memory 1005 to the CPU 1050. A basic input/output system containing the basic routines that help to transfer information between elements within the computer, such as during startup, is stored in the ROM 1015. The computer 1000 further includes a storage device 1020 for storing an operating system 1025, application/program 1030, and data.

[0026] The storage device 1020 is connected to the CPU 1050 through a storage controller (not shown) connected to the bus 1035. The storage device 1020 and its associated computer-readable media provide non-volatile storage for the computer 1000. Although the description of computer-readable media contained herein refers to a storage device, such as a hard disk or CD- ROM drive, it should be appreciated by those skilled in the art that computer-readable media can be any available media that can be accessed by the computer 1000.

[0027] By way of example, and not to be limiting, computer-readable media may comprise computer storage media. Computer storage media includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storage of information such as computer-readable instructions, data structures, program modules or other data.

Computer storage media includes, but is not limited to, RAM, ROM, EPROM, EEPROM, flash memory or other solid state memory technology, CD-ROM, DVD, or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by the computer.

[0028] According to various embodiments of the invention, the computer 1000 may operate in a networked environment using logical connections to remote computers through a network 1040, such as TCP/IP network such as the Internet or an intranet. The computer 1000 may connect to the network 1040 through a network interface unit 1045 connected to the bus 1035. It should be appreciated that the network interface unit 1045 may also be utilized to connect to other types of networks and remote computer systems.

[0029] The computer 1000 may also include an input/output controller 1055 for receiving and processing input from a number of input/output devices 1060, including a keyboard, a mouse, a touchscreen, a camera, a microphone, a controller, a joystick, or other type of input device. Similarly, the input/output controller 1055 may provide output to a display screen, a printer, a speaker, or other type of output device. The computer 1000 can connect to the input/output device 1060 via a wired connection including, but not limited to, fiber optic, Ethernet, or copper wire or wireless means including, but not limited to, Wi-Fi, Bluetooth, Near-Field Communication (NFC), infrared, or other suitable wired or wireless connections.

[0030] As mentioned briefly above, a number of program modules and data files may be stored in the storage device 1020 and/or RAM 1010 of the computer 1000, including an operating system 1025 suitable for controlling the operation of a networked computer. The storage device 1020 and RAM 1010 may also store one or more applications/programs 1030. In particular, the storage device 1020 and RAM 1010 may store an application/program 1030 for providing a variety of functionalities to a user. For instance, the application/program 1030 may comprise many types of programs such as a word processing application, a spreadsheet application, a desktop publishing application, a database application, a gaming application, internet browsing application, electronic mail application, messaging application, and the like. According to an embodiment of the present invention, the application/program 1030 comprises a multiple functionality software application for providing word processing functionality, slide presentation functionality, spreadsheet functionality, database functionality and the like.

[0031] The computer 1000 in some embodiments can include a variety of sensors 1065 for monitoring the environment surrounding and the environment internal to the computer 1000. These sensors 1065 can include a Global Positioning System (GPS) sensor, a photosensitive sensor, a gyroscope, a magnetometer, thermometer, a proximity sensor, an accelerometer, a microphone, biometric sensor, barometer, humidity sensor, radiation sensor, or any other suitable sensor. [0032] Aspects of the present invention relate to methods and systems for tracing single cells and/or populations of cells to predict or infer future and past cellular states. In some embodiments, the present invention provides a method for predicting gene expression values continuously. In some examples, the present invention is referred to herein as “TrajectoryNet.” In some embodiments, the present invention provides a method (i.e. a pipeline based on concepts from the field of information theory) that leverages these continuous gene expression features to build genomic regulatory networks. In some embodiments, these genomic networks may be used to identify therapeutic targets.

[0033] The disclosed method (i.e. TrajectoryNet) can continuously identify gene expression trends from disconnected single cell data. The disclosed method can then build genomics regulatory networks from this information.

[0034] In some embodiments, the present invention describes one or more cells, or a population of cells that comprises a dynamic molecular program. As used herein, the dynamic molecular program of one or more cells includes, but is not limited to, a genetic pathway that a given cell or cells may follow either in isolation or within a population of cells. In some embodiments, the dynamic molecular program is a gene regulatory network. The modulation of the molecular program may, in some examples, cause the expression or inhibition of genetic information within a cell or population of cells over a period of time. Tn some embodiments, the modulation of the molecular program may drive morphological and/or physiological changes of a cell or population of cells over time.

[0035] In one aspect, the present invention includes a method of estimating a dynamic molecular program of a population of cells comprising the steps of providing a set of at least two static snapshots of a population of cells undergoing a state transition at a corresponding set of at least two time indices to a neural network, calculating a set of possible population flows between the at least two time indices based on the at least two static snapshots, negatively weighting any of the set of population flows which are unrealistic, and inferring and/or calculating an estimated population flow of the population of cells between the set of static snapshot data by selecting a population flow from the set of possible population flows with the neural network. [0036] Tn some embodiments, the neural network is an ordinary differential equation (ODE) network. In some embodiments, the method further includes an Ordinary Differential Equation (ODE) solver to form a neural ODE system. In some embodiments, the state transition is a mesenchymal-to-epithelial transition (MET), or epithelial-to-mesenchymal (EMT). In some embodiments, the calculating step further comprises the steps of approximating a time-varying derivative parametrized by a set of network weights and biases and integrating the time-varying derivative across the at least two time indices to calculate a possible population flow of the set of possible population flows. In some embodiments, the dynamic molecular program is a gene regulatory network. In some embodiments, the set of population flows which are unrealistic comprise energy inefficient biological pathways.

[0037] In some embodiments, the method further comprises the step of limiting a magnitude of the time-varying derivative across the at least two time indices to a predetermined threshold. In some embodiments, the method further comprises the steps of withholding a third static snapshot of the population of cells at a corresponding third time index between the at least two time indices from the calculating step, and adjusting one or more weights of the neural network to incentivize approximation of a function inferred from the at least two static snapshots which intersects the third static snapshot.

[0038] Tn some embodiments, the set of at least two static snapshots comprise three-dimensional tumorsphere data. In some embodiments, the method further comprises the steps of identifying shared and trajectory-specific transcription factor programs using a time-lapsed causality analysis. In some embodiments, the at least two time indices of the at least two static snapshots are separated by at least 0.1 hour, 0.2 hour, 0.4 hour, 0.6 hour, 0.8 hour, 1 hour, 2 hours, 3 hours, 4 hours, 5 hours, 6 hours, 8 hours, 10 hours, 12 hours, 14 hours, 16 hours, 18 hours, 20 hours, 22 hours, 24 hours. In some embodiments, the method further comprises computing an epithelial- mesenchymal transition (EMT) signature score based on an average expression of genes known to play a role in an epithelial-mesenchymal process. In some embodiments, the method further comprises the step of calculating an interpolated snapshot of the cell population at a third time index based on the inferred population flow. [0039] The present invention provides methods for preventing a cell from undergoing a transition. The present invention also provides methods for directing a cell to undergo a transition. In either of these embodiments, the transition is a mesenchymal-to-epithelial transition (MET). In either of these embodiments, the method of preventing a cell from undergoing MET comprises contacting a cell with one or more modulators described elsewhere herein. In either of these embodiments, the one or more molecules in the MET transcriptional network are one or more selected from the group consisting estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAI1 (SNAI1), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1).

[0040] Aspects of the present invention relate to a method of directing one or more cells in a population of cells to a transition state, having the steps of obtaining a population of cells from the subject, identifying a target state transition for the population of cells to undergo, administering at least one modulator of at least one molecule within a MET transcriptional network of the population of cells thereby directing the population of cells to the targeted state transition. In some embodiments, the one or more molecules in the MET transcriptional network are one or more selected from the group consisting estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAI1 (SNAI1), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1). In some embodiments, the targeted state transition is energy optimal.

[0041] The present invention provides methods of preventing a cell from undergoing a transition. The present invention also provides methods for directing a cell to undergo a transition. In some embodiments, the transition is a mesenchymal-to-epithelial transition (MET). In some embodiments, the method of preventing or directing a cell from undergoing MET comprises contacting a cell with one or more modulators described elsewhere herein. In some embodiments, the one or more molecules in the MET transcriptional network are one or more selected from the group consisting estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAU (SNAU), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1). In some embodiments, the method comprises contacting the cell with one or more inhibitors of ESRRA, AHR, and CDH1. In some embodiments, the method comprises contacting the cell with one or more activators of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and SNAI2.

[0042] Modulators

[0043] In various embodiments, the present invention relates to modulators of one or more molecules involved in a cell transition. In some embodiments, the cell transition is a mesenchymal-to-epithelial transition (MET). In some embodiments, the modulator prevents a cell from undergoing MET. In some embodiments, the present invention comprises a composition comprising one or more modulators of one or more molecules in the MET transcriptional network. In some embodiments, the one or more molecules in the MET transcriptional network are one or more selected from the group consisting of estrogen related receptor alpha (ESRRA), a nucleic acid encoding ESRRA (e.g., mRNA encoding ESRRA, the ESRRA gene, etc.), aryl hydrocarbon receptor (AHR), a nucleic acid encoding AHR, aryl hydrocarbon receptor nuclear translocator (ARNT), a nucleic acid encoding ARNT, estrogen receptor 1 (ESRI), transcription factor Jun (JUN), a nucleic acid encoding JUN, androgen receptor (AR), a nucleic acid encoding AR, zinc finger E-box binding homeobox 1 (ZEB1), a nucleic acid encoding ZEB1, zinc finger protein SNAU (SNAU), a nucleic acid encoding SNAU, zinc finger protein SNAI2 (SNAI2), a nucleic acid encoding SNAI2, and cadherin 1 (CDH1).

[0044] In various embodiments, the modulator alters the amount of a protein, the turnover of a protein, the activity of a protein, the phosphorylation of a protein, the acetylation of a protein, the amount of an mRNA encoding a protein, the stability of an mRNA encoding a protein, the translation of an mRNA encoding a protein, the transcription of a gene encoding a protein, or a combination thereof.

[0045] In various embodiments, the modulator is an inhibitor. In various embodiments, the inhibitor decreases the amount of a protein, the stability of a protein, the activity of a protein, the phosphorylation of a protein, the acetylation of a protein, the amount of an mRNA encoding a protein, the stability of an mRNA encoding a protein, the translation of an mRNA encoding a protein, the transcription of a gene encoding a protein, or a combination thereof.

[0046] In various embodiments, the modulator is an activator. In various embodiments, the activator increases the amount of a protein, the stability of a protein, the activity of a protein, the phosphorylation of a protein, the acetylation of a protein, the amount of an mRNA encoding a protein, the stability of an mRNA encoding a protein, the translation of an mRNA encoding a protein, the transcription of a gene encoding a protein, or a combination thereof.

[0047] Inhibitors

[0048] In various embodiments, the present invention relates to modulators of one or more molecules involved in a cell transition (e.g. provides a composition for altering the MET of a cell). In certain embodiments, the composition inhibits one or more proteins or nucleic acids involved in MET. In various embodiments, the modulator is an inhibitor.

[0049] In some embodiments, the composition of the invention comprises an inhibitor of estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), cadherin 1 (CDH1), or a combination thereof. An inhibitor ESRRA, AHR, or CDH1 is any compound, molecule, or agent that reduces, inhibits, or prevents the function of an ESRRA, AHR, or CDH1. For example, an inhibitor of ESRRA, AHR, or CDH1 is any compound, molecule, or agent that decreases the amount, stability, or activity of ESRRA, AHR, or CDH1, decreases the amount, stability, or translation of an mRNA encoding ESRRA, AHR, or CDH1, decreases the transcription of a gene encoding ESRRA, AHR, or CDH1, or a combination thereof. In some embodiments, an inhibitor of ESRRA, AHR, or CDH1 comprises a nucleic acid, a peptide, a small molecule, a siRNA, miRNA, shRNA, a ribozyme, an anti-sense nucleic acid, an antagonist, an inverse agonist, an aptamer, a peptidomimetic, or any combination thereof.

[0050] In various embodiments, the inhibitor is a small molecule. When the inhibitor is a small molecule, a small molecule may be obtained using standard methods known to the skilled artisan. Such methods include chemical organic synthesis or biological means. Biological means include purification from a biological source, recombinant synthesis and in vitro translation systems, using methods well known in the art. In some embodiments, a small molecule inhibitor of the invention comprises an organic molecule, inorganic molecule, biomolecule, synthetic molecule, and the like.

[0051] Combinatorial libraries of molecularly diverse chemical compounds potentially useful in treating a variety of diseases and conditions are well known in the art as are method of making the libraries. The method may use a variety of techniques well-known to the skilled artisan including solid phase synthesis, solution methods, parallel synthesis of single compounds, synthesis of chemical mixtures, rigid core structures, flexible linear sequences, deconvolution strategies, tagging techniques, and generating unbiased molecular landscapes for lead discovery vs. biased structures for lead development.

[0052] In a general method for small library synthesis, an activated core molecule is condensed with a number of building blocks, resulting in a combinatorial library of covalently linked, corebuilding block ensembles. The shape and rigidity of the core determines the orientation of the building blocks in shape space. The libraries can be biased by changing the core, linkage, or building blocks to target a characterized biological structure (“focused libraries”) or synthesized with less structural bias using flexible cores.

[0053] The small molecule and small molecule compounds described herein may be present as salts even if salts are not depicted and it is understood that the invention embraces all salts and solvates of the inhibitors depicted here, as well as the non-salt and non-solvate form of the inhibitors, as is well understood by the skilled artisan. In some embodiments, the salts of the inhibitors of the invention are pharmaceutically acceptable salts.

[0054] Where tautomeric forms may be present for any of the inhibitors described herein, each and every tautomeric form is intended to be included in the present invention, even though only one or some of the tautomeric forms may be explicitly depicted. For example, when a 2- hydroxypyridyl moiety is depicted, the corresponding 2-pyridone tautomer is also intended.

[0055] The invention also includes any or all of the stereochemical forms, including any enantiomeric or diastereomeric forms of the inhibitors described. The recitation of the structure or name herein is intended to embrace all possible stereoisomers of inhibitors depicted. All forms of the inhibitors are also embraced by the invention, such as crystalline or non-crystalline forms of the inhibitors. Compositions comprising an inhibitor of the invention are also intended, such as a composition of substantially pure inhibitor, including a specific stereochemical form thereof, or a composition comprising mixtures of inhibitors of the invention in any ratio, including two or more stereochemical forms, such as in a racemic or non-racemic mixture.

[0056] In some embodiments, the small molecule inhibitor of the invention comprises an analog or derivative of an inhibitor described herein.

[0057] In some instances, small molecule inhibitors described herein are derivatized/analoged as is well known in the art of combinatorial and medicinal chemistry. The analogs or derivatives can be prepared by adding and/or substituting functional groups at various locations. As such, the small molecules described herein can be converted into derivatives/analogs using well known chemical synthesis procedures. For example, all of the hydrogen atoms or substituents can be selectively modified to generate new analogs. Also, the linking atoms or groups can be modified into longer or shorter linkers with carbon backbones or hetero atoms. Also, the ring groups can be changed so as to have a different number of atoms in the ring and/or to include hetero atoms. Moreover, aromatics can be converted to cyclic rings, and vice versa. For example, the rings may be from 5-7 atoms, and may be homocycles or heterocycles.

[0058] As used herein, the term “analog,” “analogue,” or “derivative” is meant to refer to a chemical compound or molecule made from a parent compound or molecule by one or more chemical reactions. As such, an analog can be a structure having a structure similar to that of the small molecule inhibitors described herein or can be based on a scaffold of a small molecule inhibitor described herein, but differing from it in respect to certain components or structural makeup, which may have a similar or opposite action metabolically. An analog or derivative of any of a small molecule inhibitor in accordance with the present invention can be used to reduce skin pigmentation.

[0059] In some embodiments, the small molecule inhibitors described herein can independently be derivatized/analoged by modifying hydrogen groups independently from each other into other substituents. That is, each atom on each molecule can be independently modified with respect to the other atoms on the same molecule. Any traditional modification for producing a derivative/analog can be used. For example, the atoms and substituents can be independently comprised of hydrogen, an alkyl, aliphatic, straight chain aliphatic, aliphatic having a chain hetero atom, branched aliphatic, substituted aliphatic, cyclic aliphatic, heterocyclic aliphatic having one or more hetero atoms, aromatic, heteroaromatic, polyaromatic, poly-amino acids, peptides, polypeptides, combinations thereof, halogens, halo-substituted aliphatics, and the like. Additionally, any ring group on a compound can be derivatized to increase and/or decrease ring size as well as change the backbone atoms to carbon atoms or hetero atoms.

[0060] Nucleic acid inhibitors

[0061] In other related aspects, the invention relates to modulators of one or more molecules involved in a cell transition (e.g. an isolated nucleic acid). In some instances, the inhibitor is an siRNA, shRNA miRNA, or anti-sense oligonucleotide, which inhibits ESRRA, AHR, CDH1, or a combination thereof. In some embodiments, the nucleic acid comprises a promoter/regulatory sequence such that the nucleic acid is capable of directing expression of the nucleic acid. Thus, the invention encompasses expression vectors and methods for the introduction of exogenous DNA into cells with concomitant expression of the exogenous DNA in the cells such as those described, for example, in Sambrook et al. (2012, Molecular Cloning: A Laboratory Manual, Cold Spring Harbor Laboratory, New York), and in Ausubel et al. (1997, Current Protocols in Molecular Biology, John Wiley & Sons, New York) and as described elsewhere herein.

[0062] In some embodiments, siRNA, shRNA, miRNA, or an anti-sense oligonucleotide is used to decrease the level of ESRRA, AHR, CDH1 protein, or a combination thereof. RNA interference (RNAi) is a phenomenon in which the introduction of double-stranded RNA (dsRNA) into a diverse range of organisms and cell types causes degradation of the complementary mRNA. In the cell, long dsRNAs are cleaved into short 21-25 nucleotide small interfering RNAs, or siRNAs, by a ribonuclease known as Dicer. The siRNAs subsequently assemble with protein components into an RNA-induced silencing complex (RISC), unwinding in the process. Activated RISC then binds to complementary transcript by base pairing interactions between the siRNA anti-sense strand and the mRNA. The bound mRNA is cleaved and sequence specific degradation of mRNA results in gene silencing. See, for example, U.S. Patent No. 6,506,559; Fire et al., 1998, Nature 391(19):306-311; Timmons et al., 1998, Nature 395:854; Montgomery et al., 1998, TIG 14 (7):255-258; David R. Engelke, Ed., RNA Interference (RNAi) Nuts & Bolts of RNAi Technology, DNA Press, Eagleville, PA (2003); and Gregory J. Hannon, Ed., RNAi A Guide to Gene Silencing, Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY (2003). Soutschek et al. (2004, Nature 432:173-178) describe a chemical modification to siRNAs that aids in intravenous systemic delivery. Optimizing siRNAs involves consideration of overall G/C content, C/T content at the termini, Tm and the nucleotide content of the 3’ overhang. See, for instance, Schwartz et al., 2003, Cell, 115: 199-208 and Khvorova et al., 2003, Cell 115:209-216. Therefore, the present invention also includes methods of decreasing levels of ESRRA, AHR, and CDH1 using RNAi technology.

[0063] In another aspect, the invention includes a vector comprising an siRNA, miRNA, or antisense oligonucleotide. In some embodiments, the siRNA, miRNA, or anti-sense polynucleotide is capable of inhibiting the expression of a target polypeptide, wherein the target polypeptide is selected from the group consisting of ESRRA, AHR, and CDH1. The incorporation of a desired polynucleotide into a vector and the choice of vectors is well-known in the art as described in, for example, Sambrook et al. (2012), and in Ausubel et al. (1997), and elsewhere herein.

[0064] In certain embodiments, the expression vectors described herein encode a short hairpin RNA (shRNA) inhibitor. shRNA inhibitors are well known in the art and are directed against the mRNA of a target, thereby decreasing the expression of the target. In certain embodiments, the encoded shRNA is expressed by a cell, and is then processed into siRNA. For example, in certain instances, the cell possesses native enzymes (e.g., dicer) that cleaves the shRNA to form siRNA. [0065] The siRNA, miRNA, shRNA, or anti-sense oligonucleotide can be cloned into a number of types of vectors as described elsewhere herein. For expression of the siRNA, miRNA, or antisense polynucleotide, at least one module in each promoter functions to position the start site for RNA synthesis.

[0066] In order to assess the expression of the siRNA, miRNA, shRNA, or anti-sense polynucleotide, the expression vector to be introduced into a cell can also contain either a selectable marker gene or a reporter gene or both to facilitate identification and selection of expressing cells from the population of cells sought to be transfected or infected using a viral vector. In other embodiments, the selectable marker may be carried on a separate piece of DNA and used in a co-transfection procedure. Both selectable markers and reporter genes may be flanked with appropriate regulatory sequences to enable expression in the host cells. Useful selectable markers are known in the art and include, for example, antibiotic-resistance genes, such as neomycin resistance and the like.

[0067] Therefore, in another aspect, the invention relates to a vector, comprising the nucleotide sequence of the invention or the construct of the invention. The choice of the vector will depend on the host cell in which it is to be subsequently introduced. In a particular embodiment, the vector of the invention is an expression vector. Suitable host cells include a wide variety of prokaryotic and eukaryotic host cells. In specific embodiments, the expression vector is selected from the group consisting of a viral vector, a bacterial vector and a mammalian cell vector. Prokaryote- and/or eukaryote-vector based systems can be employed for use with the present invention to produce polynucleotides, or their cognate polypeptides. Many such systems are commercially and widely available.

[0068] Further, the expression vector may be provided to a cell in the form of a viral vector. Viral vector technology is well known in the art and is described, for example, in Sambrook et al. (2012), and in Ausubel et al. (1997), and in other virology and molecular biology manuals. Viruses, which are useful as vectors include, but are not limited to, retroviruses, adenoviruses, adeno-associated viruses, herpes viruses, and lentiviruses. In general, a suitable vector contains an origin of replication functional in at least one organism, a promoter sequence, convenient restriction endonuclease sites, and one or more selectable markers. (See, e.g., WO 01/96584;

WO 01/29058; and U.S. Pat. No. 6,326,193.

[0069] By way of illustration, the vector in which the nucleic acid sequence is introduced can be a plasmid which is or is not integrated in the genome of a host cell when it is introduced in the cell. Illustrative, non-limiting examples of vectors in which the nucleotide sequence of the invention or the gene construct of the invention can be inserted include a tet-on inducible vector for expression in eukaryote cells.

[0070] The vector may be obtained by conventional methods known by persons skilled in the art (Sambrook et al., 2012). In a particular embodiment, the vector is a vector useful for transforming animal cells.

[0071] In some embodiments, the recombinant expression vectors may also contain nucleic acid molecules which encode a peptide or peptidomimetic inhibitor of invention, described elsewhere herein.

[0072] A promoter may be one naturally associated with a gene or polynucleotide sequence, as may be obtained by isolating the 5' non-coding sequences located upstream of the coding segment and/or exon. Such a promoter can be referred to as “endogenous.” Similarly, an enhancer may be one naturally associated with a polynucleotide sequence, located either downstream or upstream of that sequence. Alternatively, certain advantages will be gained by positioning the coding polynucleotide segment under the control of a recombinant or heterologous promoter, which refers to a promoter that is not normally associated with a polynucleotide sequence in its natural environment. A recombinant or heterologous enhancer refers also to an enhancer not normally associated with a polynucleotide sequence in its natural environment. Such promoters or enhancers may include promoters or enhancers of other genes, and promoters or enhancers isolated from any other prokaryotic, viral, or eukaryotic cell, and promoters or enhancers not “naturally occurring,” i.e., containing different elements of different transcriptional regulatory regions, and/or mutations that alter expression. In addition to producing nucleic acid sequences of promoters and enhancers synthetically, sequences may be produced using recombinant cloning and/or nucleic acid amplification technology, including PCR™, in connection with the compositions disclosed herein (U.S. Patent 4,683,202, U.S. Patent 5,928,906). Furthermore, it is contemplated the control sequences that direct transcription and/or expression of sequences within non-nuclear organelles such as mitochondria, chloroplasts, and the like, can be employed as well.

[0073] Naturally, it will be important to employ a promoter and/or enhancer that effectively directs the expression of the DNA segment in the cell type, organelle, and organism chosen for expression. Those of skill in the art of molecular biology generally know how to use promoters, enhancers, and cell type combinations for protein expression, for example, see Sambrook et al. (2012). The promoters employed may be constitutive, tissue-specific, inducible, and/or useful under the appropriate conditions to direct high-level expression of the introduced DNA segment, such as is advantageous in the large-scale production of recombinant proteins and/or peptides. The promoter may be heterologous or endogenous.

[0074] The recombinant expression vectors may also contain a selectable marker gene which facilitates the selection of transformed or transfected host cells. Suitable selectable marker genes are genes encoding proteins such as G418 and hygromycin which confer resistance to certain drugs, P-galactosidase, chloramphenicol acetyltransferase, firefly luciferase, or an immunoglobulin or portion thereof such as the Fc portion of an immunoglobulin, for example, IgG. The selectable markers may be introduced on a separate vector from the nucleic acid of interest.

[0075] Following the generation of the siRNA, miRNA, or anti-sense oligonucleotide, a skilled artisan will understand that the siRNA, miRNA, or anti-sense oligonucleotide will have certain characteristics that can be modified to improve the siRNA, miRNA, or anti-sense oligonucleotide as a therapeutic compound. Therefore, the siRNA, miRNA, or anti-sense oligonucleotide may be further designed to resist degradation by modifying it to include phosphorothioate, or other linkages, methyl phosphonate, sulfone, sulfate, ketyl, phosphorodithioate, phosphoramidate, phosphate esters, and the like (see, e.g., Agrwal et al., 1987, Tetrahedron Lett. 28:3539-3542; Stec et al., 1985 Tetrahedron Lett. 26:2191-2194; Moody et al., 1989 Nucleic Acids Res. 12:4769-4782; Eckstein, 1989 Trends Biol. Sci. 14:97-100; Stein, In: Oligodeoxynucleotides. Anti-sense Inhibitors of Gene Expression, Cohen, ed., Macmillan Press, London, pp. 97-117 (1989)).

[0076] Any oligonucleotide may be further modified to increase its stability in vivo. Possible modifications include, but are not limited to, the addition of flanking sequences at the 5' and/or 3' ends; the use of phosphorothioate or 2' O-methyl rather than phosphodiester linkages in the backbone; and/or the inclusion of nontraditional bases such as inosine, queosine, and wybutosine and the like, as well as acetyl- methyl-, thio- and other modified forms of adenine, cytidine, guanine, thymine, and uridine.

[0077] In some embodiments of the invention, an anti-sense nucleic acid sequence which is expressed by a plasmid vector is used to inhibit ESRRA, AHR, and/or CDH1 protein expression. The anti-sense expressing vector is used to transfect a mammalian cell or the mammal itself, thereby causing reduced endogenous expression of ESRRA, AHR, and/or CDH1 protein.

[0078] Anti-sense molecules and their use for inhibiting gene expression are well known in the art (see, e.g., Cohen, 1989, In: Oligodeoxyribonucleotides, Antisense Inhibitors of Gene Expression, CRC Press). Anti-sense nucleic acids are DNA or RNA molecules that are complementary, as that term is defined elsewhere herein, to at least a portion of a specific mRNA molecule (Weintraub, 1990, Scientific American 262:40). In the cell, anti-sense nucleic acids hybridize to the corresponding mRNA, forming a double- stranded molecule thereby inhibiting the translation of genes.

[0079] The use of anti-sense methods to inhibit the translation of genes is known in the art, and is described, for example, in Marcus-Sakura (1988, Anal. Biochem. 172:289). Such anti-sense molecules may be provided to the cell via genetic expression using DNA encoding the anti-sense molecule as taught by Inoue, 1993, U.S. Patent No. 5,190,931.

[0080] Alternatively, anti-sense molecules of the invention may be made synthetically and then provided to the cell. In some embodiments, anti-sense oligomers may have between about 10 to about 30 nucleotides. Tn some embodiments, anti-sense oligomers may have about 15 nucleotides. In some embodiments, anti-sense oligomers having 10-30 nucleotides are easily synthesized and introduced into a target cell. Synthetic anti-sense molecules contemplated by the invention include oligonucleotide derivatives known in the art which have improved biological activity compared to unmodified oligonucleotides (see U.S. Patent No. 5,023,243).

[0081] In some embodiments of the invention, a ribozyme is used to inhibit ESRRA, AHR, and/or CDH1 protein expression. Ribozymes useful for inhibiting the expression of a target molecule may be designed by incorporating target sequences into the basic ribozyme structure which are complementary, for example, to the mRNA sequence encoding ESRRA, AHR, and/or CDH1. Ribozymes targeting ESRRA, AHR, and/or CDH1 may be synthesized using commercially available reagents (Applied Biosystems, Inc., Foster City, CA) or they may be genetically expressed from DNA encoding them.

[0082] In some embodiments, the inhibitor of ESRRA, AHR, and/or CDH1 may comprise one or more components of a CRISPR-Cas system, where a guide RNA (gRNA) targeted to a gene encoding ESRRA, AHR, and/or CDH1, and a CRISPR-associated (Cas) peptide form a complex to induce mutations within the targeted gene. In some embodiments, the inhibitor comprises a gRNA or a nucleic acid molecule encoding a gRNA. In some embodiments, the inhibitor comprises a Cas peptide or a nucleic acid molecule encoding a Cas peptide.

[0083] Polypeptide inhibitors

[0084] In other related aspects, the invention relates to modulators of one or more molecules involved in a cell transition (e.g. an isolated peptide inhibitor that inhibits ESRRA, AHR, and/or CDH1). For example, in some embodiments, the peptide inhibitor of the invention inhibits ESRRA, AHR, and/or CDH1 directly by binding to ESRRA, AHR, and/or CDH1 thereby preventing the normal functional activity of ESRRA, AHR, and/or CDHI. In another embodiment, the peptide inhibitor of the invention inhibits ESRRA, AHR, and/or CDHI by competing with endogenous ESRRA, AHR, and/or CDHI. In yet another embodiment, the peptide inhibitor of the invention inhibits the activity of ESRRA, AHR, and/or CDHI by acting as a transdominant negative mutant.

[0085] The variants of the polypeptides according to the present invention may be (i) one in which one or more of the amino acid residues are substituted with a conserved or non-conserved amino acid residue and such substituted amino acid residue may or may not be one encoded by the genetic code, (ii) one in which there are one or more modified amino acid residues, e.g., residues that are modified by the attachment of substituent groups, (iii) one in which the polypeptide is an alternative splice variant of the polypeptide of the present invention, (iv) fragments of the polypeptides and/or (v) one in which the polypeptide is fused with another polypeptide, such as a leader or secretory sequence or a sequence which is employed for purification (for example, His-tag) or for detection (for example, Sv5 epitope tag). The fragments include polypeptides generated via proteolytic cleavage (including multi-site proteolysis) of an original sequence. Variants may be post-translationally, or chemically modified. Such variants are deemed to be within the scope of those skilled in the art from the teaching herein.

[0086] Antibody inhibitors

[0087] The invention also relates to modulators of one or more molecules involved in a cell transition (e.g. an inhibitor of ESRRA, AHR, and/or CDH1 comprising an antibody, or antibody fragment, specific for ESRRA, AHR, and/or CDH1). That is, the antibody can inhibit ESRRA, AHR, and/or CDH1 to provide a beneficial effect.

[0088] The antibodies may be intact monoclonal or polyclonal antibodies, and immunologically active fragments (e.g., a Fab or (Fab)2 fragment), an antibody heavy chain, an antibody light chain, humanized antibodies, a genetically engineered single chain Fv molecule (Ladner et al, U.S. Pat. No. 4,946,778), or a chimeric antibody, for example, an antibody which contains the binding specificity of a murine antibody, but in which the remaining portions are of human origin. Antibodies including monoclonal and polyclonal antibodies, fragments and chimeras, may be prepared using methods known to those skilled in the art.

[0089] Antibodies can be prepared using intact polypeptides or fragments containing an immunizing antigen of interest. The polypeptide or oligopeptide used to immunize an animal may be obtained from the translation of RNA or synthesized chemically and can be conjugated to a carrier protein, if desired. Suitable carriers that may be chemically coupled to peptides include bovine serum albumin and thyroglobulin, keyhole limpet hemocyanin. The coupled polypeptide may then be used to immunize the animal (e.g., a mouse, a rat, or a rabbit).

[0090] Activators [0091 ] Tn various embodiments, the present invention relates to modulators of one or more molecules involved in a cell transition (e.g. a composition for altering the MET of a cell). In certain embodiments, the composition activates one or more proteins or nucleic acids involved in MET. In various embodiments, the modulator is an activator.

[0092] In some embodiments, the composition of the invention comprises an activator of aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAI1 (SNAI1), zinc finger protein SNAI2 (SNAI2), or a combination thereof. An activator of ARNT, ESRI, JUN, AR, ZEB1, SNAU, or SNAI2 is any compound, molecule, or agent that increases the function or activity of an ARNT, ESRI, JUN, AR, ZEB1, SNAI1, or SNAI2. For example, an activator of ARNT, ESRI, JUN, AR, ZEB1, SNAU, or SNAI2 is any compound, molecule, or agent that increases the amount, stability, or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAU, or SNAI2, increases the amount, stability, or translation of an mRNA encoding ARNT, ESRI, JUN, AR, ZEB1, SNAI1, or SNAI2, increases the transcription of a gene encoding ARNT, ESRI, JUN, AR, ZEB1, SNAI1, or SNAI2, or a combination thereof. In some embodiments, an activator of ARNT, ESRI, JUN, AR, ZEB1, SNAU, or SNAI2 comprises a nucleic acid, a peptide, a small molecule, a siRNA, miRNA, shRNA, a ribozyme, an anti-sense nucleic acid, an agonist, a partial agonist, an aptamer, a peptidomimetic, or any combination thereof.

[0093] It will be understood by one skilled in the art, based upon the disclosure provided herein, that an increase in the level of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 encompasses the increase in ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 expression, including transcription, translation, or both. The skilled artisan will also appreciate, once armed with the teachings of the present invention, that an increase in the level of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 includes an increase in ARNT, ESRI, JUN, AR, ZEB 1, SNAII, and/or SNAI2 activity (e.g., enzymatic activity, substrate binding activity, etc.). Thus, increasing the level or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 includes, but is not limited to, increasing the amount of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 protein, and increasing transcription, translation, or both, of a nucleic acid encoding ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2; and it also includes increasing any activity of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 protein as well. [0094] The increased level or activity of ARNT, ESRI , JUN, AR, ZEB1 , SNAT1, and/or SNAT2 can be assessed using a wide variety of methods well-known in the art or to be developed in the future. That is, the routineer would appreciate, based upon the disclosure provided herein, that increasing the level or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 can be readily assessed using methods that assess the level of a nucleic acid encoding ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 (e.g., mRNA), the level of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 protein, and/or the level of ARNT, ESRI, JUN, AR, ZEB1, SNAIl, and/or SNAI2 activity in a biological sample obtained from a subject.

[0095] One of skill in the art will realize that in addition to activating ARNT, ESRI, JUN, AR, ZEB1, SNAIl, and/or SNAI2 directly, diminishing the amount or activity of a molecule that itself diminishes the amount or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 can serve to increase the amount or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2. Thus, an ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 activator can include, but should not be construed as being limited to, a nucleic acid, a peptide, a small molecule, a siRNA, miRNA, shRNA, a ribozyme, an anti-sense nucleic acid, an agonist, a partial agonist, an aptamer, a peptidomimetic, or any combination thereof.

[0096] Small molecule activators

[0097] In various embodiments, the present invention relates to modulators of one or more molecules involved in a cell transition (e.g. an activator that a small molecule). When the activator is a small molecule, a small molecule may be obtained using standard methods known to the skilled artisan. Such methods include chemical organic synthesis or biological means. Biological means include purification from a biological source, recombinant synthesis and in vitro translation systems, using methods well known in the art. In some embodiments, a small molecule activator of the invention comprises an organic molecule, inorganic molecule, biomolecule, synthetic molecule, and the like.

[0098] Combinatorial libraries of molecularly diverse chemical compounds potentially useful in treating a variety of diseases and conditions are well known in the art as are method of making the libraries. The method may use a variety of techniques well-known to the skilled artisan including solid phase synthesis, solution methods, parallel synthesis of single compounds, synthesis of chemical mixtures, rigid core structures, flexible linear sequences, deconvolution strategies, tagging techniques, and generating unbiased molecular landscapes for lead discovery vs. biased structures for lead development.

[0099] In a general method for small library synthesis, an activated core molecule is condensed with a number of building blocks, resulting in a combinatorial library of covalently linked, corebuilding block ensembles. The shape and rigidity of the core determines the orientation of the building blocks in shape space. The libraries can be biased by changing the core, linkage, or building blocks to target a characterized biological structure (“focused libraries”) or synthesized with less structural bias using flexible cores.

[0100] The small molecule and small molecule compounds described herein may be present as salts even if salts are not depicted and it is understood that the invention embraces all salts and solvates of the activators depicted here, as well as the non-salt and non-solvate form of the activators, as is well understood by the skilled artisan. In some embodiments, the salts of the activators of the invention are pharmaceutically acceptable salts.

[0101] Where tautomeric forms may be present for any of the activators described herein, each and every tautomeric form is intended to be included in the present invention, even though only one or some of the tautomeric forms may be explicitly depicted. For example, when a 2- hydroxypyridyl moiety is depicted, the corresponding 2-pyridone tautomer is also intended. [0102] The invention also includes any or all of the stereochemical forms, including any enantiomeric or diastereomeric forms of the activators described. The recitation of the structure or name herein is intended to embrace all possible stereoisomers of activators depicted. All forms of the activators are also embraced by the invention, such as crystalline or non-crystalline forms of the activators. Compositions comprising an activator of the invention are also intended, such as a composition of substantially pure activator, including a specific stereochemical form thereof, or a composition comprising mixtures of activators of the invention in any ratio, including two or more stereochemical forms, such as in a racemic or non-racemic mixture.

[0103] In some embodiments, the small molecule activator of the invention comprises an analog or derivative of an activator described herein.

[0104] In some embodiments, the small molecules described herein are candidates for derivatization. As such, in certain instances, the analogs of the small molecules described herein that have modulated potency, selectivity, and solubility are included herein and provide useful leads for drug discovery and drug development. Thus, in certain instances, during optimization new analogs are designed considering issues of drug delivery, metabolism, novelty, and safety.

[0105] In some instances, small molecule activators described herein are derivatized/analoged as is well known in the art of combinatorial and medicinal chemistry. The analogs or derivatives can be prepared by adding and/or substituting functional groups at various locations. As such, the small molecules described herein can be converted into derivatives/analogs using well known chemical synthesis procedures. For example, all of the hydrogen atoms or substituents can be selectively modified to generate new analogs. Also, the linking atoms or groups can be modified into longer or shorter linkers with carbon backbones or hetero atoms. Also, the ring groups can be changed so as to have a different number of atoms in the ring and/or to include hetero atoms. Moreover, aromatics can be converted to cyclic rings, and vice versa. For example, the rings may be from 5-7 atoms, and may be homocycles or heterocycles.

[0106] As used herein, the term “analog,” “analogue,” or “derivative” is meant to refer to a chemical compound or molecule made from a parent compound or molecule by one or more chemical reactions. As such, an analog can be a structure having a structure similar to that of the small molecule activators described herein or can be based on a scaffold of a small molecule activator described herein, but differing from it in respect to certain components or structural makeup, which may have a similar or opposite action metabolically. An analog or derivative of any of a small molecule activator in accordance with the present invention can be used to reduce skin pigmentation.

[0107] In some embodiments, the small molecule activators described herein can independently be derivatized/analoged by modifying hydrogen groups independently from each other into other substituents. That is, each atom on each molecule can be independently modified with respect to the other atoms on the same molecule. Any traditional modification for producing a derivative/analog can be used. For example, the atoms and substituents can be independently comprised of hydrogen, an alkyl, aliphatic, straight chain aliphatic, aliphatic having a chain hetero atom, branched aliphatic, substituted aliphatic, cyclic aliphatic, heterocyclic aliphatic having one or more hetero atoms, aromatic, heteroaromatic, polyaromatic, poly-amino acids, peptides, polypeptides, combinations thereof, halogens, halo-substituted aliphatics, and the like. Additionally, any ring group on a compound can be derivatized to increase and/or increase ring size as well as change the backbone atoms to carbon atoms or hetero atoms.

[0108] Nucleic acid activators

[0109] In other related aspects, the invention relates to modulators of one or more molecules involved in a cell transition (e.g. an isolated nucleic acid acting as an activator). In some instances, the activator is an mRNA, ssDNA, dsDNA, or combination thereof, which encodes ARNT, ESRI, JUN, AR, ZEB1, SNAI1, SNAI2, or a combination thereof. In some embodiments, the nucleic acid comprises a promoter/regulatory sequence such that the nucleic acid is capable of directing expression of the nucleic acid. Thus, the invention encompasses expression vectors and methods for the introduction of exogenous nucleic acids into cells with concomitant expression of the exogenous nucleic acids in the cells such as those described, for example, in Sambrook et al. (2012, Molecular Cloning: A Laboratory Manual, Cold Spring Harbor Laboratory, New York), and in Ausubel et al. (1997, Current Protocols in Molecular Biology, John Wiley & Sons, New York) and as described elsewhere herein.

[0110] In some embodiments, mRNA, ssDNA, dsDNA, or a combination thereof is used to increase the level of ARNT, ESRI, JUN, AR, ZEB1, SNAI1, SNAI2 protein, or a combination thereof. Said nucleic acid molecule can directly be translated into protein, can be transcribed into mRNA for translation into protein, or inserted into a genome to replace or supplement the native gene.

[0111] In some aspects, the invention relates to a vector, comprising the nucleotide sequence of the invention or the construct of the invention. The choice of the vector will depend on the host cell in which it is to be subsequently introduced. In a particular embodiment, the vector of the invention is an expression vector. Suitable host cells include a wide variety of prokaryotic and eukaryotic host cells. In specific embodiments, the expression vector is selected from the group consisting of a viral vector, a bacterial vector and a mammalian cell vector. Prokaryote- and/or eukaryote-vector based systems can be employed for use with the present invention to produce polynucleotides, or their cognate polypeptides. Many such systems are commercially and widely available.

[0112] Further, the expression vector may be provided to a cell in the form of a viral vector. Viral vector technology is well known in the art and is described, for example, in Sambrook et al. (2012), and in Ausubel et al. (1997), and in other virology and molecular biology manuals. Viruses, which are useful as vectors include, but are not limited to, retroviruses, adenoviruses, adeno-associated viruses, herpes viruses, and lentiviruses. In general, a suitable vector contains an origin of replication functional in at least one organism, a promoter sequence, convenient restriction endonuclease sites, and one or more selectable markers. (See, e.g., WO 01/96584; WO 01/29058; and U.S. Pat. No. 6,326,193.

[0113] By way of illustration, the vector in which the nucleic acid sequence is introduced can be a plasmid which is or is not integrated in the genome of a host cell when it is introduced in the cell. Illustrative, non-limiting examples of vectors in which the nucleotide sequence of the invention or the gene construct of the invention can be inserted include a tet-on inducible vector for expression in eukaryote cells.

[0114] The vector may be obtained by conventional methods known by persons skilled in the art (Sambrook et al., 2012). In a particular embodiment, the vector is a vector useful for transforming animal cells.

[0115] In some embodiments, the recombinant expression vectors may also contain nucleic acid molecules which encode a peptide or peptidomimetic activator of invention, described elsewhere herein.

[0116] A promoter may be one naturally associated with a gene or polynucleotide sequence, as may be obtained by isolating the 5' non-coding sequences located upstream of the coding segment and/or exon. Such a promoter can be referred to as “endogenous.” Similarly, an enhancer may be one naturally associated with a polynucleotide sequence, located either downstream or upstream of that sequence. Alternatively, certain advantages will be gained by positioning the coding polynucleotide segment under the control of a recombinant or heterologous promoter, which refers to a promoter that is not normally associated with a polynucleotide sequence in its natural environment. A recombinant or heterologous enhancer refers also to an enhancer not normally associated with a polynucleotide sequence in its natural environment. Such promoters or enhancers may include promoters or enhancers of other genes, and promoters or enhancers isolated from any other prokaryotic, viral, or eukaryotic cell, and promoters or enhancers not “naturally occurring,” i.e., containing different elements of different transcriptional regulatory regions, and/or mutations that alter expression. In addition to producing nucleic acid sequences of promoters and enhancers synthetically, sequences may be produced using recombinant cloning and/or nucleic acid amplification technology, including PCR, in connection with the compositions disclosed herein (U.S. Patent 4,683,202, U.S. Patent 5,928,906). Furthermore, it is contemplated the control sequences that direct transcription and/or expression of sequences within non-nuclear organelles such as mitochondria, chloroplasts, and the like, can be employed as well.

[0117] Naturally, it will be important to employ a promoter and/or enhancer that effectively directs the expression of the DNA segment in the cell type, organelle, and organism chosen for expression. Those of skill in the art of molecular biology generally know how to use promoters, enhancers, and cell type combinations for protein expression, for example, see Sambrook et al. (2012). The promoters employed may be constitutive, tissue-specific, inducible, and/or useful under the appropriate conditions to direct high-level expression of the introduced DNA segment, such as is advantageous in the large-scale production of recombinant proteins and/or peptides. The promoter may be heterologous or endogenous.

[0118] The recombinant expression vectors may also contain a selectable marker gene which facilitates the selection of transformed or transfected host cells. Suitable selectable marker genes are genes encoding proteins such as G418 and hygromycin which confer resistance to certain drugs, P-galactosidase, chloramphenicol acetyltransferase, firefly luciferase, or an immunoglobulin or portion thereof such as the Fc portion of an immunoglobulin, for example, IgG. The selectable markers may be introduced on a separate vector from the nucleic acid of interest.

[0119] One of skill in the art will realize that diminishing the amount or activity of a molecule that itself diminishes the amount or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2 can serve to increase the amount or activity of ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2. RNA interference (RNAi) is a phenomenon in which the introduction of doublestranded RNA (dsRNA) into a diverse range of organisms and cell types causes degradation of the complementary mRNA. In the cell, long dsRNAs are cleaved into short 21-25 nucleotide small interfering RNAs, or siRNAs, by a ribonuclease known as Dicer. The siRNAs subsequently assemble with protein components into an RNA-induced silencing complex (RISC), unwinding in the process. Activated RISC then binds to complementary transcript by base pairing interactions between the siRNA anti-sense strand and the mRNA. The bound mRNA is cleaved and sequence specific degradation of mRNA results in gene silencing. See, for example, U.S. Patent No. 6,506,559; Fire et al., 1998, Nature 391 (19):306-311 ; Timmons et al., 1998, Nature 395:854; Montgomery et al., 1998, TIG 14 (7):255-258; David R. Engelke, Ed., RNA Interference (RNAi) Nuts & Bolts of RNAi Technology, DNA Press, Eagleville, PA (2003); and Gregory J. Hannon, Ed., RNAi A Guide to Gene Silencing, Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY (2003). Soutschek et al. (2004, Nature 432: 173-178) describe a chemical modification to siRNAs that aids in intravenous systemic delivery. Optimizing siRNAs involves consideration of overall G/C content, C/T content at the termini, Tm and the nucleotide content of the 3’ overhang. See, for instance, Schwartz et al., 2003, Cell, 115: 199-208 and Khvorova et al., 2003, Cell 115:209-216. Therefore, the present invention also includes methods of increasing levels of ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2using RNAi technology.

[0120] In another aspect, the invention includes a vector comprising an siRNA, miRNA, or antisense oligonucleotide. In some embodiments, the siRNA, miRNA, or anti-sense polynucleotide is capable of inhibiting the expression of a target polypeptide, wherein the target polypeptide is an inhibitor of ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2. The incorporation of a desired polynucleotide into a vector and the choice of vectors is well-known in the art as described in, for example, Sambrook et al. (2012), and in Ausubel et al. (1997), and elsewhere herein.

[0121] In certain embodiments, the expression vectors described herein encode a short hairpin RNA (shRNA) activator. shRNA activators are well known in the art and are directed against the mRNA of a target, thereby decreasing the expression of the target. In certain embodiments, the encoded shRNA is expressed by a cell, and is then processed into siRNA. For example, in certain instances, the cell possesses native enzymes (e.g., dicer) that cleave the shRNA to form siRNA. [0122] The siRNA, miRNA, shRNA, or anti-sense oligonucleotide can be cloned into a number of types of vectors as described elsewhere herein. For expression of the siRNA, miRNA, or antisense polynucleotide, at least one module in each promoter functions to position the start site for RNA synthesis.

[0123] In order to assess the expression of the siRNA, miRNA, shRNA, or anti-sense polynucleotide, the expression vector to be introduced into a cell can also contain either a selectable marker gene or a reporter gene or both to facilitate identification and selection of expressing cells from the population of cells sought to be transfected or infected using a viral vector. Tn other embodiments, the selectable marker may be carried on a separate piece of DNA and used in a co-transfection procedure. Both selectable markers and reporter genes may be flanked with appropriate regulatory sequences to enable expression in the host cells. Useful selectable markers are known in the art and include, for example, antibiotic-resistance genes, such as neomycin resistance and the like.

[0124] Therefore, in another aspect, the invention relates to a vector, comprising the nucleotide sequence of the invention or the construct of the invention. The choice of the vector will depend on the host cell in which it is to be subsequently introduced. In a particular embodiment, the vector of the invention is an expression vector. Suitable host cells include a wide variety of prokaryotic and eukaryotic host cells. In specific embodiments, the expression vector is selected from the group consisting of a viral vector, a bacterial vector and a mammalian cell vector. Prokaryote- and/or eukaryote-vector based systems can be employed for use with the present invention to produce polynucleotides, or their cognate polypeptides. Many such systems are commercially and widely available.

[0125] Further, the expression vector may be provided to a cell in the form of a viral vector. Viral vector technology is well known in the art and is described, for example, in Sambrook et al. (2012), and in Ausubel et al. (1997), and in other virology and molecular biology manuals. Viruses, which are useful as vectors include, but are not limited to, retroviruses, adenoviruses, adeno-associated viruses, herpes viruses, and lentiviruses. In general, a suitable vector contains an origin of replication functional in at least one organism, a promoter sequence, convenient restriction endonuclease sites, and one or more selectable markers. (See, e.g., WO 01/96584;

WO 01/29058; and U.S. Pat. No. 6,326,193.

[0126] By way of illustration, the vector in which the nucleic acid sequence is introduced can be a plasmid which is or is not integrated in the genome of a host cell when it is introduced in the cell. Illustrative, non-limiting examples of vectors in which the nucleotide sequence of the invention or the gene construct of the invention can be inserted include a tet-on inducible vector for expression in eukaryote cells.

[0127] The vector may be obtained by conventional methods known by persons skilled in the art (Sambrook et al., 2012). In a particular embodiment, the vector is a vector useful for transforming animal cells. [0128] Tn some embodiments, the recombinant expression vectors may also contain nucleic acid molecules which encode a peptide or peptidomimetic activator of invention, described elsewhere herein.

[0129] A promoter may be one naturally associated with a gene or polynucleotide sequence, as may be obtained by isolating the 5' non-coding sequences located upstream of the coding segment and/or exon. Such a promoter can be referred to as “endogenous.” Similarly, an enhancer may be one naturally associated with a polynucleotide sequence, located either downstream or upstream of that sequence. Alternatively, certain advantages will be gained by positioning the coding polynucleotide segment under the control of a recombinant or heterologous promoter, which refers to a promoter that is not normally associated with a polynucleotide sequence in its natural environment. A recombinant or heterologous enhancer refers also to an enhancer not normally associated with a polynucleotide sequence in its natural environment. Such promoters or enhancers may include promoters or enhancers of other genes, and promoters or enhancers isolated from any other prokaryotic, viral, or eukaryotic cell, and promoters or enhancers not “naturally occurring,” i.e., containing different elements of different transcriptional regulatory regions, and/or mutations that alter expression. In addition to producing nucleic acid sequences of promoters and enhancers synthetically, sequences may be produced using recombinant cloning and/or nucleic acid amplification technology, including PCR™, in connection with the compositions disclosed herein (U.S. Patent 4,683,202, U.S. Patent 5,928,906). Furthermore, it is contemplated the control sequences that direct transcription and/or expression of sequences within non-nuclear organelles such as mitochondria, chloroplasts, and the like, can be employed as well.

[0130] Naturally, it will be important to employ a promoter and/or enhancer that effectively directs the expression of the DNA segment in the cell type, organelle, and organism chosen for expression. Those of skill in the art of molecular biology generally know how to use promoters, enhancers, and cell type combinations for protein expression, for example, see Sambrook et al. (2012). The promoters employed may be constitutive, tissue-specific, inducible, and/or useful under the appropriate conditions to direct high-level expression of the introduced DNA segment, such as is advantageous in the large-scale production of recombinant proteins and/or peptides. The promoter may be heterologous or endogenous. [0131 ] The recombinant expression vectors may also contain a selectable marker gene which facilitates the selection of transformed or transfected host cells. Suitable selectable marker genes are genes encoding proteins such as G418 and hygromycin which confer resistance to certain drugs, P-galactosidase, chloramphenicol acetyltransferase, firefly luciferase, or an immunoglobulin or portion thereof such as the Fc portion of an immunoglobulin, for example, IgG. The selectable markers may be introduced on a separate vector from the nucleic acid of interest.

[0132] Following the generation of the siRNA, miRNA, or anti-sense oligonucleotide, a skilled artisan will understand that the siRNA, miRNA, or anti-sense oligonucleotide will have certain characteristics that can be modified to improve the siRNA, miRNA, or anti-sense oligonucleotide as a therapeutic compound. Therefore, the siRNA, miRNA, or anti-sense oligonucleotide may be further designed to resist degradation by modifying it to include phosphorothioate, or other linkages, methyl phosphonate, sulfone, sulfate, ketyl, phosphorodithioate, phosphoramidate, phosphate esters, and the like (see, e.g., Agrwal et al., 1987, Tetrahedron Lett. 28:3539-3542; Stec et al., 1985 Tetrahedron Lett. 26:2191-2194; Moody et al., 1989 Nucleic Acids Res. 12:4769-4782; Eckstein, 1989 Trends Biol. Sci. 14:97-100; Stein, In: Oligodeoxynucleotides. Anti-sense Activators of Gene Expression, Cohen, ed., Macmillan Press, London, pp. 97-117 (1989)).

[0133] Any oligonucleotide may be further modified to increase its stability in vivo. Possible modifications include, but are not limited to, the addition of flanking sequences at the 5' and/or 3' ends; the use of phosphorothioate or 2' O-methyl rather than phosphodiester linkages in the backbone; and/or the inclusion of nontraditional bases such as inosine, queosine, and wybutosine and the like, as well as acetyl- methyl-, thio- and other modified forms of adenine, cytidine, guanine, thymine, and uridine.

[0134] In some embodiments of the invention, an anti-sense nucleic acid sequence which is expressed by a plasmid vector is used to inhibit ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2 protein expression. The anti-sense expressing vector is used to transfect a mammalian cell or the mammal itself, thereby causing reduced endogenous expression of ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2 protein.

[0135] Anti-sense molecules and their use for inhibiting gene expression are well known in the art (see, e.g., Cohen, 1989, In: Oligodeoxyribonucleotides, Antisense Activators of Gene Expression, CRC Press). Anti-sense nucleic acids are DNA or RNA molecules that are complementary, as that term is defined elsewhere herein, to at least a portion of a specific mRNA molecule (Weintraub, 1990, Scientific American 262:40). In the cell, anti-sense nucleic acids hybridize to the corresponding mRNA, forming a double- stranded molecule thereby inhibiting the translation of genes.

[0136] The use of anti-sense methods to inhibit the translation of genes is known in the art, and is described, for example, in Marcus- Sakura (1988, Anal. Biochem. 172:289). Such anti-sense molecules may be provided to the cell via genetic expression using DNA encoding the anti-sense molecule as taught by Inoue, 1993, U.S. Patent No. 5,190,931.

[0137] Alternatively, anti-sense molecules of the invention may be made synthetically and then provided to the cell. In some embodiments, anti-sense oligomers may have between about 10 to about 30 nucleotides. In some embodiments, anti-sense oligomers may have about 15 nucleotides. In some embodiments, anti-sense oligomers having 10-30 nucleotides are easily synthesized and introduced into a target cell. Synthetic anti-sense molecules contemplated by the invention include oligonucleotide derivatives known in the art which have improved biological activity compared to unmodified oligonucleotides (see U.S. Patent No. 5,023,243).

[0138] In some embodiments of the invention, a ribozyme is used to activate ARNT, ESRI, JUN, AR, ZEB1, SNAII, and/or SNAI2 protein expression by inhibiting a molecule that inhibits ARNT, ESRI, JUN, AR, ZEB1, SNAII, and/or SNAI2. Ribozymes useful for inhibiting the expression of a target molecule may be designed by incorporating target sequences into the basic ribozyme structure which are complementary, for example, to the mRNA sequence encoding an inhibitor of ARNT, ESRI, JUN, AR, ZEB1, SNAII, and/or SNAI2. Ribozymes may be synthesized using commercially available reagents (Applied Biosystems, Inc., Foster City, CA) or they may be genetically expressed from DNA encoding them.

[0139] In some embodiments, the activator of ARNT, ESRI, JUN, AR, ZEB1, SNAII, and/or SNAI2 may comprise one or more components of a CRISPR-Cas system, where a guide RNA (gRNA) targeted to a gene encoding ARNT, ESRI, JUN, AR, ZEB1, SNAII, and/or SNAI2, and a CRISPR-associated (Cas) peptide form a complex to repair a mutation within the targeted gene. In some embodiments, the activator comprises a gRNA or a nucleic acid molecule encoding a gRNA. In some embodiments, the activator comprises a Cas peptide or a nucleic acid molecule encoding a Cas peptide. [0140] Polypeptide activators

[0141] In other related aspects, the invention relates to modulators of one or more molecules involved in a cell transition (e.g. an isolated peptide activator that activates ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2). For example, in some embodiments, the peptide activator of the is a functional ARNT, ESRI, JUN, AR, ZEB1, SNAI1, and/or SNAI2 protein thereby increasing the level and/or functional activity of ARNT, ESRI, JUN, AR, ZEB1, SNAIl, and/or SNAI2.

[0142] The variants of the polypeptides according to the present invention may be (i) one in which one or more of the amino acid residues are substituted with a conserved or non-conserved amino acid residue and such substituted amino acid residue may or may not be one encoded by the genetic code, (ii) one in which there are one or more modified amino acid residues, e.g., residues that are modified by the attachment of substituent groups, (iii) one in which the polypeptide is an alternative splice variant of the polypeptide of the present invention, (iv) fragments of the polypeptides and/or (v) one in which the polypeptide is fused with another polypeptide, such as a leader or secretory sequence or a sequence which is employed for purification (for example, His-tag) or for detection (for example, Sv5 epitope tag). The fragments include polypeptides generated via proteolytic cleavage (including multi-site proteolysis) of an original sequence. Variants may be post-translationally, or chemically modified. Such variants are deemed to be within the scope of those skilled in the art from the teaching herein.

[0143] Antibody activators

[0144] The invention also relates to modulators of one or more molecules involved in a cell transition (e.g. an activator of ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 comprising an antibody, or antibody fragment, specific for ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2). That is, the antibody can activate ARNT, ESRI, JUN, AR, ZEB1, SNAIl, and/or SNAI2 or inhibit a molecule that inhibits ARNT, ESRI, JUN, AR, ZEB1, SNAU, and/or SNAI2 to provide a beneficial effect.

[0145] The antibodies may be intact monoclonal or polyclonal antibodies, and immunologically active fragments (e.g., a Fab or (Fab)2 fragment), an antibody heavy chain, an antibody light chain, humanized antibodies, a genetically engineered single chain Fy molecule (Ladner et al, U.S Pat. No. 4,946,778), or a chimeric antibody, for example, an antibody which contains the binding specificity of a murine antibody, but in which the remaining portions are of human origin. Antibodies including monoclonal and polyclonal antibodies, fragments and chimeras, may be prepared using methods known to those skilled in the art.

[0146] Antibodies can be prepared using intact polypeptides or fragments containing an immunizing antigen of interest. The polypeptide or oligopeptide used to immunize an animal may be obtained from the translation of RNA or synthesized chemically and can be conjugated to a carrier protein, if desired. Suitable carriers that may be chemically coupled to peptides include bovine serum albumin and thyroglobulin, keyhole limpet hemocyanin. The coupled polypeptide may then be used to immunize the animal (e.g., a mouse, a rat, or a rabbit).

EXPERIMENTAL EXAMPLES

[0147] The invention is further described in detail by reference to the following experimental examples. These examples are provided for purposes of illustration only and are not intended to be limiting unless otherwise specified. Thus, the invention should in no way be construed as being limited to the following examples, but rather, should be construed to encompass any and all variations which become evident as a result of the teaching provided herein.

[0148] Without further description, it is believed that one of ordinary skill in the art may, using the preceding description and the following illustrative examples, make and utilize the compounds of the present invention and practice the claimed methods. The following working examples therefore specifically point out exemplary embodiments of the present invention and are not to be construed as limiting in any way the remainder of the disclosure.

Example 1 : Learning transcriptional and regulatory dynamics driving cancer cell plasticity using neural ODE-based optimal transport

[0149] While single-cell technologies have allowed scientists to characterize cell states that emerge during cancer progression through temporal sampling, connecting these samples over time and inferring gene-gene relationships that promote cancer plasticity remains a challenge. To address these challenges, TrajectoryNet was developed, a neural ordinary differential equation network that learns continuous dynamics via interpolation of population flows between sampled timepoints. By running causality analysis on the output of TrajectoryNet, rich and complex genegene networks were computed that drive pathogenic trajectories forward. Applying this pipeline to scRNAseq data generated from in vitro models of breast cancer, a refined CD44^A'EPCAM⁺CAV1⁺ marker profile was identified and validated that improves the identification and isolation of cancer stem cells (CSCs) from bulk cell populations. Studying the cell plasticity trajectories emerging from this population, comprehensive temporal regulatory networks were identified that drive cell fate decisions between an epithelial-to-mesenchymal (EMT) trajectory, and a mesenchymal-to-epithelial (MET) trajectory. Through these studies, estrogen related receptor alpha was identified and validated as a critical mediator of CSC plasticity. TrajectoryNet was further applied to an in vivo xenograft model and demonstrates its ability to elucidate trajectories governing primary tumor metastasis to the lung, identifying a dominant EMT trajectory that includes elements of the disclosed newly defined temporal EMT regulatory network. Demonstrated here in cancer, the TrajectoryNet pipeline is a transformative approach to uncovering temporal molecular programs operating in dynamic cell systems from static single-cell data.

[0150] Cell state plasticity provides a mechanism for cancer cells to rapidly and dynamically evolve in a manner that facilitates primary tumor growth, metastasis and the development of therapy -resistant disease. While the advent of single-cell technologies has allowed detailed characterization of static cell states within tumors, further technological development is required to elucidate the mechanisms governing dynamic cell state transitions that may not only span days, months or years, but also create a variety of cell states shaping the tumor heterogeneity that drives disease progression. Furthermore, the transcriptional networks used by individual cells to undergo dynamic functional changes have been difficult to dissect due to the high dimensional nature of the data, as well as the computational challenge of resolving cellular trajectories over extended periods of time from static snapshot single-cell data. If these issues were addressed, it would be possible to use longitudinal patient samples to gain an unprecedented insight into the mechanisms governing metastasis and therapy-resistant disease, both of which have eluded scientists for decades.

[0151] Accordingly, a cell dynamics pipeline centered around TrajectoryNet was developed

[Tong, A., Huang, J., Wolf, G., van Dijk, D. & Krishnaswamy, S. Trajectory net: A dynamic optimal transport network for modeling cellular dynamics. Tn Proceedings of the 37th International Conference on Machine Learning (2020)], an algorithm for interpolating continuous dynamics from static snapshot data. Here, the ability of TrajectoryNet to learn cellular trajectories with enhanced capabilities of modeling cellular proliferation/death with an auxiliary proliferation network is not only improved, but tools to identify the gene networks underlying these trajectories are also added. TrajectoryNet is built on the new deep learning paradigm of neural ordinary differential equations (ODEs). In this paradigm, the neural network learns a time-varying derivative parameterized by the network weights and biases, and uses an ODE solver to integrate the derivative to compute the function at various points in time. Unlike the discrete dynamic models that are created by recurrent neural networks, neural ODEs can be trained to learn and interpolate dynamics continuously in time when given time course training data. Thus, these models combine the natural universality and trainability of neural networks with the capability of the ODE to model dynamics. When analyzing single cell data, however, one does not actually have cellular dynamics on any individual cell, since they are destroyed by the scRNAseq experiment; thus TrajectoryNet was trained to match populations over time using continuous normalizing flows. As continuous normalizing flows are not properly constrained to generate biologically plausible dynamics, regularizations were added that bias the network towards more energy-efficient and plausible paths. While this regularization helps TrajectoryNet match population distributions across time, in practice it actually gives continuously normalizing flows of each individual trajectory, thus recovering individual cellular trajectories, spanning long ranges in time and gaps in the cellular manifolds between samples. Furthermore, since cells divide and die, TrajectoryNet was equipped with an auxiliary network that learns proliferation rates of cells in order to reflect realistic cellular dynamics. As TrajectoryNet trains a model that learns cellular evolution over time, single cells can be projected into future time points, effectively imputing continuous gene expression profdes. Most critically, continuous single cell gene expression dynamics were combined with causality analysis and public gene regulatory databases to create trajectory-specific transcriptional networks that reveal the key transcription factors and genes that drive a trajectory forward. Through rigorous comparisons, it was shown that TrajectoryNet not only learns cellular trajectories better than other techniques, but also more accurately infers known gene regulatory relationships than established methods that do not incorporate cellular dynamics. CD44^/"EPCAM⁺CAV1 ⁺ [0152] The disclosed method was showcased on time-lapsed single cell data generated in in vitro and in vivo models to explore cancer cell state plasticity. First, TrajectoryNet was applied to a 3D in vitro tumorsphere model system [Chaffer, C. L. et al. Normal and neoplastic nonstem cells can spontaneously convert to a stem-like state. Proceedings of the National Academy of Sciences 108, 7950-7955 (2011)] to resolve temporal transcriptional dynamics that govern cancer stem cell (CSC) fate decisions. Using TrajectoryNet’ s computed proliferation rate, a new subset of CD44^A'EPCAM⁺CAV1⁺ CSCS was identified and enriched for tumor-initiating potential, demonstrating the power of TrajectoryNet to accurately identify specialized cells within a heterogeneous population. Tracking these cells over time, an epithelial-to-mesenchymal trajectory (EMT) was defined that drives CSCs towards a more mesenchymal cell state, and a mesenchymal -to-epithelial trajectory (MET) that drives CSCs into an epithelial non-CSC state. Using the disclosed causality framework, the first comprehensive temporally regulated transcriptional networks underlying each of these trajectories was built.

[0153] In cancer, the MET regulatory network is poorly understood. Hence, the ability of the TrajectoryNet pipeline to reveal novel gene regulatory networks by identifying estrogen receptor related alpha (ESRRA), an orphan nuclear receptor, as a key checkpoint in CSC fate determination and driver of the MET in cancer was showcased. It was shown that modulating ESRRA alters the transcriptional circuitry defining the epithelial and mesenchymal cell states and highlights its dynamic regulation throughout the in vitro tumorsphere-forming assay. It was further demonstrated that the temporally regulated gene regulatory networks identified by the TrajectoryNet pipeline accurately predict gene expression dynamics in scRNAseq data capturing primary tumor metastasis to the lung over a 16-week period in vivo. Here, TrajectoryNet reveals a dominant EMT trajectory in the in vivo metastasis data, with dynamic regulation of the newly defined core EMT regulatory network. As demonstrated herein, TrajectoryNet is an effective tool to model dynamic cell state transitions from snapshot single cell data and can identify genes critical in cell state plasticity.

[0154] First, a cell dynamics pipeline is presented featuring the TrajectoryNet Neural ODE network for inferring cellular — and associated gene — dynamics from single-cell data, and Granger causality analysis for building networks from these dynamics. Next, through rigorous comparisons, it was demonstrated that TrajectoryNet is not only able to infer trajectories significantly better than other methods [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928-943. e22 (2019); La Manno, G. et al. RNA velocity of single cells. Nature 560, 494- 498 (2018); Haghverdi, L , Btittner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016)], but can also more accurately infer gene regulatory relationships than methods that do not account for cellular dynamics [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in single-cell data. Science 346, 1250689-1250689 (2014)]. Next, introduced was a unique dataset consisting of a 5-timepoint in vitro tumorsphere assay that begins with CD44^hi cells containing a CSC population, measured with scRNA-seq, and tracks their evolution over 30 days. Next, the disclosed cell dynamics pipeline was applied to this dataset to derive temporal trajectories. From these trajectories 5 clusters of co-expressed genes were identified based on their temporal dynamics. The Traj ectoryNet-derived cellular trajectories allow us to derive a refined CSC marker set which includes EPC AM and CAV1 in addition to previously known CD44. 6) Next, two types of cellular trajectories were compared that emerge from these CSCs, one set of trajectories that is called MET trajectories as they have a more epithelial fate, while the other set, is called EMT trajectories as they have a more mesenchymal fate. 7) The disclosed causality tools were then applied to gene dynamics from these trajectories to derive the core gene regulatory circuitry of both the EMT and MET. 8) From these trajectories ESRRA was identified as a major driver of the MET transition, which could be responsible for creating the cell state heterogeneity required for secondary tumor formation. Experiments were conducted using siRNA and ESRRA antagonist treatment of 2D cultured CD44^hi cells, showing that they indeed lead to an increase in the epithelial marker CDH1 expression. 9) Finally, the relevance of the disclosed findings was tested in an in vivo system using a xenograft model where primary tumors were measured and matched lung metastases with single cell and spatial RNA-seq technologies. This experiment replicates the disclosed refined stem cell signature as well as elements of the EMT regulatory circuit.

Overview of the Cell Dynamics Pipeline

[0155] Time lapsed single-cell data (Figure 1A) poses significant challenges from a computational perspective. In cancer for example, cells undergo significant and rapid transcriptional changes in response to natural tumor evolution, the stressors of the metastatic cascade or therapeutic insults, where each individual timepoint represents a distribution of cells that largely do not overlap in cellular state with previous or subsequent timepoints. Under these circumstances, previously presented techniques for inferring dynamics such as RNA velocity [La Manno, G. et al. RNA velocity of single cells. Nature 560, 494-498 (2018)] and pseudotime [Haghverdi, L., Biittner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016)] fail as these approaches rely on well-sampled manifolds without large gaps in cell state space or time. Most critically, however, these approaches do not help users understand the gene regulatory dynamics necessary to drive these trajectories forward. Extracting longitudinal dynamics and understanding their underlying transcriptional drivers from disconnected static snapshot measurements can be challenging as there are few methods of interpolation between irregularly spaced or distant distributions and learn gene-gene regulatory dynamics. To address this knowledge-gap, TrajectoryNet was developed (Figure IB), a tool for studying cellular differentiation dynamics from time-lapsed single cell transcriptomic data. Notably, TrajectoryNet produces trajectories for each individual cell instead of a single averaged trajectory, as is the case with pseudotime, further allowing user to build sophisticated transcriptional networks by leveraging causality metrics and known gene-gene relationships. Here, a pipeline was presented that allows users to: i) learn single cell trajectories from time lapsed single cell data (Figure ICi), ii) interpolate continuous gene trends for each trajectory (Figure ICii) and iii) learn complex transcriptional networks from trajectories of interest (Figure ICiii).

[0156] TrajectoryNet learns the dynamics of cell state from cellular snapshot data (Figure 1A). Here its key features were highlighted and adaptations for use in this context and elaborated on herein. The heart of the model is a single neural network f with parameters 9 which takes as input a current cell state z(l) and the current time t and outputs the instantaneous change in cell

state with respect to time . This network /tells us how each cell evolves instantaneously and

continuously. To track its progress over longer time scales, standard ordinary differential equation (ODE) solvers can be used to integrate its position over time. Thus, the neural ODE system is an ensemble system consisting of a neural network computing a derivative and an ODE solver that integrates the derivative to compute a per cell time-trajectory. While the original neural ODE paper made the key contribution of showing how to train such a system, by using a method called adjoint sensitivity, here this framework was utilized to learn cellular population dynamics.

[0157] If longitudinal data was available — data that tracked a particular cell over time (Figure

1 A) — then a matching loss could be directly applied between the estimated cell at time

and the true cell position and apply gradient descent on However,

TrajectoryNet only assumes access to snapshot data, and so does not measure the cell z(to) at ti. Therefore, a distribution level loss was needed, for which the framework of normalizing flows was considered.

[0158] In a normalizing flow, a random variable zo is transformed via a bijective function zi = f(z) then the change in probability of is computed through a change of variables formula

[0159] In the original formulations of normalizing flows, discrete neural network layers would transform probability distributions from zo to zr, to Z2 and so forth, gradually transforming one distribution of cells to another via reversible operations. Neural ODEs allow a continuous version of this kind of normalizing flow using an instantaneous change of variables formula. If z(t) is a continuous random variable that is transformed via differential equation d

t then the change in probability is computed as This calculated

probability at each tk, where data was available is penalized by a KL divergence penalty for known q.

[0160] While continuous normalizing flows provide a way to transform one probability distribution to another, there is no reason that such flows need to be biologically plausible, as is needed to simulate the cellular differentiation process. Indeed, flows can be convoluted ways of altering one distribution via a series of circuitous paths. In order to constrain this to the realm of biological plausibility, a penalty was added that was believed to be true of biophysical systems: energy efficiency. In other words, it was believed that when cellular populations transform, they do so in an efficient way [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928-943. e22 (2019)]. This penalty is enforced by penalizing the magnitude of the derivative computed by the neural network at each measured timestep of the cellular differentiation process.

Interestingly, this regularization allows the neural network to perform a dynamic optimal transport, i.e., an optimal transport between the populations where the paths of transport (rather than just the final displacements) are energy optimal. This neural ODE design of dynamic optimal transport is indeed the key contribution of TrajectoryNet.

[0161] A penalty was also added to constrain the flow to the observable cellular manifold. This penalty restricts cellular paths to flow along regions of cellular density, as determined by the distance to the K-nearest neighboring cell. Finally, because cells are not simply transported but also proliferate and die, TrajectoryNet was altered to perform an unbalanced optimal transport as discussed herein.

[0162] The objective of TrajectoryNet can be decomposed into four parts:

[0163] The first term is a distribution matching term, which minimizes the divergence between the predicted distribution at time T and the true distribution at time T. This ensures that the model matches the data at each measured timepoint. [0164] The second term, energy integrates over the length of the path that the cell took from time 0 to time T to ensure that cellular paths are energy efficient. This is a key penalty in the disclosed framework that provably renders TrajectoryNet into a dynamic optimal transport framework. Of all the possible set of paths that link distributions together the one that minimizes the total path energy was chosen. Here, the energy was parameterized as the sum of squared path lengths:

[0165] Considering just these first two terms the objective is now equivalent to a relaxed optimal transport problem where Eq. 4 was minimized subject to the continuity equation dtp + 7- (pf) = 0 enforcing mass is neither created nor destroyed, and the boundary distributions match with p(x, 0) = p and a relaxed p(x, 1) ~ v where p(x, t) is the density at time t.

[0166] The third and fourth components are responsible for modeling cell proliferation and cell death help train a proliferation network g to transform the balanced OT problem in Eq. 4 into an unbalanced one. Intuitively, Lbaiance makes sure that the total mass of the system is not changing so that and Lproiiferation penalizes the (squared) cost of mass creation

and destruction enforcing how close to a

balanced optimal transport the solution is. These first four terms can be seen as solving a dynamic optimal transport problem:

[0167] which can be restated in the dynamic formulation for a vector flow field f and a scalar growth field g, as minimizing:

Equation 6

[0168] subject to the continuity equation which enforces that mass is

neither created or destroyed after taking into account a scalar addition pg, and subject to =

[0169] The final term ensures that the flow remains on the cellular manifold. While generally it was assumed that cells on average take short paths, cells are also thought to transition over a low dimensional manifold, an assumption which is also made in graph-based methods where paths are restricted to transitions between observed cells. In contrast, the flows are continuous in gene space and time, nevertheless it was still found useful to encourage paths that are near data. A thresholded k-nearest-neighbors penalty was used which regularizes the sum of distances to the k nearest neighbors:

[0170] Using these four terms yields an optimization similar to dynamic unbalanced optimal transport through direct optimization of the drift and proliferation fields. This direct optimization allows us to impose additional biological priors on the cell paths. Additional descriptions of these terms, the TrajectoiyNet algorithm and comparisons to other techniques is discussed herein. As disclosed below, rigorous benchmarking studies comparing TrajectoiyNet to other techniques on synthetic (Figure 7) and real single cell data (Figure 8) is described. .

Generating continuous gene expression data from trajectories

[0171] Once a trained mode was established as a dynamic optimal transport on single-cell

RNA-sequencing timecourse data, this model can be used to predict the path of a cell. Given a cell with state x at time the state forward and backwards in time can be integrated according to the ODE defined by to get the corresponding cell state at arbitrary

timepoints. Specifically, the cell state can be estimated according to the disclosed model

at any time

[0172] Furthermore, this equation can be numerically approximated using an ODE solver for each cell individually or parallelized over multiple cells. In practical terms, this allows users to infer the continuous gene expression dynamics for a given cell or cluster of cells based on its state x(ti) at its measured time ti.

Learning transcriptional networks with Granger causality

[0173] Continuous gene expression trends were combined with novel time lapsed causality analysis as well as public databases to learn the transcriptional networks underlying cellular trajectories. While many methods that infer gene-gene networks, like correlation or mutual information methods [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in single-cell data. Science 346, 1250689 (2014)], assume that the underlying cellular states are static, it is believed that incorporating cellular dynamics information is critical to inferring the true underlying transcriptional network. As TrajectoryNet outputs complete transcriptional trends for each cell across time, better gene-gene relationships can be inferred across time using Granger causality. Granger causality goes beyond correlations in time series by testing whether a variable (or set of variables) X forecasts a variable Y. More specifically, a variable X Granger causes a variable Y if predictions of Y based on its past values and on the past values of X are better than the predictions of Y based on only its own past values. To determine the likely gene regulatory structure, Granger causality analysis on the average TrajectoryNet trajectories was used for groups of cells. It was determined that there is a directed edge from a gene X to gene Y if gene X Granger-causes Y. This edge was scored based on the Granger p-value signed by the direction of influence of A on Y. This is determined by the sign of the coefficient on X in the linear regression of Y based on the past values of X and Y. If this coefficient is positive, this is taken to mean that Ahas a likely up-regulatory effect on Y. Conversely, if the linear regression coefficient of Ais negative, this is taken to mean Ahas a down-regulatory effect on Y (Figure 9A). By taking the log transformed Granger test statistic and adding a sign for direction of effect between a transcription factor and a target gene, an interaction strength score was computed that determines whether a transcription factor X may have an effect on a particular target gene Y. By aggregating signed Granger analysis results across large numbers of transcription factors and genes, transcriptional programs responsible for defining key trajectories can be identified.

[0174] It is noted that Granger causality is not the same as “true causality” as it is only evidence of X preceding Y, and there are many instances where this precedence may be incidental where changing X does not change Y. Therefore, these results are combined with known gene-gene relationships found in the Transcriptional Regulatory Relationships Unraveled by Sentencebased Text mining (TRRUST v2) database [Han, H. et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Research 46, D380-D386 (2017)]. While these databases are also imperfect, the integration of TrajectoryNet, causality analysis and a public gene regulatory network better approximates the underlying transcriptional network. In this instance, the key driver transcription factors were identified using a total Granger causal score (TGCS):

[0175] with the idea that key drivers will strongly Granger-cause many other genes. It was shown that TGCS outperforms other methods which assume static cellular states at inferring gene regulatory interactions from time series data (Figures 9B-D). Following TGCS scoring the TRUSST v2 database was then used to generate a seed network (using Cytoscape 3.9.0). This seed network was then used to extract the first-degree nodes of the top drivers of the trajectory (based on TGCS score) and their adjacent edges. The edges of this network were further pruned and weighed with Granger scores to create higher fidelity networks. It was however, cautioned that Granger scores do not prove causality and that these should be regarded as candidates for further validation. For a more in-depth explanation of the disclosed gene network inference comparisons on synthetic data see discussion herein.

Time-resolved single cell data for learning molecular drivers of cell fate decisions [0176] The TrajectoryNet pipeline was chosen to apply to study CSC plasticity, a major problem in cancer biology today. CSCs are highly specialized cells that drive tumor initiation and metastasis [Al-Hajj, M., Wicha, M. S., Benito-Hernandez, A., Morrison, S. J. & Clarke, M. F. Prospective identification of tumorigenic breast cancer cells. Proceedings of the National Academy of Sciences 100, 3983-3988 (2003)]. One of their key biological features is to use cell state plasticity to self-renew or differentiate into various non-CSC progeny to create tumor heterogeneity. To date, the temporal regulation and biological networks underlying CSC plasticity have not been resolved.

[0177] To investigate the cellular and transcriptional dynamics behind CSC plasticity and potentially identify new strategies to control it, the CD44^hi CSCs was suspended in a 3- dimensional tumorsphere assay. This population of heterogeneous CD44^/;' cells are enriched for CSCs that reside in a hybrid epithelial-mesenchymal state. In this assay, CD44^?" CSCs traverse at least two distinct trajectories: one where they self-renew and/or move towards a more mesenchymal CSC state via an EMT, and one where they differentiate into epithelial CD44^to non-CSC progeny via an MET [Chaffer, C. L. et al. Normal and neoplastic nonstem cells can spontaneously convert to a stem-like state. Proceedings of the National Academy of Sciences 108, 7950-7955 (2011)]. Note that both the trajectories were termed as EMT and MET started from the CSC state and went towards a more epithelial or more mesenchymal state (Figure 2A). To measure the range of cellular states created, tumorspheres were collected at days 2, 12, 18, and 30 and single cell RNA sequencing was performed (scRNA-seq) (Figure 2B). ScRNA-seq were also measured on the unsorted day 0 population. This dataset, comprising of 16,983 genes measured in 17,983 cells across 5 timepoints, creates a unique opportunity to study CSC transitions and better understand when and how cell fate decisions are made.

TrajectoryNet learns transcriptional dynamics driving cancer cell plasticity

[0178] TrajectoryNet was applied to the time-lapsed single cell tumorsphere measurements. By projecting TrajectoryNet inferred trajectories on top of a PHATE visualization [Moon, K. R. et al. Visualizing structure and transitions in high-dimensional biological data. Nature Biotechnology 37, 1482-1492 (2019)] of the time-lapsed single cell data, a diverse set of trajectories were observed starting at day 0 and ending at day 30 (Figure 2C). Visualizing TrajectoryNet’s inferred growth or proliferation rate, subsets of cells were identified that were disproportionately more likely to contribute to the following time point (Figure 2D). This inferred proliferation rate identifies rapidly dividing cells that contribute to the emergence of distinct phenotypic cell states during tumorsphere development. Next, an EMT signature score was computed based on the average expression of genes known to play a role in EMT [Chakraborty, P., George, J. T., Tripathi, S., Levine, H. & Jolly, M. K. Comparative study of transcriptomics-based scoring metrics for the epithelial-hybrid-mesenchymal spectrum. Frontiers in Bioengineering and Biotechnology 8 (2020)] and visualize this score on the disclosed PHATE embedding. Using this analysis, a heterogeneous range of epithelial non-CSCs (red) and hybrid epithelial-mesenchymal CSCs (blue) were clearly identified across timepoints (Figure 2E). When this score was taken into account and analyzed at the day 0 timepoint, which is the single cell data corresponding to the unsorted 2D cultured HCC38 cells, an epithelial cell fraction representing CD44^/o cells (red dots) and the hybrid epithelial-mesenchymal CD44^/!' fraction (blue dots) was observed. Critically, the trajectories defining cells that progress through to day 30 originate in cells within the CD44^hi fraction of the day 0 population, thus demonstrating the ability of TrajectoryNet to correctly identify tumorsphere cells-of-origin.

[0179] As TrajectoryNet learns the continuous gene expression dynamics, the complete expression dynamics of individual genes across the tumorsphere time-course can be interpolated. Visualizing these genes shows clear co-expression dynamics, allowing us to cluster the interpolated expression-vs-time trends of all trajectories into 5 co-expression clusters using k- means (Figures 2F-G). These gene co-expression clusters are ordered temporally, with Cluster 1 describing genes that have highest expression early in tumorsphere formation during day 0, and Cluster 4 describing genes that have highest expression late in tumorsphere formation during day 30. Finally, Cluster 5 describes a set of genes that have high expression both early and late in the transition (Figures 2F-G). A well-defined core EMT circuitry is critical to CSC plasticity [Yang, J. et al. Guidelines and definitions for research on epithelial-mesenchymal transition. Nature Reviews Molecular Cell Biology 21, 341-352 (2020)], yet little is known about the temporal regulation of these transcription factors in dynamic systems. To demonstrate that TrajectoryNet can resolve gene-specific dynamics, it was shown that the core EMT transcription factors were not simply co-expressed and co-regulated at the same time, rather, they exhibited distinct temporal expression patterns where SNAI1 peaked in Cluster 1, SNAI2 and TWIST 1 peaked in Cluster 3, while ZEB1 and ZEB2 in Cluster 5 have peak expression early and late across the tumorsphere time course. These data demonstrate that TrajectoryNet can delineate the complex temporal regulatory interactions that initiate and maintain cancer plasticity dynamics.

A refined CSC marker profde to improve identification of the tumorsphere cell-of-origin

[0180] Cell surface markers are used to enrich for breast CSCs from heterogeneous cell populations, e.g., combinations of CD44^Z", EPCAM^/o, CD104+ and CD24+ [14-16], yet only a subpopulation of those isolated cells actually contain tumor-initiating properties. The tumorsphere assay is an in vitro surrogate assay that measures tumor initiation capacity, where only CSCs are capable of seeding a tumorsphere when placed as single cells in suspension culture [Chaffer, C. L. et al. Normal and neoplastic nonstem cells can spontaneously convert to a stem-like state. Proceedings of the National Academy of Sciences 108, 7950-7955 (2011)]. To improve the identification of the tumorsphere cell-of-origin the output of TrajectoryNet’ s novel proliferation model was first used to identify the subpopulation of CD44^fe cancer stems cells that disproportionately initiate tumorspheres. At the day 0 time point comprising bulk 2D cultured cells CD44^;" and CD44^/o cell populations were observed as expected (Figure 3A). Using TrajectoryNet, a CD44*' cell population was identified with a high proliferation rate and computationally predicted their cell cycle phase [Tirosh, I. et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189-196 (2016).] (Figure 3B). This analysis reveals a significant overlap between TrajectoryNet’ s inferred proliferation rate and cells in the G2/M and S cell cycle stage, suggesting that those cells differentially contribute to tumorsphere development (Figure 3B). To test this hypothesis, CD44^Z" cells were isolated from each cell cycle phase by flow cytometry using the FUCCI cell cycle sensor system [Sakaue-Sawano, A. et al. Visualizing spatiotemporal dynamics of multicellular cell-cycle progression. Cell 132, 487-498 (2008)] and tested each population for tumorsphere-forming potential in vitro. As predicted, CD44^Z" cells that reside in the S and G2/M phase are significantly enriched for tumorsphere-forming potential compared to unsorted or G1 phase CD44^hl cells (Figure 3C).

[0181] TrajectoryNet was then used to identify new cell surface markers that improve and refine CSC isolation. Differential expression analysis was performed between CD44^hi and CD44^/o cells and identify surface markers highly expressed in the CD44^hi population whose expression overlaps with S/G2 phase cells including PTN, EPCAM, CAV1, MMP7, VCAN and ANAX5 (Figure 3D). Using Density Re-sampled Estimation of Mutual Information (DREMI) [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in single-cell data. Science 346, 1250689-1250689 (2014)] to measure non-linear associations between the expression of these markers and proliferation rate, it was confirmed that each gene associates either positively or negatively with proliferation (Figure 3D). This analysis shows that EPC AM (epithelial cell adhesion molecule) is, as predicted from previous literature, enriched in this population, and newly, that CAV1 (caveolin 1) positively associates with the TrajectoryNet predicted proliferation rate within the CD44'" day 0 subpopulation (Figure 3D).

[0182] To validate that EPCAM and CAV1 further enrich for tumorsphere-forming cells within the CD44*' fraction, single EPCAM⁺ cells, single CAV1⁺ cells, and double EPCAM⁺CAV1⁺ cells were isolated from the CD44^hi day 0 subpopulation using flow cytometry. After measuring tumorsphere-forming potential, it was found that neither EPCAM nor CAV1 alone further enriched for tumorsphere-forming potential, however, double EPCAM⁺CAV1⁺ cells formed significantly more tumorspheres (1.7 fold; p<0.05), while EPCAM⁺CAV1‘ cells showed significantly reduced (10 fold; p<0.001) than unsorted CD44^hi cells (Figure 3E). Thus, TrajectoryNet enabled the identification of biological properties (cell cycle) and a refined CD44^ft'EPCAM⁺CAVl⁺ CSC marker profile that improve the identification and isolation of CSCs from heterogeneous cell populations.

Tracing CSC plasticity along the EMT and MET trajectories

[0183] To define the transcriptional trajectories that drive CSCs through the EMT towards a mesenchymal cell state or the MET trajectory towards an epithelial cell state, the day 30 sample was zoomed into. Louvain clustering [Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment 2008, Pl 0008 (2008)] identified three cell clusters - one apoptotic cell cluster (defined by high expression of mitochondrial genes), one epithelial, and one mesenchymal cell cluster (defined by the EMT signature score) (Figure 4A) [Chakraborty, P., George, J. T., Tripathi, S., Levine, H. & Jolly, M. K. Comparative study of transcriptomics-based scoring metrics for the epi th elial -hybrid-mesenchymal spectrum. Frontiers in Bioengineering and Biotechnology 8 (2020)]. The existence of the epithelial and mesenchymal cell states were confirmed within the tumorspheres by immunofluorescent co-staining for ZEB1 (mesenchymal) or CDH1 (epithelial) protein expression at day 7 and day 28 (Figure 4B). This data shows fairly uniform ZEBl*'CDHl^/o expressing cells at day 7, consistent with expansion of the CSC pool at this early timepoint. In contrast, a striking dichotomy of cells emerging by day 28 was seen, where some cells retain the mesenchymal ZEBE'CDHE⁰ state, while others have transitioned to the ZEB UCDH l⁷" epithelial cell state. These data indicate that cell fate decision between the EMT and MET trajectories emerges after day 7.

[0184] To differentiate between the two cell fates, one more epithelial and one more mesenchymal, the cellular trajectories were traced starting from the two fates backward. One set of trajectories is called EMT trajectories, and the other MET trajectories, noting as before that they both start from the same hybrid CSC state (Figure 4C). Starting from CSCs, divergence of the EMT and MET trajectories was shown between day 12-30 (Figure 4D). To elucidate the gene regulatory network underlying these trajectories, the top 5,273 most highly and variably expressed genes were identified (genes that were above the fiftieth percentile in their dispersion and expression scores) as well as the top 461 most highly and variably expressed transcription factors, which were chosen by the same metrics. Causality analysis was performed using the disclosed signed Granger framework as described previously for every transcription factor-gene combination in the unpartitioned trajectory (Figure 4E), as well as the individual EMT and MET trajectories (Figure 10). Visualizing these signed Granger values across Combined, MET and EMT trajectories, enhancing regulatory relationships were clearly seen between transcription factors expressed in one gene cluster with genes expressed in the same and the following clusters. Furthermore, concomitant repressive regulatory relationships were identified with genes expressed in previous clusters. For instance, transcription factors in Cluster 2 have a strong enhancing regulatory relationships with genes in Cluster 2 and 3 while having repressive regulatory relationships with genes in Clusters 1 (Figures 4E).

[0185] Using the TGCS framework for each of the combined, EMT and MET trajectories, top 100 transcription factors that had the most regulatory effect within their respective trajectory were found. From these analyses the top 80 highly regulatory transcription factors shared between the combined, EMT and MET trajectories were found (Figure 4F), suggesting that these factors are required for survival within the tumorsphere regardless of cell fate. Thirty-seven (37) highly regulatory core transcription factors unique to the EMT trajectory and 23 highly regulatory core transcriptional factors unique to the MET trajectory were found (Figure 4F and Supplementary Figures 11A, 1 IB). For example, in the EMT trajectory, HES1 is down- regulated, while F0X03 and DDIT3 expression remain high (Figure 4G). In the MET trajectory, ARNT and ZEB1 remain low, while ESRRA levels follow an arc trajectory that increases early in the trajectory, then decreases at the end point when the epithelial state emerges (Figure 4G). Interestingly, 18 of 37 EMT-specific genes appear in Cluster 1, while the majority of MET- specific genes (11 of 23) appear later in Cluster 4 (Figure 4F). This data confirms the earlier observation that commitment to the MET trajectory occurs later in tumorsphere development (Figure 4B).

[0186] The EMT database was used (dbEMT 2.0 [Zhao, M., Liu, Y., Zheng, C. & Qu, H. dbEMT 2.0: An updated database for epithelial-mesenchymal transition genes with experimentally verified information and precalculated regulation information for cancer metastasis. Journal of Genetics and Genomics 46, 595-597 (2019)]) to determine if the core EMT and MET transcription factors identified by the TrajectoryNet pipeline have been implicated in the EMT. Interestingly, only 12 of the 37 EMT transcription factors (SKIL, F0X03, ELK3, SOX9, NCOA3, DNMTI, ETS2, SNAH, RUNX3, HES1, SNAJ2, TBX2\ and only 4 of the 23 MET transcription factors (ESRRA, ETV1, TRPS1, ZEBE), are in the database. These data show that the TrajectoryNet pipeline has uncovered a majority of new transcriptions yet to be associated with these cancer cell plasticity programs.

[0187] Gene regulatory networks were therefore built around the EMT and MET core transcription factors using the TRRUSTv2 database (Figure 4H and Supplementary Figure 12). Of note, the EMT core hub (green) regulates a unique set of downstream EMT interactors (purple), while the MET core hub (orange) regulates a different set of unique downstream MET interactors (yellow). There is also a common set of downstream interactors (grey node) between the EMT and MET core hubs, suggesting that those genes lie at the interface of the two transcriptional circuits and may therefore play a causal role in cell-fate decisions between the two trajectories. Interestingly, there are minimal direct regulatory interactions between the core EMT and MET transcription factors themselves. Notable exceptions are the MET transcription factor ATF3 regulation of DDIT3, and the MET transcription factor ESRRA regulation of the SNAI1 and SNAI2 transcription factors (Figure 4H).

[0188] Analysis of the top significantly enriched pathways (p < 0.05) within the MET subcluster (comprising 98 genes made up of the core MET transcription factors, unique downstream MET interactors, and, common interactors) were API pathway, AHR pathway, and regulation of miRNA metabolic process. Top pathways (p < 0.05) enriched in the EMT subcluster (comprising 251 genes made up of the core EMT transcription factors, downstream EMT interactors, and, common interactors) were hemopoiesis, regulation of cell death and transcriptional misregulation in cancer pathways (Figure 4H and Supplementary Figure 12). These analyses resolve, for the first time, the distinct gene regulatory networks and signaling pathways required to drive CSC plasticity along the EMT versus MET trajectory.

Validating the temporal MET gene regulatory network

[0189] The MET is critical for driving CSC differentiation into epithelial non-CSC progeny to create the heterogeneity required for robust tumor growth [Castano, Z. et al. Il- 1 P inflammatory response driven by primary breast cancer prevents metastasis-initiating cell colonization. Nature cell biology 20, 1084 (2018)]. Yet to date, there is no comprehensive regulatory network that describes the MET gene regulatory network. TrajectoryNet and the Granger Causality analysis was therefore used to build and validate one based on the disclosed data. Briefly, the TRRUST v2 database was used to extract and visualize the gene regulatory relationships of the 23 core transcriptional factors regulating the MET trajectory (Figure 5 A) [Shannon, P. et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Research 13, 2498-2504 (2003)]. In this network, the core transcription factors are annotated by solid rectangles (Figure 5B). Using the EMT database (dbEMT 2.0 [Zhao, M., Liu, Y., Zheng, C. & Qu, H. dbEMT 2.0: An updated database for epithelial-mesenchymal transition genes with experimentally verified information and precalculated regulation information for cancer metastasis. Journal of Genetics and Genomics 46, 595-597 (2019)]), it was shown that the core MET transcription factors not only regulate known EMT genes (circles), but also novel genes yet to be associated with the EMT (diamonds) (Figure 5B). [0190] By color-coding and organizing the core transcription factors and their associated gene regulatory networks by Gene Cluster number, the temporal cross-talk was visualized between early time point Clusters 2 (green) and 3 (blue) and late time point Clusters 4 (orange) and 5 (red) (Figure 5B). Notably, ATF3 is involved at the earliest time point of the MET initiation (Gene Cluster 2) followed by ESRRA, ETV1, NFATC2 (Gene Cluster 3), NFAT5, ASH1L (Gene Cluster 4), and then ZEB1, ARNT, TRSPI and ZNF350 (Gene Cluster 5). It was also shown that the Gene Cluster nodes commonly interact with four factors (JUN, AHR, AR, ESRF) (white circles), suggesting that these genes may serve as intermediary signaling nodes between the gene clusters (Figure 5B).

[0191] Uniquely, ESRRA and PAX9 (Gene Cluster 3) expression peak midway through the MET trajectory, then decrease as the epithelial cell state emerges, whereas all other core MET transcription factors are down-regulated (Figure 11 A). Due to poorly defined regulatory relationships of PAX9 in public data, its direct regulatory impact on the MET network was unable to be determined. Clear interactions for ESRRA, however, were observed. Thus, the fact that ESRRA expression peaks early in the MET trajectory (Cluster 3), and that ESRRA directly suppresses SNA11 and SNA12, suggests that ESRRA may have two critical functions in the MET trajectory: (i) to initiate the MET transcriptional program through its direct downstream interactors in the MET circuitry, and (ii) directly suppress the core EMT transcriptional circuitry via down -regulation of SNAH and SNAJ2. Accordingly, the regulatory role of ESRRA as a putative instigator of the MET trajectory was further delved into (Figure 5C).

[0192] While initially regarded as an orphan nuclear receptor, recent work has shown that ESRRA is closely related to the estrogen receptor and is commonly associated with poor outcome and poor prognosis in triple-negative breast cancer [Berman, A. Y. et al. ERRci regulates the growth of triple-negative breast cancer cells via s6kl -dependent mechanism. Signal Transduction and Targeted Therapy 2 (2017)]. Reducing ESRRA expression via antagonists or gene knockdown decreases cell proliferation and tumorigenicity [Berman, A. Y. et al. ERRa regulates the growth of triple-negative breast cancer cells via s6kl -dependent mechanism. Signal Transduction and Targeted Therapy 2 (2017); Ma, J.-H. et al. STAT3 targets ERR-a to promote epithelial-mesenchymal transition, migration, and invasion in triple-negative breast cancer cells. Mol. Cancer Res. 17, 2184-2195 (2019)]. Importantly, ESRRA directly alters the expression of EMT transcription factors SNAI1/2XO regulate the EMT in triple-negative breast cancers [De Luca, A. et al. Mitochondrial biogenesis is required for the anchorage-independent survival and propagation of stem-like cancer cells. Oncotarget 6, 14777-14795 (2015); Wu, Y.-M. et al. Inhibition of ERRa suppresses epithelial mesenchymal transition of triple negative breast cancer cells by directly targeting fibronectin. Oncotarget 6, 25588-25601 (2015); Huang, J.-W. et al. Effects of estrogen-related receptor alpha (ERRa) on proliferation and metastasis of human lung cancer A549 cells. J. Huazhong Univ. Sci. Technolog. Med. Sci. 34, 875-881 (2014); Chen, S., Ye, J., Kijima, I., Kinoshita, Y. & Zhou, D. Positive and negative transcriptional regulation of aromatase expression in human breast cancer tissue. J. Steroid Biochem. Mol. Biol. 95, 17-23 (2005)]. Accordingly, based on published interactions and gene-expression trends within the MET network, a module was identified (Figure 5C) wherein ESRRA directly (SNAI1, SNAI2) and indirectly (ZEB1') suppresses core EMT transcription factors to enable the upregulation of CDH1 and emergence of the epithelial cell state, whilst also modulating AHR- and AR- dependent pathways that may be critical for the emergence of the epithelial cell state. Indeed it was shown that AR is critical for induction of the EMT that creates new CSC populations in response to chemotherapeutic insults [San Juan, B. P. et al. Targeting phenotypic plasticity prevents metastasis and the development of chemotherapy-resistant disease. medRxiv (2022)]. Furthermore, ESRRA interacts with several core intermediary nodes, JUN, AR, ESRI, and with core transcription factors in Cluster 5 (AHR, ARNT) (Figure 5C). Together, these data suggest the ESRRA sub-network is a regulatory circuitry traversing temporal gene clusters to drive the MET trajectory and fortify the emergence of the epithelial cell state.

[0193] To test this hypothesis, ESRRA protein expression across the 3D tumorsphere timecourse by immunofluorescence (IF) was analyzed. Protein levels of ZEB 1, to mark the mesenchymal cell state, and CDHI, to mark epithelial cells was also analyzed (Figure 5D). These studies confirm critical temporal regulation of ESRRA during the tumorsphere timecourse. First it was noted that ESRRA is co-expressed with ZEB1 at day 7 demonstrating that it is present in the CSC population. However, a marked increase in ESRRA expression from day 7 to day 14 was observed. Interestingly, ZEB1 down-regulation follows peak ESRRA expression. ESRRA then decreases for the remainder of the time-course, such that at day 28 ESRRA and ZEB1 are down-regulated, while CDHI is up-regulated marking the presence of epithelial cell state. These protein expression findings are confirmed when comparing TrajectoryNet inferred expression dynamics with the protein (IF) trends (Figure 5E). Together, these data confirm unique temporal gene expression patterns not previously identified, whereby ESRRA expression is required for the CSC state, yet an increase in ESRRA expression triggers a down-regulation of EMT transcription factors SNAI1, SNAI2, and ZEB1, thereafter, ESRRA itself is downregulated to enable CDH1 upregulation and the emergence of the epithelial cell state. To confirm a causal relationship of ESRRA expression with the emergence of the epithelial cell state as observed at day 28, it was demonstrated that siRNA and or ESRRA antagonist treatment of 2D cultured CD44^Z" cells lead to a significant increase in CDH1 expression (p < 0.05) (Figure 5F-G and Figure 13).

[0194] Thus, TrajectoryNet enabled identification of a core regulatory unit defining the first comprehensive temporal MET regulatory network in breast cancer with ESRRA confirmed as a CSC marker and early key regulator of cell fate decisions towards the epithelial cell state.

Validating TrajectoryNet gene expression dynamics from primary tumors to lung metastases

[0195] In this disclosure, the processes of EMT and MET were followed using in vitro 5- timepoint data sampling. While this allowed following the dynamics of cancer cell plasticity carefully in time, to establish the relevance of these dynamics to in vivo systems was also desired. To do this, a xenograft model of metastatic triple-negative breast cancer was used. ScRNAseq and spatial transcriptomics was performed on primary tumors and matched lung metastases from four mice bearing the MDA-MB-231 xenograft. Following orthotopic implantation in the mammary fat pad, primary tumors were resected at 12 weeks, lung metastases at 16 weeks, and both samples were analyzed thereafter by scRNAseq (Chromium 10X) and spatial transcriptomics (Visium 10X, primary tumors only) (Figure 6A). The scRNAseq data was embedded into 3 dimensions using PCA and applied TrajectoryNet to map two trajectories. Firstly, the differentiation of CSCs within a primary tumor were mapped, and secondly, the trajectories of CSCs originating in the primary tumor that migrate to the lung to form a metastasis were mapped (Figure 6B).

[0196] To validate the refined CD44^A'EPCAM⁺CAV1⁺ CSC marker profile, regions of the primary tumors were first visualized and lung metastases expressing high levels of established CSC markers CD44, EPCAM and ZEB1 (dotted circles) (Figure 6C). Supporting the validation of CAV1 as a new cell surface marker to refine identification of the CSC state, it was confirmed that CAV1 is positively correlated in both samples with CSC markers CD44 (primary p=0.19, lung p=0.40), EPCAM (primary p=0.86, lung p=0.97), and ZEB1 (primary p=0.92, lung p=0.48) (Figure 6C). To further validate these findings, the primary tumor scRNAseq data was integrated with matched spatial transcriptomic data using scMMGAN [Amodio, M. et al. Single-cell multimodal GAN reveals spatial patterns in single-cell data from triple-negative breast cancer. Patterns (N Y) 3, 100577 (2022)], and show their spatial overlap in primary tumors (Figure 6D). The co-localization of these CSC markers supports the identification of regions of primary and lung metastases enriched for cells residing in the refined CSC state (Figure 6E) that was identified by the prior cell-of-origin analysis (Figure 3D).

[0197] Next, the EMT score was applied [Chakraborty, P., George, J. T., Tripathi, S., Levine, H. & Jolly, M. K. Comparative study of transcriptomics-based scoring metrics for the epithelial- hybrid-mesenchymal spectrum. Frontiers in Bioengineering and Biotechnology 8 (2020)] to identify regions of the tumors enriched for the epithelial or mesenchymal cell state (Figure 6F). This analysis shows that CSCs reside in areas of the tumor enriched for epithelial markers, reflecting their hybrid epithelial/mesenchymal phenotype. Radiating out from the CSC-enriched area, a predominant emergence of mesenchymal cells was observed, suggesting that CSCs in this in vivo model predominantly follow an EMT trajectory. This is not surprising given that MDA- MB-231 cells retain their CD44^hi state in vivo, which is associated with a more mesenchymal phenotype [Bierie, B. et al. Integrin-P4 identifies cancer stem cell-enriched populations of partially mesenchymal carcinoma cells. Proceedings of the National Academy of Sciences 114, E2337-E2346 (2017)]. Unbiased clustering was then performed to identify the region of the primary tumors and lung metastases associated with the most mesenchymal cell state (green) to highlight the emergent EMT trajectory (Figure 6G).

[0198] Finally, the dynamics within the in vivo EMT trajectories of genes associated with components of the newly-defined core EMT hub from Figure 4H were examined, including TBX2, DDIT3, MXI1, DNMT1, ETS2 and FLYWCH1. Complementary temporal gene expression patterns associated with the EMT trajectory within the primary tumor were observed (and highlighted their unique spatial distributions within the primary tumor; Figure 14), and within the primary to lung metastasis trajectory (Figure 6H). For example, TBX2 is linearly up-regulated in the primary tumor trajectory, whereas, in the primary to lung metastasis trajectory, its expression is initially flat at the beginning, reflecting the dynamics of the primary CSC to lung CSC portion of the trajectory, and thereafter linearly increases reflecting the lung CSC to lung mesenchymal portion of the trajectory. These dynamics are also reflected in the continuous expression of DDIT3, MXI1, DNMT1 and ETS2 (Figure 6H). These analyses demonstrate that the refined EMT core circuitry, that was identified in an in vitro model system, is conserved across distinct triplenegative breast cancer cell lines (HCC38 versus MDA-MB-231), from in vitro to in vivo systems, and across tumors growing in vastly different tissues (mammary fat pad versus lung tissue).

[0199] Together these analyses show that the TrajectoryNet pipeline identifies critical markers enabling improved classification and identification of the CSC state, and, defines biological trajectories that define distinct CSC trajectories and ultimately, cell fates, across disparate timepoints and organs in in vitro and in vivo model systems. Accordingly, the ability to characterize and quantify a tumor’s CSC component as a measure of tumor aggressiveness, coupled with the discovery of temporal gene regulatory strategies to drive CSCs towards different cell fates, is a critical advance for the discovery of temporal gene regulatory networks underlying these critical biological processes.

Discussion

[0200] Dynamic cell state evolution is a hallmark of many diseases. In cancer, dynamic transitions manifest as changes from relatively benign cell states to highly aggressive, proliferative and therapy-resistant cell states. Until now, it has been difficult to study these dynamics comprehensively because single-cell and other high throughput data modalities have largely captured static snapshot measurements. While learning dynamics has often been a holy grail for this type of data, previous efforts have significant shortcomings. Pseudotime [Haghverdi, L., Biittner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016)] and RNA velocity [La Manno, G. et al. RNA velocity of single cells. Nature 560, 494-498 (2018)] methods are useful in certain contexts, however they are not able to compute trajectories at the single cell level, interpolate between distant or disconnected distributions of cells, or learn transcriptional networks underlying cellular state transitions. To address the knowledge-gap in learning data dynamics, the TrajectoryNet pipeline was developed, which implements a breakthrough neural network paradigm called neural ODE to learn a dynamic optimal transport between time-lapsed single cell populations. The gene dynamics information from TrajectoryNet was then used to perform time-lagged Granger causality analysis to learn dynamic transcriptional networks that drive trajectories forward. It was demonstrated that TrajectoryNet is not only able to infer trajectories significantly better than other OT [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928-943. e22 (2019)], splicing [La Manno, G. et al. RNA velocity of single cells. Nature 560, 494-498 (2018)] and graph based [Haghverdi, L., Biittner, M., Wolf, F. A., Buettner, F. & Theis, F. I. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016)] methods, but also can more accurately infer gene regulatory relationships than methods that do not account for cellular dynamics, like correlation and mutual information [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in single-cell data. Science 346, 1250689-1250689 (2014)] based techniques.

[0201] Here, TrajectoryNet was applied to identify the cellular initiators and molecular drivers of CSC plasticity using a triple-negative breast cancer model. First, it was shown that TrajectoryNet enables improved characterization and isolation of cells residing in the CSC state. Using the tumorsphere assay, an in vitro surrogate of tumor-initiating potential, cells at the end of the trajectory were traced back to their cell-of-origin within the starting population of cells, thereby identifying the proliferative subpopulation of cells most likely to represent the CSC population that initiates tumorspheres. Applying this strategy, TrajectoryNet reveals a refined CD44*'EPCAM⁺CAV1⁺ CSC marker profile that was validated to identify cells significantly enriched for tumorsphere-forming potential above the CD44^luEPCAM⁺ markers previously described [Al-Hajj, M., Wicha, M. S., Benito-Hernandez, A., Morrison, S. I. & Clarke, M. F. Prospective identification of tumorigenic breast cancer cells. Proceedings of the National Academy of Sciences 100, 3983-3988 (2003)]. This strategy can now be applied in other settings, for example, to define cell populations that drive site-specific metastasis, or to identify cells that emerge in therapy-resistant cell states following cytotoxic treatments. These types of analyses will have profound impact on the future development of therapeutic strategies to specifically identify and target aggressive subpopulations of cancer cells. [0202] Importantly, the ability of the TrajectoryNet pipeline to define detailed temporal analysis of transcription factors and their associated targets regulating dynamic CSC plasticity programs in 3D in vitro tumorsphere models and in in vivo models of breast cancer metastasis was demonstrated. Accordingly, the first comprehensive gene regulatory networks defining the EMT trajectory that drive CSCs towards a mesenchymal cell fate, and the MET trajectory that drives CSC differentiation into epithelial non-CSC progeny were identified. Given that the details of the MET program were largely unknown, an MET subnetwork for further validation was identified, with ESRRA as a potential upstream initiator. ESRRA is present in CSCs was shown, and that its expression peaks to initiate the MET then regresses to enable the emergence of the non-CSC state. Indeed, the ESRRA regulatory network comprises layered and complex gene regulation acting to both initiate an epithelial program was shown, potentially through regulation of AHR:ARNT signaling pathway [Mulero-Navarro, S. & Fernandez-Salguero, P. M. New trends in aryl hydrocarbon receptor biology. Frontiers in Cell and Developmental Biology 4 (2016)], whilst simultaneously suppressing the mesenchymal program through direct and indirect regulation of the core EMT circuitry (SNAI1, SNAI2, ZEB1 and CDH I ). Consistent with the disclosed findings, previous studies have shown that ESRRA silencing results in the upregulation of CDH1 in MDA-MB-231 cell lines (TNBC model) and HEC-1A cells (human endometrial adenocarcinoma model) [Yoriki, K. et al. Estrogen-related receptor alpha induces epithelial-mesenchymal transition through cancer-stromal interactions in endometrial cancer. Scientific Reports 9 (2019)].

[0203] Thus, the identification of the ESRRA subnetwork as a strategy to regulate CSC plasticity could lead to the validation of novel targets that prevent metastatic outgrowth by driving CSC differentiation into a non-CSC state that is sensitive to chemotherapy. Future studies validating the role of the ESRRA network in metastasis and chemotherapy-resistance will shed light on the importance of these findings in the clinical setting [Berman, A. Y. et al. ERRa regulates the growth of triple-negative breast cancer cells via S6K1 -dependent mechanism. Signal Transduct. Target. Ther. 2, 17035 (2017); Carey, L. A. et al. The triple negative paradox: primary tumor chemosensitivity of breast cancer subtypes. Clin. Cancer Res. 13, 2329-2334 (2007); Bray, F. et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 68, 394-424 (2018)]. [0204] Finally, the ability of the TrajectoryNet pipeline to reveal gene regulatory networks in complex cellular transitional models was highlighted, demonstrated here spanning a 16 week in vivo time course. Temporal regulation of a subset of newly-defined regulatory components of the EMT was shown, highlighting complementary gene expression dynamics that hold true across the in vitro and in vivo models. The TrajectoryNet pipeline is a breakthrough method that can unravel complex temporal genomic networks governing cell state dynamics. Here the ability of TrajectoryNet to reveal novel dynamic cancer gene networks was demonstrated, however, the disclosed pipeline can be used to study any type of dynamic transition captured via single cell technology over multiple timepoints including stem cell differentiation, response to therapeutic interventions, or infections. As longitudinal data becomes increasingly available and critical in the study of dynamic systems, TrajectoryNet will be applied to a broader set of datasets and types.

Overview of TrajectoryNet

[0205] TrajectoryNet [Tong, A., Huang, J., Wolf, G., van Dijk, D. & Krishnaswamy, S. Trajectory net: A dynamic optimal transport network for modeling cellular dynamics. In Proceedings of the 37th International Conference on Machine Learning (2020)] learns a continuous normalizing flow between cross-sectional snapshot measurements in such a way that the flow represents biologically plausible paths of development and differentiation. The key to obtaining paths that adhere to previous knowledge of cell development was to penalize a continuous normalizing flow [Grathwohl, W., Chen, R. T. Q., Bettencourt, J., Sutskever, I. & Duvenaud, D. FFJORD: Free-form Continuous Dynamics for Scalable Reversible Generative Models. In ICLR (2019). 1810.01367], implemented as a neural ODE, such that it performs dynamic optimal transport [Benamou, I.-D. & Brenier, Y. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numerische Mathematik 84, 375-393 (2000)]. This enables us to efficiently solve the problem of cells evolving over time as a transport problem which is penalized over path lengths of the cells (See Eq. 3).

Neural Ordinary Differential Equation Model

[0206] TrajectoryNet is built on a neural ordinary differential equation (neural ODE) framework

[Chen, R. T. Q., Rubanova, Y ., Bettencourt, J. & Duvenaud, D. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31 (2018). 1806.07366], As previously defined, a neural ODE defines a neural network /with parameters 0 which defines a time derivative for the flow at every point and time. Thus for any time t the flow is defined by ome initial time to, for any x(t), it was that

. This integral could be computed using Euler integration where at a

fixed set of times was calculated that

incrementally. However, this can accumulate significant error depending on the dynamics, therefore ODE solvers were used with lower error and or adaptive step sizes. An important question is how to compute the loss with respect to the parameters 0. For this, the adjoint method was used. Given a loss function L

( ) then, the gradient was taken with respect to the parameters 9 using the adjoint dynamics By running the adjoint dynamics,

was computed without having to store the intermediate steps of the forward-time ODE

solver, which greatly reduces memory usage.

[0207] An important note here is that F is approximately invertible due to the deterministic dynamics. However, the time derivative determined by the network / need not be invertible, and is representable by any standard neural network architecture. This is in contrast with the architectures used in invertible neural networks that directly represent an invertible function, and must be severely restricted for fast Jacobian computation [Rezende, D. J. & Mohamed, S. Variational Inference with Normalizing Flows. In Proceedings of the 32nd International Conference on Machine Learning, vol. 37, 1530-1538 (2015). 1505.05770],

Continuous Normalizing Flows

[0208] While TrajectoryNet models the flow of cells over time, it is trained from cross-sectional population data. At a set of timepoints there was a corresponding set of

discrete samples X from the population of viable cells at that time.

Importantly, as with other destructive measurement techniques, a cell present at time t, is destroyed and is cannot be measured at time ti+i. Instead, the most likely dynamics of single cells from the population that is present at each timepoint must be inferred. To accomplish this, the optimization as a continuous normalizing flow (CNF) was formulated, which transforms an initial known at time to to match an empirical distribution at time ti. Since

approximately invertible by definition, and one can efficiently calculate the divergence of the probability distribution then given a probability distribution P(to) one can find P(ti) = Fn(p(to)) where F denotes the pushforward operator of for any subset

A continuous normalizing flow also defines a likelihood at every point x_t.. One can calculate log for any x_ti using the change of variables, and in particular the instantaneous

change of variables as derived in [Chen, R. T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31 (2018). 1806.07366],

[0209] where is the density of a standard Gaussian distribution. This log

likelihood definition allows the CNF to be trained using maximum likelihood, i.e. trained using the loss

Continuous Normalizing Flows with Multiple Target Timepoints

[0210] Previously continuous normalizing flows were applied to a single timepoint with a known distribution at time tO and a single empirical distribution at time tl. This can be used for generative modeling of complex distributions such as distributions of images [Chen, R. T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31 (2018). 1806.07366] or medical charts [Rubanova, Y ., Chen, R. T. Q. & Duvenaud, D. Latent ODEs for Irregularly-Sampled Time Series. arXiv: 1907.03907 [cs, stat] (2019). 1907.03907], Here, to model a series of distributions over time with one smooth time varying flow fe was desired. Tn this work, a few tricks to train this longer flow over time were applied. One way to accomplish this would be to integrate points sampled from and so on backwards to and at each time step apply a

distribution matching loss between empirical distributions

Such a loss is used in [Yang, K. D. & Uhler, C. Scalable Unbalanced Optimal Transport Using Generative Adversarial Networks. In 7th International Conference on Learning Representations, 20 (2019); Hashimoto, T. B., Gifford, D. K. & Jaakkola, T. S. Learning Population-Level Diffusions with Generative Recurrent Networks. In Proceedings of the 33rd International Conference on Machine Learning, 2417-2426 (2016)]. However, distribution type losses are difficult to calculate and to optimize, and while there is significant ongoing work improving the stability, efficiency, and flexibility of such losses, there are still significant drawbacks.

[0211] Instead, this was tackled by integrating samples from all timepoints in parallel backwards to a timepoint to where one can define a continuous probability density function With a

single integration,

was calculated, and therefor maximized.

Optimal Transport

[0212] Informally, optimal transport considers the problem of moving one pile of mass to another at minimum cost. Optimal transport considers a ground distance defined between points and lifts this to distances between measures. Define the metric measure space

then the squared 2-Wasserstein distance between two measures defined on

[0213] where II(p, v) is the set of joint distributions with marginals p and v respectively and The optimal transport distance can be thought of the minimal total cost of

moving mass piled at p to mass piled at v. This optimization is in general difficult to compute for measures over a continuous space such as . It is shown how to approximately solve a

dynamic version of the optimization with TrajectoryNet in the following sections.

Unbalanced Optimal Transport [0214] A common variant of optimal transport is unbalanced optimal transport where the mass equality constraint is relaxed. Intuitively, there may be cases where for some additional cost, in addition to transporting mass, one can “teleport” mass for some additional cost. The unbalanced transport problem balances between this transport and teleportation cost for each mass. The unbalanced optimal transport problem is defined with respect to some ^-divergence

often taken to be the Kullback-Leibler divergence (DKL), but in principle can be used with more general cp divergences. Given

and a regularization parameters the unbalanced optimal

transport problem was defined as

[0215] where

controls the relative cost of transportation vs. teleportation. To model cell division and cell death, unbalanced transport was used.

Dynamic Optimal Transport

[0216] Dynamic optimal transport adds a time component to the static version of optimal transport. Adding a time component is useful as it gives a way to link the theory of optimal transport to that of dynamical systems and in particular fluid dynamics [Benamou, J.-D. & Brenier, Y. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numerische Mathematik 84, 375-393 (2000)]. Instead of optimizing over the matching between distributions, the dynamic formulation considers an optimization over time- parameterized paths. For a given interval [to, ti\ with a source distribution u. at time to and a target distribution v at time ti, a time-dependent probability distribution pt(x) and a time dependent vector field f(x, t) was defined, such that if the probability distribution evolves according to the continuity equation

[0219] In other words, a velocity field f(x, t) with minimum L² norm that transports mass at p to mass at v when integrated over the time interval is the optimal plan for an L² Wasserstein distance. This can be seen by noticing that the optimal paths for each point pair (xo, xi) are geodesics and that inf/ This problem is in general challenging to

solve, particularly in the high dimensional case where a discretization of space-time (such as used in [Papadakis, N., Peyre, G. & Oudet, E. Optimal Transport with Proximal Splitting. SIAM Journal on Imaging Sciences 7, 212-238 (2014). 1304.5784] is impractical. CNFs were used to parameterize and learn the vector field f which approximates dynamic optimal transport.

Dynamic Optimal Transport with a Continuous Normalizing Flow

[0220] Continuous normalizing flows (CNFs) traditionally model an unknown data density by mapping it to a normal distribution using a continuous-time invertible function. While these are traditionally used only to “normalize” a single distribution, there is nothing in principle stopping one from mapping from L — > to — > tnorm where one has a known normal density at tnorm. By adding these additional steps and an energy regularization, it can be shown that TrajectoryNet can approximate dynamic OT. Instead of the hard boundary constraint p(-, to) = p and p(-, ti) = v as in Dynamic OT, a hard constraint at to exists, but relax the constraint at ti minimizing DKL(P(-, ti), v. Using the standard Lagrangian formulation, for large enough penalization these forms are equivalent. This is formalized in the following Theorem.

[0221] Theorem 1 [0222] (Theorem 4.1 [Tong, A., Huang, J., Wolf, G , van Dijk, D. & Krishnaswamy, S. Trajectory net: A dynamic optimal transport network for modeling cellular dynamics. In Proceedings of the 37th International Conference on Machine Learning (2020)]). With time varying field such that and subject to the continuity Eq. 13. There exists a sufficiently

large I such that

[0223] By parameterizing the field f with a neural network, this theorem allows us to optimize for the dynamic optimal transport flows with deep learning techniques. For a full proof see [Tong, A., Huang, J., Wolf, G., van Dijk, D. & Krishnaswamy, S. Trajectory net: A dynamic optimal transport network for modeling cellular dynamics. In Proceedings of the 37th International Conference on Machine Learning (2020)], a proof sketch was provided here. This result can be seen by viewing X as a Lagrange multiplier over the relaxing the constraint p(-, ti) = v. Here, this constraint can be relaxed with any proper divergence, but chose the KL-Divergence because of its link to maximum likelihood. As approximately satisfying the

equality constraint. See [Tong, A., Huang, J., Wolf, G., van Dijk, D. & Krishnaswamy, S. Trajectory net: A dynamic optimal transport network for modeling cellular dynamics. In Proceedings of the 37th International Conference on Machine Learning (2020)] for a more detailed analysis.

Modelling continuous trajectories with TrajectoryNet

[0224] Cells were modelled as points in a continuous state-space over a time interval

We assume that cells evolve via a drift through state space described by an

ordinary differential equation (ODE):

[0225] where x(t) G C denotes the state of a cell at time t and /is a time varying vector field. Additionally, cells may grow or die according to a time varying proliferation rate modeled as

represents a state where on average cells are proliferating and

conversely g(x, t) < 1 denote states and times where cells are dying. At a population level, this can be modeled according to the continuity equation:

[0226] which holds at every point

[0227] To learn the underlying dynamics from a cross-sectional observation model was desired.

A set of discrete samples

were available of the measurements from the continuous population of cells p(tt) present at k discrete timepoints in where

X(ti) consists of n observations Here, because of the destructive

nature of the measurements, it was assumed if a cell is observed at t, it is not observed at any

[0228] Nevertheless, the goal of TrajectoryNet is to learn the underlying dynamics of this space described by / and g in equation Eq. 17. To accomplish this,/ and g were parameterized by neural networks fe and ge, and then trained using a KL-divergence objective implemented using maximum likelihood in the style of continuous normalizing flows [Chen, R. T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31 (2018). 1806.07366],

[0229] The full objective of TrajectoryNet can be decomposed into four parts as appears in Eq. 3 repeated here for clarity:

Equation 18 where

q

Performing Unbalanced Dynamic Optimal Transport Using A Proliferation Rate Neural Network

[0230] To model cell proliferation and death, the system was modelled as an unbalanced optimal transport problem. The standard optimal transport problem is balanced in that all mass in the source distribution is moved to match the target distribution, however this does not model cell proliferation and cell death. The unbalanced optimal transport problem allows for the creation and destruction of mass (at a cost) thereby modelling cell proliferation and cell death. For a source measure // and a target measure v over a domain denote the set of all

joint measures on whose marginals are n and v. Then the squared Wasserstein-Fischer-

Rao (WFR) metric is f [

[0231] where p(x, t) is a time-dependent density that interpolates between u and v,f and g are time-dependent vector and scalar fields respectively modeling the drift, proliferation and destruction of mass, and p,f and g satisfy the continuity equation

q

[0232] which matches Eq. 17 substituting the general measure M for a density p. This enforces balance of mass at a particular point in space time (x, f) i.e. that the change in mass dtp balances with the amount of mass leaving the point V • (pf) and the amount of mass growth at that location pg. The WFR is related to static formulations of unbalanced optimal transport in [Liero, M., Mielke, A. & Savare, G. Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures. Inventiones mathematicae 211, 969-1117 (2018); Chizat, L., Peyre, G., Schmitzer, B. & Vialard, F.-X. Unbalanced optimal transport: Dynamic and Kantorovich formulations. Journal of Functional Analysis 274, 3090-3123 (2018)]. While the dynamic formulation is interesting theoretically, it presents numerical challenges, and thus most work has focused on discrete and static unbalanced transport. Previous work has explored the benefits of the unbalanced formulation in terms of robustness, generalization, and more accurate modeling of cell trajectories [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928- 943. e22 (2019)], however this method could only be applied to Euclidean geodesics, where more complex geodesics based on other knowledge about the system (Ldenstty) were learned. These works focus on the discrete case, where the optimization is over a fixed number of cells and cannot be easily generalized to new cells. While the ancestors, descendants, and proliferation rate for measured cells are modeled, these parameters cannot be easily extended to any unmeasured cell states. TrajectoryNet learns continuously parameterized neural networks which allows one to optimize for the proliferation rate directly with the unbalanced transport portion of the loss in Eq.

18.

[0233] The distribution matching, energy, balance, and proliferation losses approximate the WFR unbalanced dynamic optimal transport problem were now described. The TrajectoryNet unbalanced OT loss can be thought of as a relaxation of the WFR optimization where the two constraints p(x, T) = v(x), and are relaxed as regularizations. As the KL

divergence term — ► 0 the population distribution at time T approaches the observed

density v. Similarly, as the Imbalance term goes to zero, the constraint is satisfied.

The WFR distance can be formulated in dynamic terms following [43, 44] as

[0234] such that This can be thought of as minimizing the

cost following the average single point over time rather than minimizing the cost with respect to a fixed position in space time. From Eq. 26 it is easier to see the connection to the TrajectoryNet loss. L_energy and L_{proliferation} make up the main terms within the integral and and

^balance enforce the constraints. Since the loss is evaluated by drawing samples from p and integrating backwards in time if all of these losses are minimized with the correct relative X regularizations weights such that and L_balance go to zero while ^energy and

^proliferation are minimized, then the WFR distance is computed.

Practical Considerations for Efficient Computation

[0235] One can approximate the first part of this continuous time equation using a Riemann sum as

Equation 27

[0236] This requires a forward integration using a standard ODE solver to compute as shown in [Chen, R. T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31 (2018). 1806.07366], If the case where the divergence is small is considered, then this can be combined with the standard backwards pass for even less added computation. Instead of penalizing

on a forward pass, the same quantity on a backwards pass was penalized.

[0237] One can also compute the change in the data log-likelihood using an augmented neural ODE. Here, an additional dimension was added that represents the relative change in log probability over time.

[0238] This means that the TrajectoryNet loss can be computed in a single integration of an ODE solver in dimensions.

[0239] For the energy regularization L_nergy, in practice it was found that both a penalty on the Jacobian or additional training noise helped to get straight paths with a lower energy regularization similarly to [Finlay, C., Jacobsen, J.-H., Nurbekyan, L. & Oberman, A. M. How to train your neural ODE: The world of Jacobian and kinetic regularization. ICML (2020).

2002.02798], It was found that a value of A- large enough to encourage straight paths, but unsurprisingly also shortens the paths undershooting the target distribution. To counteract this, a penalty was added on the norm of the Jacobian of f as used in [Vincent, P , Larochelle, H , Lajoie, I., Bengio, Y. & Manzagol, P.-A. Stacked Denoising Autoencoders: Learning Useful Representations in a Deep Network with a Local Denoising Criterion. Journal of Machine Learning Research 3371-3408 (2010); Rifai, S., Vincent, P., Muller, X., Glorot, X. & Bengio, Y. Contractive Auto-Encoders: Explicit Invariance During Feature Extraction. In Proceedings of the 29th International Conference on Machine Learning, 833-840 (201 1)]. Since /represents the derivative of the path, this discourages paths with high local curvature, and can be thought of as penalizing the second derivative (acceleration) of the flow. The energy loss is then

[0240] wher is the Frobenius norm of the Jacobian of f. Without energy

regularization TrajectoryNet paths are unconstrained. However, with energy regularization the paths of the optimal map was approached. The energy loss gives control over how much to penalize indirect, high energy paths.

[0241] Optimal transport is traditionally performed between a source and target distribution. Extensions to a series of distributions is normally done by performing optimal transport between successive pairs of distributions as in [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928-943. e22 (2019)]. This creates flows that have discontinuities at the sampled times, which may be undesirable when the underlying system is smooth in time as in biological systems. The dynamic model approximates dynamic OT for two timepoints, but by using a single smooth function to model the whole series the flow becomes the minimal cost smooth flow over time.

Enforcing Transport on a Manifold

[0242] Methods that display or perform computations on the cellular manifold often include an implicit or explicit way of normalizing for density and model data geometry. PHATE [Moon, K. R. et al. Visualizing structure and transitions in high-dimensional biological data. Nature Biotechnology 37, 1482-1492 (2019)] uses scatter plots and adaptive kernels to display the geometry. To constrain the flows to the manifold geometry is desirable, but not to its density. The flow was penalized such that it is always close to at least a few measured points across all timepoints through the following function:

[0243] This can be thought of as a loss that penalizes points until they are within h Euclidean distance of their k nearest neighbors. The values h = 0.1 and k = 5 were used in all of the disclosed experiments. The L_density was evaluated on an interpolated tim every batch.

Training

[0244] For simplicity, the neural network architecture of TrajectoryNet consists of three fully connected layers of 64 nodes with leaky ReLU activations. It takes as input a cell state and time and outputs the derivative of state with respect to time at that point. To train a continuous normalizing flow, access was needed to the density function of the source distribution. Since this is not accessible for an empirical distribution, an additional Gaussian at to was used, defining the standard Gaussian distribution, where Pt(x) is the density function at time t.

[0245] For a training step, samples

were drawn and the loss calculated with a single backwards integration of the ODE. While there are a number of ways to computationally approximate these quantities, a parallel method was used to iteratively calculate the log P based on log - To make a backward pass through all timepoints, the final

timepoint was the start, and integrated the batch to the second to last timepoint, concatenated these points to the samples from the second to last timepoint, and continued till to, where the density was known for each sample. It was noted that this can compound the error especially for later timepoints if k is large or if the learned system is stiff, but gave significant speedup during training. To sample from was sampled then the adjoint method was used to perform

the integration This is similar to other continuous normalizing

flows where the model is trained backwards in time, but the invertible dynamics of the ODE to sample was used.

Pre-training of the Growth Network [0246] To pretrain the growth network, a simple and computationally efficient method was used that adapts discrete static unbalanced optimal transport to the disclosed framework in the continuous setting. A network was trained, which takes as input a cell

state and time pair and produces a proliferation rate of a cell at that time. This is trained to match the result from discrete optimal transport. It was noted that adding proliferation rate regularization in this way does not guarantee conservation of mass. One could normalize o

be a probability distribution during training, e g., as However, this

now requires an integration over

which is too computationally costly. Instead, the equivalence of the maximum likelihood formulation over a fixed proliferation function g was used and normalized it after the network was trained.

Comparison of TrajectoryNet Algorithm

[0247] Next, the disclosed TrajectoryNet algorithm was benchmarked on both simulated and real single-cell data. TrajectoryNet was compared to RNA Velocity [La Manno, G. et al. RNA velocity of single cells. Nature 560, 494-498 (2018); Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. BioRxiv 820936 (2019)] and Diffusion Pseudotime [Haghverdi, L., Biittner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016)].

[0248] Existing methods in cellular trajectory learning attempt to infer a trajectory within one timepoint [La Manno, G. et al. RNA velocity of single cells. Nature 560, 494-498 (2018); Haghverdi, L., Biittner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016); Saelens, W., Cannoodt, R., Todorov, H. & Saeys, Y. A comparison of single-cell trajectory inference methods. Nature Biotechnology 37, 547-554 (2019)], or interpolate linearly between two timepoints [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928-943. e22 (2019); Yang, K. D. & Uhler, C. Scalable Unbalanced Optimal Transport Using Generative Adversarial Networks. In 7th International Conference on Learning Representations, 20 (2019)], but TrajectoryNet can interpolate non-linearly using information from more than two timepoints. TrajectoryNet has advantages over existing methods in that it: can interpolate by following the manifold of observed entities between measured timepoints, thereby solving the static-snapshot problem, and can create continuous-time trajectories of individual entities, giving researchers the ability to follow an entity in time.

Comparing TrajectoryNet to other trajectory inference tools on synthetic data

[0249] To test the ability of TrajectoryNet to follow a manifold, two toy datasets were generated, an arch dataset representing a single trajectory and a branching tree dataset representing a diverging structure (Figure 7A). In both cases the datasets are designed to represent lowdimensional (ID) curved manifolds in a higher dimensional ambient space (2D). Both datasets are split into three timepoints supplied as training data and

data ? is held out as test data. In Figure 7B, the EMD was compared between the predicted data a with the ground truth data X1/2. Here the

TrajectoryNet approach was compared to six baseline methods. The previous baseline, which predicts the previous timepoint as The next baseline, which predicts the next

timepoint as t_1/2i'.e The optimal transport interpolant (OT), which predicts the

McCann interpolant of the exact optimal transport plan between

as used in [Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming. Cell 176, 928-943. e22 (2019)]. The random transport interpolant, which predicts the McCann interpolant of a random transport plan between Xo and Xi. An RNA-velocity method RNA_v, which takes the ground truth instantaneous velocity dXo of data at Ao, finds the optima

which minimizes mint E then

predicts Finally, TrajectoryNet was compared to a standard continuous

normalizing flow model [Chen, R. T. Q., Rubanova, Y., Bettencourt, J. & Duvenaud, D. Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31 (2018). 1806.07366] CNF with energy regularization but without additional manifold regularization. TrajectoryNet was found to perform the best in terms EMD (A-^2, £1/2) followed by the random and optimal transport McCann interpolations. This is because TrajectoryNet was able to follow the manifold while other methods were not. [0250] Next, the learned trajectories and vector fields were evaluated qualitatively to check the ability of various methods to follow the cell manifold. As seen in Figure 7C the paths for TrajectoryNet are able to follow the curve of the manifold as compared to an energy regularized CNF, which follows the straight line paths i.e. the dynamic optimal transport paths for a Euclidean geometry. While these methods define an interpretable vector field everywhere, the other baselines do not. In Figure 7C, the imputed distribution of standard OF was visualized in red, a random interpolation in purple, and RNA-velocity in brown.

Comparing TrajectoryNet to RNA-velocity and Pseudotime on real scRNAseq data

[0251] TrajectoryNet along with other optimal transport methods contrast in a few notable ways from RNA-velocity or Pseudotime-based approaches to inferring cell time. In Figure 8, the differences in these approaches were

shown on the disclosed tumorsphere dataset. TrajectoryNet learns a distribution flow p(x, f) where the marginal p(x, t = 4)

approximates the observed distribution at time t. This contrasts to a Pseudotime approach which assigns a Pseudotime (or developmental time to each cell. While there are a large variety of Pseudotime inference algorithms [Saelens, W ., Cannoodt, R., Todorov, H. & Saeys, Y. A comparison of single-cell trajectory inference methods. Nature Biotechnology 37, 547-554 (2019)], many of them do not take the observed cell time into account. Here the disclosed results were compared against the popular Diffusion Pseudotime algorithm [Haghverdi, L , Buttner, M., Wolf, F. A., Buettner, F. & Theis, F. I. Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods 13, 845-848 (2016)]. In Figure 8 A, the inferred pseudotime of different algorithms were qualitatively evaluated by coloring the PHATE embedding of the disclosed tumorsphere data by the inferred cell time (normalized to between lie between zero and one) for TrajectoryNet, scVelo, and Diffusion Pseudotime. Interestingly, scVelo reverses the ordering of the cells, assigning cells at day 0 a later Pseudotime than all other timepoints on average. It was suspected this was because the days here are relatively far apart, and RNA-velocity analysis can only predict cell trajectories over shorter time scales. It was also seen that Diffusion Pseudotime misorders the timepoints. While cells for day 0, 2,12, and 18 are correctly ordered, cells from day 30 are assigned pseudotime values that are on average earlier than those in days 12 and 18. TrajectoryNet in contrast is able to infer trajectories over longer time scale by taking observed cell time into account and (optionally) RNA-velocity estimates at each cell. [0252] This was quantitatively evaluated next in Figure 8B. Here, the distribution of inferred psuedotime was plotted for each day for different methods. TrajectoryNet orders the days correctly while RNAVelocity and Diffusion Pseudotime methods do not order the days monotonically.

[0253] In this dataset, RNA velocity data was not incorporated into the TrajectoryNet loss because it was found to be unreliable. In Figure 8C, it was shown the scVelo inferred velocity flow plotted separately on each day on three different 2D projections: A principal component analysis (PCA) projection, a UMAP [Mclnnes, L., Healy, J. & Melville, J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (2020). 1802.03426] projection, and a PHATE [Moon, K. R. et al. Visualizing structure and transitions in highdimensional biological data. Nature Biotechnology 37, 1482-1492 (2019)] projection. Cells were colored by cell phase for reference. It was found that depending on the projection, the inferred order of cells is different. This suggests RNA-velocity type analysis may not be reliable in this dataset.

Comparing TrajectoryNet Pipeline and Total Granger Causal Score inference against static gene network inference methods on synthetic data

[0254] The disclosed gene network inference platform first infers the trajectories of individual cells via TrajectoryNet then performs total granger causality score (TGCS) analysis on an average trajectory for each endpoint. The belief is that accurately inferring the dynamics information with TrajectoryNet will more accurately predict gene-gene interactions than approaches that do not include this information. To establish that this approach can recover the ground truth structure on simulated single-cell data with a known regulatory structure, TGCS analysis was compared to three other static methods for inferring regulatory structure from single-cell data, DREMI [Krishnaswamy, S. et al. Conditional density-based analysis of t cell signaling in single-cell data. Science 346, 1250689 (2014)], Pearson correlation, and Spearman correlation across genes. The BoolODE package was used [Pratapa, A., Jalihal, A. P., Law, J. N., Bharadwaj, A. & Murali, T. M. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature Methods 17, 147-154 (2020)] to simulate multiple dataset structures including linear, cyclic, and bifurcating structures (Figure 9B). The predicted regulatory structures were tested to see how well they match the true regulatory structure of the simulation across experiments using area under the receiver operator characteristic curve (AUCROC). For these experiments, evaluation of the diagonal was excluded, assuming that there were self loops in a time varying system. Furthermore, whether there was an edge was evaluated, and did not differentiate between positive and negative regulatory relationships. Since the performance of the disclosed TGCS was of interest here, the ground truth simulated trajectories were used and ignored the post-processing steps that add additional noise to the data. Across unique runs of each algorithm on 100 randomly generated data of each type, it was found that the disclosed TGCS analysis consistently infers gene-gene regulatory relationships more accurately than other approaches (Figures 9C-D).

Other Computational Methods Details

[0255] scRNA-seq data from 5 samples, were processed with 10X and CellRanger pipeline according to the following steps. Sample demultiplexing and read alignment to the NCBI reference GRCh38 was completed to map reads using CellRanger. Prefdtering was performed using parameters in scprep vl.0.3. Cells that contained at least 1,000 unique transcripts were kept for further analysis to generate a cell by gene matrix containing 17,983 cells and 16,983 genes. Normalization was performed using default parameters with LI normalization, adjusting total library side of each cell to 10,000. Any cell expressing mitochondrial genes greater than 10% of their overall transcriptome were removed. Raw data fdes for scRNA-seq data will be available for download through GEO under an accession number to be assigned with no restrictions on publication.

In vivo scRNA-seq sequencing and pre-processing

[0256] scRNA-seq data from four replicates for the primary tumor and lung metastasis were aligned to GRCh38 and read mapping was completed with CellRanger. Cells were retained that mapped strongly to the human genome and expressed at least 1000 genes. Genes expressed in at least 20 cells were then retained and filtered the data to remove cells with fewer than 2000 counts and more than 5000 total counts. Finally, library size was normalized and square-root transformed the data. For each tissue, an MNN kernel was used to build a cellular graph with batch correction between the replicates. MAGIC was then used [van Dijk, D. et al. Recovering gene interactions from single-cell data using data diffusion. Cell 174, 716-729. e27 (2018)] to transform the MNN graph into the data space and ran TruncatedSVD to reduce the dimensionality to 100.

[0257] For calculating the epithelial-to-mesenchymal score, the approach was adapted from [Chakraborty, P., George, J. T., Tripathi, S., Levine, H. & Jolly, M. K. Comparative study of transcriptomics-based scoring metrics for the epithelial-hybrid-mesenchymal spectrum. Frontiers in Bioengineering and Biotechnology 8 (2020)] to calculate gene correlation with the mean expression of epithelial markers EPCAM, CDH3, CTNNB1, KRT8, KRT18, as CDH1 was filtered out in above quality control steps.

In vivo spatial sequencing and pre-processing

[0258] Spatial transcriptomics was conducted using 10X Visium Spatial Gene Expression Slide and Reagent Kit, 16 rxns (PN- 1000184), according to the protocol detailed in document CG000239RevD for the TNBCs and CG000239RevE for the Xenografts, available in 10X Genomics demonstrated protocols.

[0259] Next, scMMGAN was leveraged [Amodio, M. et al. Single-cell multi-modal GAN reveals spatial patterns in single-cell data from triple-negative breast cancer. Patterns (N Y) 3, 100577 (2022)] to generate single-cell expression values for each spatial voxel in the same data space as the single-cell data. With scMMGAN, a generator was used consisting of three internal layers of 128, 256, and 512 neurons with batch norm and leaky rectified linear unit activations after each layer, and a discriminator consisting of three internal layers with 1,024, 512, and 256 neurons with the same batch norm and activations except with minibatching after the first layer. The geometry-preserving correspondence loss with a coefficient of 10, cycle-loss coefficient of 1, learning rate of 0.0001, and batch size of 256 was used.

Transcriptional network generation

[0260] Transcriptional networks were built to visualize direct and indirect regulatory interactions of the EMT and MET trajectory using core transcription factors (transcription factors) identified by TrajectoryNet. Briefly, TRUSST v2 (Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining) database was used as the canvas for all known regulatory relationships across the human genome. Transcriptional interactions were filtered for the 23-core MET and 37-core EMT transcription factors derived for Gene Clusters 1-5. The resultant network (Fig 4H) was organized to enable the visualization of (i) direct interactions between the core MET and EMT transcription factors, (ii) common nodes that share regulatory relationships with the core MET and EMT transcription factors as well as (iii) unique nodes that shared direct regulatory relationships with either the core MET or core EMT transcription factors.

[0261] To develop a detailed transcriptional network describing the MET, the TRUSST v2 regulatory canvas was used to filter direct regulatory relationships of the 23 core MET transcription factors. Nine of the 23 transcription factors had known regulatory relationships on TRUSST v2, which allowed mapping the genes regulated by these core nodes. Next, the Epithelial Mesenchymal Gene Database, dbEMT 2.0 was used to distinguish known EMT genes (diamond) from novel EMT genes (ellipse) in the network. The resultant network was then organized based on the temporal gene clusters (2-5) and their regulatory relationships to observe direct and indirect crosstalk across the clusters forming the MET network (Fig 5C). All networks were built using Cytoscape 3.9.0

Enrichment analysis

[0262] Gene lists for the MET and EMT sub-networks (Figures 14A, 14E) were extracted using Cytoscape 3.9.0. These gene lists were subjected to analysis on Metascape to identify enriched biochemical pathways [Zhou, Y. et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nature Communications 10 (2019)].

Tumor sphere assay

[0263] Single-cell suspensions were plated in ultra-low attachment 96-well plates (Corning # CLS3474, New York, USA) at low densities optimized to ensure tumorspheres arose from single anchor-independent cells. HCC38 CD44^A' cells were seeded in 100 pl at 100 cells/well. Cell-line specific serum-free media was supplemented with 1% (v/v) penicillin/streptomycin, 20 ng/ml EGF 20 ng/ml FGFb, 4 pg/mL heparin, lx B27, and 1% (v/v) methylcellulose (Sigma-Aldrich). Fresh media was topped up every 5 days by adding 50 pl per well of the appropriate tumorsphere media. Tumorspheres were counted at day 30 under 4x magnification and averaged 10 tumorspheres/well .

[0264] For 3D immunostaining, tumorspheres were fixed with 10% formalin (Australian Biostain Pty Ltd) for 1 hour at RT in a gently rocking rotator and washed in TBS (3 x 15 min). Tumorspheres were then permeabilized with 100% methanol for 10 minutes at 4°C and washed in TBS (3 x 15 min). Tumorspheres were then blocked in TBS 5% BSA, 10% Horse Serum and 0.1% Triton O/N at 4 degrees with rotation. Following incubation with primary antibodies ESRRA (Cell Signaling Technology, E1G1J, 1 :200, Cat no. 13826); ZEB1 (Santa Cruz, #H-102, 1 :200); CDH1 (BD Biosciences/(36/E-Cadherin, Cat. no. 610181, 1:200) diluted in blocking buffer O/N at 4°C in rotation, spheres were washed (4 x 30 min) in TBS and stained with the appropriate fluorophore-conjugated secondary antibodies (1 :500) (Ms Cy3 (#M30010), Rb 647 (#A32722), Ms 488 (#A11001), Rb 488 (#711546152), ThermoFisher Scientific) and DAPI O/N at 4°C. Tumorspheres were washed with TBS (4 x 30 min) then resuspended in 20 pl of mounting media (ProLong Diamond Antifade Mounting Media (ThermoFisher Scientific)) and mounted between a glass slide and a coverslip spaced by tape.

[0265] Labelled tumorspheres were imaged using confocal microscopy (Leica DMI 6000 SP8 with 40x (NA 1.3) or 63x (NA 1.4) oil objectives or a Nikon AIR confocal with 20x Plan Apochromat air objective (NA 0.75) at 2x zoom using an HD25 resonance scanner) using identical acquisition settings (optimized per protein marker) for all time points. Quantitative image analysis was performed using CellProfiler (v4.2.1, [Carpenter, A. E. et al. Cellprofiler: image analysis software for identifying and quantifying cell pheno-types. Genome Biology 7, R100 (2006)]) to segment individual cells via maximum cross-entropy -based threshold detection of nuclei (DAPI) and cell bodies (sum of all channels). Mean intensity per cell was measured for each marker, with per cell image and quantitative data integrated via a Knime software (v4.6.4, [Berthold, M. R. et al. KNIME: The konstanz information miner. In Data Analysis, Machine Learning and Applications, 319-326 (Springer Berlin Heidelberg, 2008)]) for visualization [Bryce, N. S. et al. High-content imaging of unbiased chemical perturbations reveals that the phenotypic plasticity of the actin cytoskeleton is constrained. Cell Systems 9, 496-507. e5 (2019).; Lock, I. G. et al. Visual analytics of single cell microscopy data using a collaborative immersive environment Tn Proceedings of the 16th ACM STGGRAPH International Conference on Virtual -Reality Continuum and its Applications in Industry (ACM, 2018)].

FACS analysis

[0266] Bulk HCC38 cell lines were cultured in 2D tissue culture dishes. For isolating CD44^hi cells from this bulk population, cells were trypsinized and stained with CD44 antibody (BD Biosciences anti-human CD44-PE-cy7 (1.800)) for 25 min at 4C. CD44^hi cells were sorted on BD Aria III. Sorted cells were cultured in media supplemented with 0.1% (v/v) gentamicin and 1% (v/v) antibiotic-antimycotic for at least two passages to avoid contamination. Multiple rounds of FACS enrichment were performed on these expanded cultures until pure CD44hi populations were isolated. To identify EPCAM^+/-, CAV1^+/-, EPCAM/CAV1^+/+ populations, HCC38 CD44*' were subjected to FACS sorting using Anti-CD326 (EPCAM) (Invitrogen #53-8326-42) and Anti-Caveolin 1 (BD Biosciences). Data acquisition was performed using BD Aria III and FACSDiva software (BD Biosciences and data analysis was performed using Flowjo XI 0.7.1.

ESSRA knockdown and validation

[0267] Preliminary validation of the MET network included an in-vitro knockdown of ESRRA in HCC38 CD44^ft' cells (n=3 biological replicates). Predesigned siRNA specific to ESRRA and scrambled siRNA were purchased from Integrated DNA Technologies, USA (TriFECTa® RNAi Kit, Design ID hs.Ri.ESRRA.13). HCC38 CD44^hi cells were seeded in 24 well plates at a density of 9000 cells/well and ESRRA knockdown was performed using 10 nM of pooled siRNA (hs.Ri.ESRRA.13.1, hs.Ri.ESRRA.13.2 and hs.Ri.ESRRA.13.3) using Lipofectamine RNAiMax (ThermoFisher Scientific, USA) as per the manufacturer’s protocol. Cells were harvested for protein extraction 48hrs post siRNA transfection to study the downstream effect on CDH1 expression upon ESRRA knockdown using western blot.

[0268] In parallel, HCC38 CD44^hi cells were also treated with an ESRRA inhibitor ((2- Aminophenyl)(l-(3-isopropylphenyl)-lH-l,2,3-triazol-4-yl)methanone) (BLD Pharm, China, Cat no. BD01201330) [referred to as Compound 14 (C14) here on] at 5 and 10 pM concentrations. Vehicle controls were treated with DMSO. Cells were harvested 48 hours post-treatment and processed for protein extraction which were studied for changes in ESRRA and CDH1 expression using western blot (n=3 biological replicates).

Western Bloting

[0269] Proteins were extracted from control and treated cells using ice cold modified R1PA buffer (50 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 0.1% SDS with supplemented with IX protease and phosphatase inhibitor cocktails). The lysates were sonicated using a QSONICA Q55 probe sonicator at 50 kHz for 20 seconds in an ice bath. This whole cell lysate was centrifuged at 14,000 x g for lOmin at 4°C and stored in -80°C until further use. Proteins were quantified using the Pierce™ BCA Protein Assay Kit (Cat. 23227, Thermo Fisher Scientific, USA) as per the manufacturer’s instructions. 20 pg of whole cell lysates were separated on ID SDS-PAGE using NuPAGE™ 4 to 12%, Bis-Tris (Cat. NP0322BOX, ThermoFisher Scientific, USA). Proteins separated on the gel were transferred onto nitrocellulose membrane using Trans-Turbo Transfer system (BioRad Laboratories, USA) at 1.3V for 7min. The membrane was blocked using 5% milk in TBS Ihr. The membrane was next incubated overnight at 4°C with primary antibodies against the protein of interest [1 : 1000 ESRRA (Cell Signaling Technologies, USA, Cat. no. 13826), 1 :1000 CDH1 (Cell Signaling Technologies, USA, Cat. no. 3195S), 1 :5000 GAPDH (Cell Signaling Technologies, USA, Cat. no. 97166S), 1 :5000 P-Actin (Cell Signaling Technologies, USA, Cat. no. 3700S)]. After incubation with the primary antibody, membranes were washed with IX TBST for 3 x lOmin on a rocker at room temperature. Next, membranes were incubated with appropriate HRP- conjugated secondary antibodies. Washing steps were repeated, and the membrane was developed using ECL substrate (Western Lightning™ Ultra, Perkin Elmer or Clarity Western ECL Substrate) and scanned using Fusion FX Vilber Lourmat scanner. Each Western blot experiment was performed using three biological replicates, to calculate the statistical significance (p-value) of relative fold-change in expression. GAPDH or P-Actin was used to normalize the relative fold-change expression value of ESRRA and CDH1. The signal intensity of the bands in western blots was quantified using Image Studio Lite version 5.2 (LLCOR Biotechnology, USA). [0270] The disclosures of each and every patent, patent application, and publication cited herein are hereby incorporated herein by reference in their entirety. While this invention has been disclosed with reference to specific embodiments, it is apparent that other embodiments and variations of this invention may be devised by others skilled in the art without departing from the true spirit and scope of the invention. The appended claims are intended to be construed to include all such embodiments and equivalent variations.

Claims

CLAIMS What is claimed is:

1. A method of estimating a dynamic molecular program of a population of cells, comprising: providing a set of at least two static snapshots of a population of cells undergoing a state transition at a corresponding set of at least two time indices to a neural network; calculating a set of possible population flows between the at least two time indices based on the at least two static snapshots; negatively weighting any of the set of population flows which are unrealistic; and inferring an estimated population flow of the cells between the set of static snapshot data by selecting a population flow from the set of possible population flows with the neural network.

2. The method of claim 1, further comprising an Ordinary Differential Equation (ODE) solver to form a neural ODE system.

3. The method of claim 1, wherein the state transition is a mesenchymal-to-epithelial transition (MET), or epithelial-to-mesenchymal transition (EMT).

4. The method of claim 1, wherein the calculating step further comprises the steps of: approximating a time-varying derivative parametrized by a set of network weights and biases; and integrating the time-varying derivative across the at least two time indices to calculate a possible population flow of the set of possible population flows.

5. The method of claim 4, further comprising the step of limiting a magnitude of the timevarying derivative across the at least two time indices to a predetermined threshold.

6. The method of claim 1, wherein the set of at least two static snapshots comprise three- dimensional tumorsphere data.

7. The method of claim 1 , further comprising the step of identifying shared and trajectoryspecific molecular programs using a time-lapsed causality analysis.

8. The method of claim 1, wherein the at least two time indices of the at least two static snapshots are separated by at least 24 hours.

9. The method of claim 1, wherein the dynamic molecular program is a gene regulatory network.

10. The method of claim 1, further comprising the step of calculating an interpolated snapshot of the cell population at a third time index based on the inferred population flow.

11. The method of claim 1, wherein the set of population flows which are unrealistic comprise energy inefficient biological pathways.

12. A method of preventing a mesenchymal cell from differentiating into an epithelial cell comprising contacting the cell with one or more modulators of one or more molecules involved in the mesenchymal-to-epithelial transition (MET) transcriptional network.

13. The method of claim 12, wherein the one or more molecules in the MET transcriptional network are one or more selected from the group consisting of: estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAU (SNAU), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1).

14. A method of directing one or more cells in a population of cells to a transition state, comprising the steps of: obtaining a population of cells from the subject; identifying a target state transition for the population of cells to undergo; administering at least one modulator of at least one molecule within a MET transcriptional network of the population of cells thereby directing the population of cells to the targeted state transition.

15. The method of claim 14, wherein the at least one molecule in the MET transcriptional network is selected from the group consisting of: estrogen related receptor alpha (ESRRA), aryl hydrocarbon receptor (AHR), aryl hydrocarbon receptor nuclear translocator (ARNT), estrogen receptor 1 (ESRI), transcription factor Jun (JUN), androgen receptor (AR), zinc finger E-box binding homeobox 1 (ZEB1), zinc finger protein SNAU (SNAI1), zinc finger protein SNAI2 (SNAI2), and cadherin 1 (CDH1).

16. The method of claim 14, wherein the targeted state transition is energy optimal.