Multispectral Fluorescence Microscopy Enables Accurate Single-Cell Tracking.

Compared with phase-contrast imaging of CLS structures, we found that multispectral fluorescence imaging clearly delineated individual cell boundaries, even in dense cell clusters (Fig. 1B and Movies S1, S2, and S3). We tracked four phenotypes that described cellular morphology: (i) area, (ii) perimeter, (iii) box ratio (the ratio of length and width of the smallest bounding rectangle of each cell, to indicate the extent to which the cell is elongated), and (iv) maximum radius (the maximum distance between the cell centroid and its boundary, to serve as a more sensitive measure of filopodia formation than box ratio). Cell morphology has been demonstrated to be a sensitive indicator for interactions by endothelial cells with the extracellular environment. For example, changes in the morphology of individual endothelial cells are dramatically more sensitive to hypoxia than to cell number or other common indicators of blood-vessel phenotype (22, 23). The specific parameters we used to measure cell morphology have also been used to characterize other cell types (e.g., area and box ratio for keratocytes) (21, 24).

We also tracked another three phenotypes that described cellular motility: (i) displacement (the distance traveled relative to the previous time point), (ii) accumulative distance (the sum of displacements since t = 0), and (iii) displacement vector angle change (the direction of displacement). Previous studies have established the importance of migration of endothelial cells in microvasculogenesis (25, 26).

Overall, we tracked 125, 183, and 165 cells (in three datasets) at 35 time points, and we analyzed images from 11 time points (starting at t = 0 and analyzing the image every 60 min afterward). This analysis produced a matrix of 36,421 array elements, where each element represented a quantifiable phenotype for one cell at one time point.

Microvasculogenesis Occurs in Five Distinct and Reproducible Phases According to Quantitative Analysis of Collective Multicellular Behavior.

We first studied how the population-averaged phenotypes progressed dynamically over time and whether the changes were statistically reproducible. Quantitative analysis of individual cellular behaviors revealed, at the population level, a carefully orchestrated sequence of subevents (Fig. 2 shows four key parameters in normalized units, or N.U.; Fig. S2 shows all eight parameters in absolute units):

• 1. Rearrangement and aggregation, where HUVECs migrated to spatial positions in advance of direct cell–cell contact. Cell displacement ranged from low to high. As cells migrated, they also spread, leading to increases in cell area and perimeter while maintaining relatively round shapes.

• 2. Spreading, where cells stayed in similar positions but spread. Displacement was small, whereas surface area, perimeter, and box ratio increased consistently.

• 3. Elongation and formation of cell–cell contacts, where cells elongated with cell area decreasing slightly. An increase in cell displacement as cells from adjacent nodes moved toward each other, coupled with a decrease in cell area that reflected rounding of cells, resulted in elongated protrusions as reflected by an increase in box ratio.

• 4. Plexus stabilization, with a decrease in cell displacement that reflected stabilization of the formed structure. Cell area and perimeter also decreased, consistent with an increase of thickness of cells (assuming no large changes in cell volumes); an increase in cell thickness was also confirmed by confocal imaging.

• 5. Plexus reorganization. After the initial structure of the plexus took shape, displacement exhibited a sharp increase to its maximum, whereas perimeter and area decreased gradually. These changes reflected a reorganization of polygonal structures in the vascular plexus and continued until the end of the assay. (Over extended periods of time, regression of the plexus is known to occur.)

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Dynamics of morphological and motility phenotypes during microvasculogenesis, based on averages across cell population. (A) Normalized population-averaged parameters for perimeter (red), box ratio (orange), area (purple), and displacement (green). Normalized units of 0 and 1 corresponded to the lowest and highest values, respectively, attained by each parameter over all time points. The SE is indicated for each data point (for n = 90, 112, and 132 cells that stayed in the field of view for the entire assay, respectively, for the three datasets). (B) Averages of population-averaged phenotypes across three datasets, with the five distinct temporal phases (as identified by visual inspection of changes in the slopes of phenotypes) shaded in white and gray. (C) Schematic model of CLS formation. Red arrows show the distinct feature of the population behavior at each phase (e.g., aggregation in first phase). Light blue cells show a subpopulation of cells that form the nodes. Previously unidentified phases for microvasculogenesis are indicated by asterisks.

These population-level trends were quantitatively reproducible among the three independent trials (see Fig. S2D and Movie S4, showing images over all five phases; the minimum number of cells required to represent the population-level behavior is shown in Fig. S3).

Behaviors of Most Endothelial Cells Follow a Handful of Unique Patterns.

Although useful, such dynamic bulk averages (reflecting the evolution of the multicellular population) can arise from a myriad of underlying distributions of individual cellular behaviors (4). We examined next whether there existed subgroups of endothelial cells that exhibited similar behavioral dynamics and how closely correlated such dynamics are to the population averages. We used a correlation clustering algorithm, which determines membership in a cluster based on pairwise comparisons (27); here, cells can be represented as vertices on a graph, with correlations (i.e., if Pearson's correlation coefficient of the phenotype values at different time points for the first cell vs. those for the second cell exceeded 0.6) represented as edges (Fig. 3A and Fig. S4A). We identified multiple dominant patterns of dynamic behavior for each phenotype, which were reproducible across different independent datasets (Fig. S5). We could group ~60–70% of cells into only three distinct clusters of dynamic behavior with respect to cell area (Fig. 3B). Similarly, analysis of three other phenotypes (displacement, box ratio, and perimeter) also grouped 60–70% of cells into three major clusters of phenotypic behavior (see Fig. 4A).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Clustering of cells with similar temporal variations in phenotypes. (A) Correlation clustering to identify cells with similar temporal dynamics of phenotypes. We performed this algorithm on the dataset with the highest number of completely tracked cells (132 cells). The two graphs (Left) show how, for cell area, correlations of the temporal dynamics for two pairs of cells (i.e., cell areas at different time points) were calculated; a correlation was found if the Pearson correlation coefficient was above 0.6. (Right) A whole-system graph G(V, E) for area, where each cell is represented as a vertex, V, and the presence of correlation between two cells shown as an edge, E. (B) Graphical representations for six identified clusters (with the corresponding number of cells) according to cell area. (C Upper) Segmented time-lapse images with cells false-colored according to their cell-area cluster (with three clusters of cells shown in purple, light blue, and dark blue and unclustered cells shown in white). (Lower) The graph shows the temporal progression of the cell areas for the three most abundant clusters, peaking at 60, 120, and 300 min, respectively.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Correlation of migratory dynamics to structural role in the final plexus. (A) Temporal variations in displacement for the three most abundantly populated cell clusters, relative to the population-averaged behavior. Colored lines show the dynamics of all individual cells that contribute to the cluster average (thick colored lines, with the most to least abundant clusters indicated by light blue, dark blue, and purple). (B) Temporal variations of the absolute values of displacement. Population average is shown as a thick black line, and the averages of the three most populated clusters are shown as thick colored lines. (C) Image of plexus at t = 540 min, with cells at nodes and branches indicated. Black lines indicate branches, and red circles indicate the nodes. (D) Percentage of cells in each cluster forming node cells in the final plexus. (E) In a 50:50 mixture of untreated cells and cells treated with myosin II inhibitors, the second most populated cell cluster (Left, of which 70% were drug-treated cells) showed maximum displacement at t = 240 min, with most of these cells forming nodes. The third most populated cluster (Right, of which 75% were untreated cells) showed maximum displacement at t = 0 min, with most of these cells forming branch structures.

We could also identify four “classes” of morphological behaviors, where each class exhibited the same dynamic behaviors in both degree of cell spreading (cell area) and branching morphology (box ratio). The four classes represented the behaviors of more than half (51%) of the cell population (Fig. S6B). Hence, the dynamic behaviors of most endothelial cells followed a small number of distinct patterns.

Cells in Superficially Similar Microenvironments Diverge in Behaviors and Fates.

To interpret the significance of such dynamic patterns, we tracked the behaviors of these cell clusters in Movies S1, S2, S3, S4, and S5 over the time course of tissue morphogenesis. Through segmentation of images, where cells were false-colored according to their cluster membership (Fig. 3C), the images showed how individual HUVECs, while appearing to reside initially in superficially similar microenvironments, diverged into distinct temporal patterns of cellular behaviors (Fig. 3C, Fig. S6A, and Movie S5). The first and third most populated cell clusters, for example, all spread quickly but then shrank in size; ultimately, in the final plexus structure, they assumed the structural role of connection points (nodes) for multiple CLS (Fig. 3C). By contrast, cells in the second most populated cluster, which were generally the midrange-size cells in their immediate vicinity, typically formed the branches in the final plexus structure.

We further characterized the relationship between the behavioral dynamics of each cell and their structural role in the CLS. We classified all endothelial cells in the final plexus structure as either nodes (38% of cells) or branches s2% of cells) (Fig. 4C). The strongest correlation was observed when clustered by changes in displacement (Fig. 4A). Here, 67% of the cells in the first cluster (whose displacements peaked at around t = 240 min) formed nodes (Fig. 4D); hence, the migratory dynamics of subpopulations of cells could be correlated to their structural role in the final plexus.

Initial Inhibition of Myosin II Changed Subsequent Displacement Dynamics and Fates.

In addition to the internal cytoskeletal machinery (28), endothelial cells can make use of other motility mechanisms, such as coordinated movement through neighboring cells (29, 30). We investigated whether individual cells can be biased toward one set of behavioral dynamics vs. another set by initially limiting the role of active cytoskeletal machinery. We mixed together a population of untreated endothelial cells with endothelial cells whose ability to form stress fibers and/or focal adhesions was initially disrupted (via direct inhibition of active myosin II via blebbistatin as well as inhibition of Rho-associated kinase (ROCK) via Y-27632; Fig. S7). As observed previously (31), the drug-treated cells could still change cell shape by forming filopodia and lamellipodia. Like the native population (Fig. 4A), a cell cluster emerged with a displacement curve peak at around t = 240 min (Fig. 4E Left). Interestingly, in this 50:50 mixture of drug-treated and untreated cells, most cells in this cluster were drug-treated, and most cells with these migratory dynamics (~80% in this case) again formed nodes in the final plexus, just like the native population. We note that cell–cell connections were beginning to be formed already when these cells began migrating (i.e., displacement is increasing during the “elongation” phase identified from Fig. 2); hence, coordinated movement was an available option for these cells. By contrast, another highly populated cell cluster that migrated quickly in the initial minutes consisted mainly of untreated cells and subsequently formed mostly branch structures.

Cellular Behavior Can Exhibit Large Variations with Different Types of Distributions.

The variations observed for some phenotypes (which ranged from 20% to 120% relative to mean values; Fig. 5A) are higher than those observed in other biological processes, such as gene expression [which exhibited up to ~60% variation relative to average population behavior (32)]. This comparison possibly reflects the complexity of cell–cell interactions and the diversity of phenotypic behaviors needed to form a multicellular tissue. The extent of variations also fluctuated over time: over the first phase of microvasculogenesis (aggregation), variations in displacement were large but rapidly decreased because most HUVECs settled into position for cell–cell networking. For other phenotypes (such as box ratio and perimeter), initial fluctuations among the cell population were relatively small but increased over the ~60 min of the aggregation phase.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Variations in behavior from cell to cell. (A) Normalized SDs across the cellular population for perimeter (red), box ratio (orange), area (purple), and displacement (green) at a given time point. The plot shows the degree of variability of each phenotype across the cellular population at different time points. (B) Normalized SDs of cell perimeter (red), box ratio (orange), area (purple), and displacement (green) for a given cell across all time points. The plot shows how much a given cell varies during the process of microvasculogenesis. For comparison, square brackets indicate the range of variability for a cell population at a fixed time (for the indicated parameter as obtained from A). (C) Histograms of percentage of cells vs. value of phenotype (cell area and box ratio) at three different time points. Although both phenotypes were derived from cellular morphology, they exhibited different statistical distributions throughout the time course. (D) Percentage of time points within indicated values of cell area or box ratio for three different cells that represented a range of distributions observed. A single cell can exhibit very similar behavior over all time points (e.g., cell no. 104, which exhibited a small cell area) or adopt two sets of very different behavior (e.g., cell no. 126, which was quite round or very elongated).

Different underlying distributions of phenotypes can give rise to the observed variability (Fig. S8). We plotted how two different phenotypes describing cell morphology were distributed across the population. For example, histograms enumerating cells with different cell areas showed that the endothelial cells tended to spread over time and then retract, maintaining either relatively unimodal or uniform distributions at all times (Fig. 5C Left). By contrast, histograms of box ratio values showed initially Poisson-type distribution, then relatively uniform distribution at t = 120 min, and finally bimodal distribution at the end of the assay (Fig. 5C Right). Although both phenotypes were derived from cellular morphology, they exhibited different statistical distributions throughout the time course.

Studies on other cell types have indicated that individual cellular variability over time is often less than the variability in the population (24). To examine how much a given cell varies during the process of microvasculogenesis, we examined the variability of individual cell areas across all time points (Fig. 5B). We found that, for cell area, most cells varied less over time than the population variability. (We obtained similar results for perimeter and displacement.) We also observed that some cells tend to exhibit very similar behavior over all time points (e.g., cell no. 104, which exhibited a small cell area), others vary within a range, whereas others adopt entirely different sets of behavior (e.g., cell no. 126, which was quite round or very elongated), depending on the stage of microvasculogenesis (Fig. 5D). These trends reflect the different types of variability that can take place because some cells exhibited changes within a range of bounded values, whereas others transitioned among different states altogether.

Finally, we note that the trends of the population averages were dominated by the behavior of only ~30–40% of cells and masked the behaviors of most cells in the population (displacement is shown in Fig. 4 A and B; other parameters are shown and discussed in SI Materials and Methods and Figs. S4 and S9).