Analyses of Multi-dimensional Single Cell Trajectories Quantify Transition Paths Between Nonequilibrium Steady States

A problem ubiquitous in almost all scientific areas is escape from a metastable state, or relaxation from one stationary distribution to a new one. More than a century of studies lead to celebrated theoretical and computational developments such as the transition state theory and reactive flux formulation. Modern transition path sampling and transition path theory focus on an ensemble of trajectories that connect the initial and final states in a state space . However, it is generally unfeasible to experimentally observe these trajectories in multiple dimensions and compare to theoretical results. Here we report and analyze single cell trajectories of human A549 cells undergoing TGF-β induced epithelial-to-mesenchymal transition (EMT) in a combined morphology and protein texture space obtained through time lapse imaging. From the trajectories we identify parallel reaction paths with corresponding reaction coordinates and quasi-potentials. Studying cell phenotypic transition dynamics will provide testing grounds for nonequilibrium reaction rate theories. INTRODUCTION Transition path theories and transition path sampling techniques have been developed to investigate transition dynamics between equilibrium and nonequilibrium attractors of nonlinear dynamical systems . Basically one divides the state space into regions containing the initial (A) and final (B) attractors, and an intermediate (I) region. A reactive trajectory is one that originates from region A, and enters region I then B before re-entering region A (Fig. 1A). The reactive trajectories form an ensemble of transition paths that connect regions A and B, and those that fail to arrive at B before returning to A form an ensemble of nonreactive trajectories. The ensembles of reactive and nonreactive trajectories are subject to various theoretical analyses. Central to rate theories is a one-dimensional variable called the reaction coordinate (RC, denoted as s), which is a geometric parameter describing the progression along a reaction path. Continuous efforts have been made on defining an optimal RC for a dynamical system. In practice, one often approximates the RC as a string of discrete points, and divides the generally multi-dimensional state space into an 1-D array of Voronoi cells that connects the initial and final states, so the system transits between neighboring Voronoi cells (Fig. 1B). Theoretically, for a stochastic system, one can define an action S, e.g., the Onsager-Machlup action, for a reactive trajectory, so the probability of observing this trajectory is ∝exp(−S) . Minimization of the action through variational analysis leads to one or more most probable paths (MPPs), and reactive trajectories mainly fluctuate around each MPP forming a tube in the state space that connects A and B. Correspondingly algorithms such as finite temperature string method have been developed to define a RC as the center of the tube. Figure1 Schematic diagram of defining reaction coordinate from reactive trajectory ensemble. (A) Example reactive and nonreactive trajectories in a quasi-potential system. Also shown is a valley (tube) connecting regions A and B that most reactive trajectories fall in. (B) Discrete representation of a 1-D reaction coordinate (colored dots) on the filled contour map with corresponding Voronoi cells. The cyan line is a reactive trajectory that starts from A and ends in B. Advances of the field of transition rate theories would undoubtedly benefit from direct experimental recording of the multi-dimensional trajectories, which unfortunately is technically challenging especially for a chemical system. To address this challenge, we studied instead a cellular process in which cells change between different phenotypes. A cell is a nonlinear dynamical system, and a cell phenotype is a stable or metastable state of the system. Cells of multicellular organisms can be induced to undergo cell phenotypic transitions (CPTs) such as differentiation and reprogramming. RESULTS Applying rate theories to cell phenotypic transitions first requires a mathematical representation of the cell states and cell trajectories. Mathematically one can represent a cell state by a point in a multi-dimensional space defined by gene expression or other cell properties. Noticing that cells have phenotype-specific morphological features that can be monitored even with transmission light microscopes, recently we developed a framework that defines cell states in a combined morphology and texture feature space. The framework allows one to trace individual cell trajectories during a CPT process through live-cell imaging. Figure 2 Summary of pipeline for recording and analyzing single cell trajectories. (A) Time lapse imaging of A549-Vimentin-RFP cells treated with TGF-β. (B) Deep-learning aided single cell segmentation and tracking on the acquired time lapse images. (C) Extraction of morphology and vimentin features of single cells. (D) Representation of single cell states in a multidimensional morphology/texture feature space. (E) Transition path analyses over recorded trajectories. Right: A representative single cell trajectory of EMT in the feature space. Color represents time (unit: hour). Left: the same trajectory colored by the regions in the feature space (E, I, and M) each data point resides. Reduced units are used in this and all other figures. We applied the framework to study the TGF-β induced EMT with a human A549 cell line with mCherry fused to endogenous vimentin. Vimentin is a type of intermediate filament protein commonly used as a mesenchymal marker. We performed time lapse imaging of the EMT process of A549 induced by TGF-β continuously for two days sampled every 10 minutes (Fig. 2A), and performed single cell segmentation and tracking on the acquired images (Fig. 2B). We then quantified the images and performed principal component analysis (PCA) to form a set of orthonormal basis vectors of collective variables, which include cell body shape of 296 degrees of freedom (DoF) quantified by an active shape model, and texture features of cellular vimentin distribution quantified by 13 Haralick features (Fig. 2C). Then the state of a cell at a given instant is represented as a point in the composite morphology/texture feature space (Fig. 2D), and the temporal evolution of the state forms a continuous trajectory in the space. Before TGF-β treatment, a population of cells assumes a localized stationary distribution in this 309 dimensional composite space, and most cells are epithelial. TGF-β treatment destabilizes such distribution, and the cells relax into a new stationary distribution dominated by mesenchymal cells. We recorded 196 continuous trajectories in the state space. Figure 2E (left) and Fig. S1 show such a trajectory revealing how a cell transits step-by-step from an epithelial cell with round polygon shapes and a localized vimentin distribution, to the mesenchymal phenotype with elongated spear shapes and a dispersive vimentin distribution. For transition path analyses we divided the space into epithelial (E), intermediate (I), and mesenchymal (M) regions (Fig. 2E right). With this division of space we identified a subgroup of 139 recorded single cell trajectories that form an ensemble of reactive trajectories that connect E and M by day 2. First, we set to identify the RC using a modification of the finite temperature string method 14, . The procedure starts with a trial RC to define the initial division of the Voronoi cells. One then optimizes the trial RC iteratively by minimizing the distance dispersion between the string point and sample points within each Voronoi cell. Since here we have an ensemble of continuous trajectories, we modified the iteration procedure slightly. Specifically, we minimized both the distance between the ensemble of measured reactive trajectories and the string point within each individual Voronoi cell, as well as the overall distance between each individual trajectory and the trial RC (Fig. S2, see Methods for details). The iteration procedure gives the RC of TGF-β induced EMT in A549 cells (Fig. 3A left, Fig. S3), which characterizes common features of the reactive trajectories (Fig. 3A right). Along the curve s the cell shape changes dominantly through elongation and growth (Fig. 3B top), and most of the 13 vimentin Haralick features increase or decrease monotonically and continuously over time (Fig. 3B bottom). For quantitative analysis, it is reasonable to assume that dynamics of the collective variables (x) can be described by a set of Langevin equations in the morphology/texture feature space, dx / dt = F(x)+η(t) , where F(x) is a governing vector field, and η are white Gaussian noises with zero mean. Then in principle one can reconstruct F(x) from the single cell trajectory data, F(x) = dx / dt . Specifically we restricted to reconstructing the dynamics along the RC, and the ansatz becomes a 1-D convection-diffusion process, ds / dt = −dφ / ds+η . Notice that for a 1D system even without detailed balance one can define a quasi-scalar potential φ . Figure 3C presents an enlarged view of an area in Fig. 3A, which shows the instant velocities (ds/dt) of various trajectories segments. Numerically we related the potential gradient with the mean velocity within the i-th Voronoi cell by (Fig. 3D), dφ(si ) / ds = ds / dt the i -th cell ≈ 1 NiΔt s(tα +Δt)− s(tα ) ( ) α ∑ . The sum was over all Ni recorded trajectory segments that locate within the i-th Voronoi cell at any time t (Fig. 3A & C). On the obtained curve of dφ/ds v.s. s (Fig. 3E), the zeroes correspond to stationary points of the potential. We then reconstructed the quasi-potential through integrating over s, φ(si ) = φ(s0 )+ (dφ / ds)ds s0 si ∫ (Fig. 3G). Also shown in Fig. 3G is the quasi-potential of the untreated cells, defined from the steady state distribution of untreated cells along the RC, ( ) ( ) 0 log i ss i s p s φ ∝ − (Fig. 3F).


INTRODUCTION
Transition path theories and transition path sampling techniques have been developed to investigate transition dynamics between equilibrium and nonequilibrium attractors of nonlinear dynamical systems 2,3 . Basically one divides the state space into regions containing the initial (A) and final (B) attractors, and an intermediate (I) region. A reactive trajectory is one that originates from region A, and enters region I then B before re-entering region A (Fig. 1A). The reactive trajectories form an ensemble of transition paths that connect regions A and B, and those that fail to arrive at B before returning to A form an ensemble of nonreactive trajectories. The ensembles of reactive and nonreactive trajectories are subject to various theoretical analyses.
Central to rate theories is a one-dimensional variable called the reaction coordinate (RC, denoted as s), which is a geometric parameter describing the progression along a reaction path 4 . Continuous efforts have been made on defining an optimal RC for a dynamical system 5 . In practice, one often approximates the RC as a string of discrete points, and divides the generally multi-dimensional state space into an 1-D array of Voronoi cells that connects the initial and final states, so the system transits between neighboring Voronoi cells (Fig. 1B). Theoretically, for a stochastic system, one can define an action S, e.g., the Onsager-Machlup action, for a reactive trajectory, so the probability of observing this trajectory is ∝ exp(−S) . Minimization of the action through variational analysis leads to one or more most probable paths (MPPs), and reactive trajectories mainly fluctuate around each MPP forming a tube in the state space that connects A and B. Correspondingly algorithms such as finite temperature string method have been developed to define a RC as the center of the tube 6 . Figure1 Schematic diagram of defining reaction coordinate from reactive trajectory ensemble. (A) Example reactive and nonreactive trajectories in a quasi-potential system. Also shown is a valley (tube) connecting regions A and B that most reactive trajectories fall in. (B) Discrete representation of a 1-D reaction coordinate (colored dots) on the filled contour map with corresponding Voronoi cells. The cyan line is a reactive trajectory that starts from A and ends in B.
Advances of the field of transition rate theories would undoubtedly benefit from direct experimental recording of the multi-dimensional trajectories, which unfortunately is technically challenging especially for a chemical system. To address this challenge, we studied instead a cellular process in which cells change between different phenotypes 7 . A cell is a nonlinear dynamical system, and a cell phenotype is a stable or metastable state of the system 8 . Cells of multicellular organisms can be induced to undergo cell phenotypic transitions (CPTs) such as differentiation and reprogramming.

RESULTS
Applying rate theories to cell phenotypic transitions first requires a mathematical representation of the cell states and cell trajectories. Mathematically one can represent a cell state by a point in a multi-dimensional space defined by gene expression 9 or other cell properties 10 . Noticing that cells have phenotype-specific morphological features that can be monitored even with transmission light microscopes, recently we developed a framework that defines cell states in a combined morphology and texture feature space 7 . The framework allows one to trace individual cell trajectories during a CPT process through live-cell imaging. We applied the framework to study the TGF-β induced EMT with a human A549 cell line with mCherry fused to endogenous vimentin 7 . Vimentin is a type of intermediate filament protein commonly used as a mesenchymal marker 11 . We performed time lapse imaging of the EMT process of A549 induced by TGF-β continuously for two days sampled every 10 minutes ( Fig.  2A), and performed single cell segmentation and tracking on the acquired images (Fig. 2B). We then quantified the images and performed principal component analysis (PCA) to form a set of orthonormal basis vectors of collective variables, which include cell body shape of 296 degrees of freedom (DoF) quantified by an active shape model 12 , and texture features of cellular vimentin distribution quantified by 13 Haralick features 13 (Fig. 2C). Then the state of a cell at a given instant is represented as a point in the composite morphology/texture feature space (Fig. 2D), and the temporal evolution of the state forms a continuous trajectory in the space. Before TGF-β treatment, a population of cells assumes a localized stationary distribution in this 309 dimensional composite space, and most cells are epithelial. TGF-β treatment destabilizes such distribution, and the cells relax into a new stationary distribution dominated by mesenchymal cells. We recorded 196 continuous trajectories in the state space. Figure  First, we set to identify the RC using a modification of the finite temperature string method 6,14,15 . The procedure starts with a trial RC to define the initial division of the Voronoi cells. One then optimizes the trial RC iteratively by minimizing the distance dispersion between the string point and sample points within each Voronoi cell. Since here we have an ensemble of continuous trajectories, we modified the iteration procedure slightly. Specifically, we minimized both the distance between the ensemble of measured reactive trajectories and the string point within each individual Voronoi cell, as well as the overall distance between each individual trajectory and the trial RC (Fig. S2, see Methods for details). The iteration procedure gives the RC of TGF-β induced EMT in A549 cells (Fig. 3A left, Fig. S3), which characterizes common features of the reactive trajectories (Fig. 3A right). Along the curve s the cell shape changes dominantly through elongation and growth ( Fig. 3B top), and most of the 13 vimentin Haralick features increase or decrease monotonically and continuously over time (Fig. 3B bottom).
For quantitative analysis, it is reasonable to assume that dynamics of the collective variables (x) can be described by a set of Langevin equations in the morphology/texture feature space, is a governing vector field, and η are white Gaussian noises with zero mean. Then in principle one can reconstruct F(x) from the single cell trajectory data, F(x) = dx / dt . Specifically we restricted to reconstructing the dynamics along the RC, and the ansatz becomes a 1-D convection-diffusion process, ds / dt = −dφ / ds +η . Notice that for a 1D system even without detailed balance one can define a quasi-scalar potential ϕ 16,17 . Figure 3C presents an enlarged view of an area in Fig. 3A, which shows the instant velocities (ds/dt) of various trajectories segments. Numerically we related the potential gradient with the mean velocity within the i-th Voronoi cell by (Fig. 3D), The sum was over all N i recorded trajectory segments that locate within the i-th Voronoi cell at any time t (Fig. 3A & C). On the obtained curve of dϕ/ds v.s. s (Fig. 3E), the zeroes correspond to stationary points of the potential. We then reconstructed the quasi-potential through  The quasi-potentials in Fig. 3G provide mechanistic insights on the TGF-β induced EMT. Before TGF-β induction, cells reside on the untreated cell potential centered with a potential well corresponding to the epithelial attractor. After induction, the system relaxes following the new potential into a new well corresponding to the mesenchymal attractor. Notably in the new potential the original epithelial attractor almost disappears, reflecting that the epithelial phenotype is destablized under the applied TGF-β concentration.
One might question whether by approximating a complex dynamical process such as EMT along the 1D RC one may miss some key dynamical information. To examine such possible limitation, we adopted a different strategy of dimension reduction over the original reactive trajectories using self-organizing map (SOM). SOM is an unsupervised artificial neural network that utilize neighborhood function to represent the topology structure of input data 18 . The algorithm clusters the recorded cell states into 144 discrete states (Fig. 4A), and represents the EMT process as a Markovian transition process among these states (Fig. 4B). Shortest path analysis using the Dijkstra algorithm 19 over the transition matrix reveals two groups of paths: vimentin Haralick PC1 varies first, and concerted variation of morphology PC1 and vimentin Haralick PC1, with finite probabilities of transition between the two groups (Fig. 4C). This result is consistent with our previous trajectory clustering analysis using soft dynamics time warping (DTW) distance between reactive trajectories 20 , and suggests that indeed the 1D approximation misses some dynamical features of the process. To further support this conclusion, we also examined the density of reactive trajectories, ρ R , in the plane of morphology PC1 and vimentin Haralick PC1. The contour map of ρ R shows two peaks corresponding of the two groups of shorted paths in the directed network (Fig. 4D). The peak that vimentin Haralick PC1 varies firstly is higher than the peak of concerted variation, indicating more reactive trajectories along this path. To go beyond the 1D RC formalism, we grouped the reactive trajectories based on the DTW distance, then identified the RCs for each group separately following the modified string method (Fig. 4E). The two RCs first diverge from the E region to follow two distinct paths, then converge within the M region. In one group (Fig. 4F), most of the Haralick feature changes take place before major morphology change. The two segments of the quasi-potential curve with shallow and steep slopes also reflect the two-stage dynamics. For the group with concerted dynamics (Fig. 4G), both cell shape and Haralick features vary gradually along the RC. The quasi-potential also decreases steadily, only interrupted by a short transition region. These analyses provide direct evidence that TGF-β induced EMT proceeds through parallel paths with distinct dynamical features.

DISCUSSIONS
In summary, through analyzing experimentally measured single cell trajectories with dynamical systems theories we obtained quantitative mechanistic insights of TGF-β induced EMT. The dynamics can be mapped to that of a chemical system transiting from the ground to an excited potential then relaxing to a new stable state. This work demonstrates the existence of a unified theoretical framework describing transition and relaxation dynamics in systems with and without detailed balance. CPTs, with relevant spatial and temporal scales orders of magnitude larger than typical molecular systems and thus more accessible for multiplex quantitative measurements, can serve as an ideal test ground for reaction rate theories.

1) Procedure for determining a reaction coordinate
We followed a procedure adapted from what used in the finite temperature string method for numerical searching of reaction coordinate and non-equilibrium umbrella sampling 14,15 , with a major difference that we used experimentally measured single cell trajectories (Fig. S2).
a) Identify the starting and ending points of the reaction path as the means of data points in the epithelial and mesenchymal regions, respectively. The two points are fixed in the remaining iterations. b) Construct an initial guess of the reaction path that connects the two ending points in the feature space through linear interpolation. Discretize the path with N (= 35) points (called images, and the k-th image denoted as ! with corresponding coordinate X(s k )) uniformly spaced in arc length. c) Collect all the reactive single cell trajectories that start from the epithelial region and end in the mesenchymal region.  Fig. 4. We used a larger smooth factor for the latter to avoid overfitting as there is less data in each of the two groups compared to the overall data set used in Fig. 3. e) We iterated the whole process in step 3 until there was no further change of Voronoi polyhedron assignments of the data points. f) For obtaining the quasi-potential of a larger range of s, extrapolate the obtained reaction path forward and backward by adding additional image points (5 for the single reaction path and 2 for the two parallel reaction paths) beyond the two ends of the path linearly, respectively. These new image points are also uniformly distributed along s as the old image sets do. Re-index the whole set of image points as {s 0 , s 1 ,..., s i ,..., s N , s N +1 }.

2) Reconstruction of quasi-potential along the reaction coordinate
Based on the theoretical framework in Fig. 3D, we followed the procedure below: (a) The N + 2 image points of an identified RC divide the space in N + 2 Voronoi cells that data points can assign to. Ignore the first and last Voronoi cells, and use the remaining N cells for the remaining analyses.

3) calculation of dynamics of morphology and Haralick features along reaction path
The reaction path is calculated in the principal component (PC) space of morphology PC1, vimentin Haralick PC1, PC3 and PC4. Distribution of cells show significant shift before and after TGF-β treatment in these dimensions 7 . To reconstruct dynamics in the original features space from PCs, the reaction path's coordinates on the other dimensions of PC are set as means of data in the corresponding Voronoi cell of each point on the reaction path. We obtain the reaction path in full dimension of PC space. The dynamics of morphology and Haralick features are calculated by inverse-transform of coordinates of PCs.

4) Self-organizing map and shortest transition paths in the directed network
The self-organizing map is an unsupervised machine learning method to represent the topology structure of date sets. We used a 12x12 grid (neurons) to perform space approximation of all reactive trajectories (Fig. 4A). The SOM was trained for 50 epochs on the data with Neupy (http://neupy.com/pages/home.html). We set the learning radius as 1 and standard deviation 1.
These neurons divide the data into 144 micro-clusters ( { } ψ ). With the single cell trajectory data, we counted the transition probabilities from cluster i to cluster j (including self), with p i, j j ∑ = 1.
If the transition probability is smaller than 0.01, the value is then reset as 0. With the transition probability matrix, we built a directed network of these 144 neurons (Fig. 4B). The distance (weight) of the edge between neuron i ψ and neuron j ψ is defined as the negative logarithm of transition probability ( , log i j p − ). We ignored the connection between two clusters. We set the neurons that close to the center of epithelial and mesenchymal state (sphere with radius = 0.7) as epithelial community and mesenchymal community, respectively, and used Dijkstra algorithm to find the shortest path between each pair of epithelial and mesenchymal neurons 19 with NetworkX 22 . We recorded the frequency of neurons and edge between these neurons that are past by these shortest paths.

5) Calculation of density of reactive trajectory
The density of reactive trajectory on the plane of morphology PC1 and vimentin Haralick PC1 is calculated with the following procedure: a) Divide the whole plane into 200×200 grids. b) In each grid, count the number of reactive trajectory (only the parts of each reactive trajectory that are in the intermediate region are taken into consideration) that enters and leaves it. If a reactive trajectory passes certain grid multiple times, only 1 is added in this grid's density. Thus, the density matrix is obtained.