Endothelial cells (ECs) are heterogeneous across and within tissues, reflecting distinct, specialised functions. EC heterogeneity has been proposed to underpin EC plasticity independently from vessel microenvironments. However, heterogeneity driven by contact-dependent or short-range cell–cell crosstalk cannot be evaluated with single cell transcriptomic approaches, as spatial and contextual information is lost. Nonetheless, quantification of EC heterogeneity and understanding of its molecular drivers is key to developing novel therapeutics for cancer, cardiovascular diseases and for revascularisation in regenerative medicine. Here, we developed an EC profiling tool (ECPT) to examine individual cells within intact monolayers. We used ECPT to characterise different phenotypes in arterial, venous and microvascular EC populations. In line with other studies, we measured heterogeneity in terms of cell cycle, proliferation, and junction organisation. ECPT uncovered a previously under-appreciated single-cell heterogeneity in NOTCH activation. We correlated cell proliferation with different NOTCH activation states at the single-cell and population levels. The positional and relational information extracted with our novel approach is key to elucidating the molecular mechanisms underpinning EC heterogeneity.

Endothelial cells (ECs) form the inner layer of blood and lymphatic vessels, and play major roles in tissue development, angiogenesis, inflammation and immune cell trafficking (Potente et al., 2011). ECs are functionally plastic and rapidly adapt to changes in the environment to preserve homeostasis. Local EC dysregulation is a hallmark of diseases such as atherosclerosis, ischemia and cancer (Park-Windhol and D'Amore, 2016). ECs from different tissues and vascular beds (e.g. arteries, capillaries and veins) exhibit distinct metabolism, morphology and gene expression (Augustin and Koh, 2017), and contribute in diverse ways to tissue development and regeneration (Itkin et al., 2016; Kusumbe et al., 2014). It is well established that ECs are phenotypically heterogenous not only among different tissues, reflecting specialised organ-specific functions (Rafii et al., 2016), but also within the same tissue. Maintenance of endothelial homeostasis depends on new ECs substituting senescent cells, and the role of endothelial progenitor cells with high repopulating potential has been highlighted in endothelia of large vessels (Yoder, 2018).

Inter-endothelial adherens junctions (IEJs) are dynamically regulated by VE cadherin (CDH5) shuffling between the cell membrane and intracellular compartments. This process presents variations across vascular beds (Augustin and Koh, 2017) and involves molecular mechanisms including VEGF receptors, cytoskeletal proteins and NOTCH family members. VE cadherin and NOTCH signalling are well studied in angiogenesis and development. Nonetheless, the role of these molecular mechanisms in EC monolayer maintenance is less clear and this knowledge is essential to understand vessel homeostasis in different organs in the human body. So far, a comprehensive analysis of EC cultures exploring and quantifying phenotypic variance has proven prohibitively difficult because of lack of adequate tools.

Single-cell phenotyping has identified and characterised intermediate cell states (MacLean et al., 2018; Siu et al., 2020), and demonstrated correspondence between phenotypes and function (Dueck et al., 2016). However, challenges in discriminating functional phenotypic variance from biological noise have emerged (Eling et al., 2019). Single-cell transcriptomic (sc-OMICS) data is becoming available (mostly in mouse) (Kalucka et al., 2020) and could, in principle, advance the characterisation of human ECs (Tikhonova et al., 2019) and provide an overview of molecular processes in distinct EC populations. Nonetheless, sc-OMICS data lacks spatial information, which is essential to map variable phenotypes to function and is also required to understand higher-level functions depending on cell–cell connectivity. Of note, the EC cell cycle and their rapid adaptative phenomena, such as rapid increase of endothelial permeability upon VEGF stimulation are regulated by built-in sensing mechanisms that depend on cell–cell interaction (Acar et al., 2008). Importantly, EC monolayers maintain their integrity over years while exerting a variety of system-level functions, which are emerging properties of cells in contact (McCarron et al., 2017). It has been proposed that endothelial adaptability and diversity of functions within an EC population depends on cell heterogeneity (McCarron et al., 2019).

To examine individual cell heterogeneity and extract spatial information from EC monolayers, we developed an endothelial cell profiling tool (ECPT) based on high-content image analysis (HCA). ECPT captures a wealth of cellular, subcellular and contextual information enabling extensive characterisation of the cell cycle and IEJs. This unbiased approach allows quantification of EC heterogeneity measuring feature variance (Eling et al., 2019). Taking advantage of machine learning based methods, we performed automated and accurate classification of IEJs using junctional CDH5 immunostaining, and evaluation of NOTCH activation at the single cell level. In total, we analysed data from >20,000 images across nine independent experiments detailing selected measurements from >300,000 cells including individual junction objects. Overall, we present (1) a novel tool for single EC profiling at a previously inaccessible scale, (2) the validity of this analysis to quantify previous observations and (3) new key relationships across features that are distinct between different EC types.

Challenging quantification of endothelial cell heterogeneity by high-content imaging

Although heterogeneity of ECs has been proposed, few transcriptomic studies have successfully quantified the phenotypic variance of ECs within the same vascular bed (Chavkin and Hirschi, 2020). We previously reported the value of specific vascular endothelial (VE) cadherin (CDH5) and NOTCH quantification to benchmark ECs in intact monolayers (Wiseman et al., 2019) and now extend our observations to include relational information between individual cells. These cannot currently be retrieved with single cell transcriptomic analysis. We cultured primary human umbilical vein ECs (HUVECs), human aortic ECs (HAoECs) and human pulmonary microvascular ECs (HPMECs) for 48 or 96 h in the absence or presence of VEGF. We aimed for sufficient cells to reach confluency upon adhesion and spreading (∼24 h), and we cultured for a further 24 or 72 h to enable formation of stable IEJs. ECs demonstrate a uniform cobblestone morphology under low magnification, confirming monolayer stabilisation.

Immunostaining and live imaging of EC monolayers has highlighted differences in cell morphology and junction patterns in response to stimuli or under distinct cell culture conditions (Abu Taha et al., 2014; Lampugnani et al., 1992; Seebach et al., 2016; Vestweber et al., 2009). We focussed on analysing proliferation, junctions, and NOTCH activation. We stained for CDH5, NOTCH1 (or the actin cytoskeleton), HES1 (a NOTCH target gene) and nuclei. In line with previous literature, qualitative inspection at high magnification revealed remarkable morphological differences between the cell types in culture (Fig. 1).

Fig. 1.

Overview of EC heterogeneity. Microphotographs comparing HUVECs, HAoECs and HPMECs stained for DNA (Hoechst 33342, nuclei), VE-cadherin (CDH5), NOTCH1 and HES. Images representative of nine experiments. Scale bar: 50 µm.

Fig. 1.

Overview of EC heterogeneity. Microphotographs comparing HUVECs, HAoECs and HPMECs stained for DNA (Hoechst 33342, nuclei), VE-cadherin (CDH5), NOTCH1 and HES. Images representative of nine experiments. Scale bar: 50 µm.

EC junctions have been linked to specific function and phenotypes (Abu Taha et al., 2014; Lampugnani et al., 1995). Linear and stabilised junctions are present in the mature quiescent phenotype; jagged and discontinuous junctions results from proliferation, migration or the immature/mesenchymal phenotype; reticular patterns, observed in more stabilised junctions, can be also observed as transient structures (Kim and Cooper, 2018) and are associated with immune cell transmigration (Fernández-Martín et al., 2012).

In our study, HUVECs, HAoECs and HPMECs showed different junction patterns and cell morphology in standard in vitro culture conditions. Inspection of NOTCH1 and HES1 staining intensity and localisation at the single-cell level appeared to be highly heterogeneous both intra-population and inter-population (Fig. 1). Measuring and scoring individual EC phenotypes to demonstrate, quantify and explain EC heterogeneity at single cell resolution is a prohibitively tedious and biased approach if not automated. Therefore, we set out to create a dedicated platform that can seamlessly and comprehensively quantify variance of phenotypic parameters in cultured EC monolayers.

Creation of a semi-automated image analysis pipeline for EC phenotyping

We developed a method to profile EC phenotypes at the single-cell level based on extracted high-content analysis (HCA) features: Endothelial Cell Profiling Tool (ECPT, Fig. 2 and Appendix 1 at https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers). We used open-source software to generate an image analysis workflow able to take input from fixed or live images, with the aim of characterizing EC heterogeneity at the single-cell level (Fig. 2). Features and benchmark comparisons between ECPT and available tools is presented in Table 1. Fig. 2A demonstrates the unique capability of ECPT to precisely segment cells according to junctional CDH5 staining. Accurately segmenting cells is key to all downstream analyses and to evaluate cell junctions. However, CDH5 staining is very heterogeneous across different cell populations including those of different origins (e.g. HAoECs, HPMECs and HUVECs) and treatment (e.g. untreated versus VEGF-treated ECs). In terms of image analysis, this renders standard thresholding techniques insufficient to appropriately contrast large collections of images for segmentation. To overcome this problem, we developed a workflow in Fiji/ImageJ leveraging the Weka segmentation, a powerful machine learning (ML) tool. Fig. 2Ai shows the original CDH5 image, Fig. 2Aii shows the corresponding output images following application of a Weka model based on an annotated dataset of 70 CDH5 images chosen randomly across our database of >20,000 images. Fig. 2Aiii shows an overlay of individual junction objects. Individual junctions between two cells are precisely identified via ECPT and junction features are measured. Different junction classes were obtained by applying ML-aided object classification (Cell Profiler Analyst, CPA using Fast Gentle Boosting) based on a the measurements of a junction and colour code in Fig. 2Aiii. Nuclei segmentation and downstream measurements were performed using standard Cell Profiler (CP) modules and thresholding methods (Otsu, Fig. 2Bi) to evaluate nuclear morphology, cell cycle (Fig. 2Bii) and NOTCH signalling (Fig. 2C).

Fig. 2.

Image segmentation and ECPT features. (A) Acquisition of high resolution (40× OM) images on each channel. The CDH5 channel (i) is used to segment the cells and outline the cell junctions using Weka segmentator Fiji (ii), allowing individual junction and cells segmentation (iii). (B) Nuclei are stained with Hoechst 33342 (i) and the intensity is used to segment nuclei and measure DNA intensity to perform cell cycle analysis (ii). (C) Measurements of HES1 nuclear signal intensity (i). Empirical cumulative distribution functions (ECDF) across all images by cell type (ii). Example nHES1 intensity map (iii) and frequency of cells with negative nLMi relative to HES1. (D) Single cells PCA analysis on untreated EC cultured for 96 h. Arrows indicate loading of different variables. Colours in C and D represent EC type. Scale bars: 50 µm.

Fig. 2.

Image segmentation and ECPT features. (A) Acquisition of high resolution (40× OM) images on each channel. The CDH5 channel (i) is used to segment the cells and outline the cell junctions using Weka segmentator Fiji (ii), allowing individual junction and cells segmentation (iii). (B) Nuclei are stained with Hoechst 33342 (i) and the intensity is used to segment nuclei and measure DNA intensity to perform cell cycle analysis (ii). (C) Measurements of HES1 nuclear signal intensity (i). Empirical cumulative distribution functions (ECDF) across all images by cell type (ii). Example nHES1 intensity map (iii) and frequency of cells with negative nLMi relative to HES1. (D) Single cells PCA analysis on untreated EC cultured for 96 h. Arrows indicate loading of different variables. Colours in C and D represent EC type. Scale bars: 50 µm.

Table 1.

ECPT features and comparison with other available software

ECPT features and comparison with other available software
ECPT features and comparison with other available software

Overall, extending our previous proof-of-principle study (Wiseman et al., 2019), we chose features describing cell proliferation, morphology, spatial organisation, NOTCH activation and ML-aided classification of junctions. We then applied the ECPT workflow to a dataset composed of >20,000 images from nine independent experiments obtaining a bulk of >300,000 single cells across different cell types (i.e. HAoECs, HPMECs and HUVECs) and experimental conditions (i.e. initial cell density, time in culture and VEGF treatment).

We first asked whether HAoECs, HPMECs and HUVECs could be distinguished solely using the chosen features. We set to capture variance within our multivariate datasets, to identify discrete populations among heterogeneous cells and to understand, which parameters had greater weight in defining cell phenotype. Principal components analysis (PCA) based dimensionality reduction cumulatively captured 51.1% of the variance (Fig. 2D). HAoECs, HPMECs and HUVECs did not form distinct clusters and the cell populations were partially overlapping. Nevertheless, cell populations segregated according to the first two PCA components. This indicates that our parameters of choice capture key differences in the three EC populations. To understand the origin of this heterogeneity we dissected individual sets of features dependent on the same biological mechanism. We developed a dedicated user-friendly, open-source interface using the shiny package within Rstudio (https://shiny.rstudio.com) inspired by guidelines described previously (Lord et al., 2020). This shiny application allows subsetting through interactive and iterative selection of single cells or groups of cells for further comparison and analysis within R studio (Appendix 1 at https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers). We welcome readers to further explore our dataset available at https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers. To understand certain phenotypic differences in detail, we quantified variance of a few selected traits: proliferation, junction stability and NOTCH activation.

Single-cell analysis of cell cycle

The integrated Hoechst 33342 nuclei signal intensity, commonly used as a proxy for DNA content in flow cytometry, was analysed based on previously published protocol (Roukos et al., 2015) to evaluate cell proliferation. The different stages of cell cycle (G1, S and G2/M) are thus defined via thresholding, with distinct peaks for G1 and G2/M (Fig. 3A). Late mitotic (LM) cells were detected by ML-aided object classification using features of Hoechst (DNA) and CDH5 stainings, including texture features (Fig. 3A). This automated method allowed us to overcome specific limitations of previous approaches (Jones et al., 2009; Roukos et al., 2015).

Fig. 3.

Cell cycle analysis. (A) Histogram plots of DNA intensities by cell type (HUVECs, HAoECs and HPMECds) and treatments [control (CTRL) and VEGF]. Cells are classified in cell cycle phases depending on the signal intensity. Inset displays LM cells engaged in mitosis which were detected and quantified using CPA. (B) Number of nuclei per image (cell density) in sub-confluent (black) or confluent (red) experiments for the three different EC types (HUVECs, HAoECs and HPMECs) and treatments (CTRL and VEGF). (C) Percentage of dividing cells (S, G2/M, LM) per image in HAoEC, HPMEC and HUVEC with or without treatment with 50 ng/ml VEGF in subconfluent or confluent conditions. The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each statistical comparison indicated as annotations in individual plots. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 3.

Cell cycle analysis. (A) Histogram plots of DNA intensities by cell type (HUVECs, HAoECs and HPMECds) and treatments [control (CTRL) and VEGF]. Cells are classified in cell cycle phases depending on the signal intensity. Inset displays LM cells engaged in mitosis which were detected and quantified using CPA. (B) Number of nuclei per image (cell density) in sub-confluent (black) or confluent (red) experiments for the three different EC types (HUVECs, HAoECs and HPMECs) and treatments (CTRL and VEGF). (C) Percentage of dividing cells (S, G2/M, LM) per image in HAoEC, HPMEC and HUVEC with or without treatment with 50 ng/ml VEGF in subconfluent or confluent conditions. The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each statistical comparison indicated as annotations in individual plots. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Cell density is a well-established determinant of EC proliferation (Andriopoulou et al., 1999; Bazzoni and Dejana, 2004); therefore, we first set out to establish whether cell densities were homogeneous across cell types and different experiments.

By subculturing HUVECs, HAoECs and HPMECs under the same standardised conditions, we noticed that these EC types present different intrinsic proliferation rates. We thus quantified proliferating cells and compared subconfluent and confluent cells, which were also cultured for 48 h or 96 h upon seeding, as a longer time in culture has been shown to further stabilise the EC monolayer (Bazzoni and Dejana, 2004).

Proliferation was higher in subconfluent HAoECs and in HUVECs at baseline (average 16.9% and 13.8% of dividing cells per field) and lower in HPMECs (average 7.8%). VEGF induced a small detectable increase in proliferation in all subconfluent ECs (18.9%, 9.6% and 16.2% in HAoECs, HPMECs and HUVECs, respectively). In subconfluent cultures, an average of 3.5% to 6.7% of cells were LM, demonstrating cells were effectively dividing at the time of the experiment.

In confluent cultures, the proliferation was lower in all ECs except HPMECs, which maintained values as in subconfluent conditions. In confluent cultures, the percentage of cells in cell cycle ranged from 9.8% to 12.7% in all ECs with little differences between control and VEGF-treated conditions. Longer culture conditions also had little effect on proliferation rates. In all confluent cultures examined, the average percentage of LM cells dropped to values between 1.2% and 2%. Overall, our results are consistent with previous observations and demonstrate that EC proliferation is affected by cell density. Furthermore, confluent cultures under our conditions have low proliferation rates, consistent with that observed in vivo (∼4.8% of cells in S, G2 or M phases as measured by scRNA-seq of mouse aortic endothelium; Lukowski et al., 2019).

Together these data validate the use of ECPT to characterise EC cell cycle under different experimental conditions. These are key quality controls when analysing cell phenotype and dynamic structures, such as inter-EC junctions.

Analysis of morphology and cytoskeleton arrangement

Cell morphology statistics were obtained directly from cell objects segmentation. Fig. 4A shows cell area comparisons across EC types and cell cycle stage in either subconfluent (red, green and blue traces) or confluent (dashed grey traces) conditions. Overall, the three cell types had different average area and, as expected in subconfluent cultures, cell area increased with phases of cell cycle (between G0/G1 and G2/M phases) (Ginzberg et al., 2018). In confluent cultures, cell area did not change across phases of cell cycle and confluent cells had lower proliferation rates. VEGF treatment did not induce significant differences in cell area in these conditions.

Fig. 4.

Analysis of cell morphology and cytoskeleton. (A) Scaled density plots showing distribution of cell areas across cell types (HAoECs, HPMECs and HUVECs) and treatment [control (CTRL) and VEGF] in sub-confluent (coloured traces, colour code as shown in legend) or confluent (dashed grey traces) conditions. (B) Width-to-length ratio (WLR) of cell areas across cell types (HAoECs, HPMECs and HUVECs) and treatment (CTRL and VEGF). (C) Analysis of stress fibres across cell types (HAoECs, HPMECs and HUVECs) and treatment (CTRL and VEGF). The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each representation or statistical comparison indicated as annotations in individual plots. ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 4.

Analysis of cell morphology and cytoskeleton. (A) Scaled density plots showing distribution of cell areas across cell types (HAoECs, HPMECs and HUVECs) and treatment [control (CTRL) and VEGF] in sub-confluent (coloured traces, colour code as shown in legend) or confluent (dashed grey traces) conditions. (B) Width-to-length ratio (WLR) of cell areas across cell types (HAoECs, HPMECs and HUVECs) and treatment (CTRL and VEGF). (C) Analysis of stress fibres across cell types (HAoECs, HPMECs and HUVECs) and treatment (CTRL and VEGF). The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each representation or statistical comparison indicated as annotations in individual plots. ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Analysis of the width-to-length ratio (WLR) in confluent cultures (Fig. 4B) shows a clear reduction in WLR in VEGF-treated cells, demonstrating cell elongation and confirming effectiveness of VEGF treatment. No difference was noted in WLR across different EC types or cell densities. By qualitative inspection of phalloidin-stained cells, we noticed that different ECs had qualitatively different distribution of actin cytoskeleton (see Abu Taha et al., 2014; Millán et al., 2010). In our subconfluent conditions, HAoECs, HPMECs and HUVECs demonstrated different frequencies of cells with stress fibres (Fig. 4C). Data analysis showed that HAoECs and HUVECs had the highest proportion of cells with stress fibres, whereas HPMEC had the lowest frequency. VEGF treatment induced a significant increase of cells with stress fibres for HAoECs but a decrease in HUVECs, whereas HPMECs maintained a low level - like baseline.

Overall, the EC types analysed were remarkably different in terms of cell area and arrangement of cytoskeleton. As expected, cell area increased with progression of cell cycle in subconfluent cultures, but this effect was less prominent in confluent conditions where cells had lower, yet appreciable, proliferation rates. As expected, VEGF treatment clearly induced elongation in all cell types and induced cytoskeleton rearrangement in HAoECs and HUVECs.

Analysis of junction heterogeneity

The description of EC junction dynamics and morphology has been of great interest due to the link with vascular permeability and angiogenesis. In general, non-proliferating ECs are also those with continuous (‘stabilised’) IEJs, while migrating or proliferating EC demonstrate jagged or discontinuous junctions, which are rapidly remodelling (Fernández-Martín et al., 2012; Millán et al., 2010). Average population measurements have been used by us (Veschini et al., 2011, 2007; Wiseman et al., 2019) and others to assess stability of EC monolayers, which in turn correlate with functions such as trans-endothelial permeability to large proteins (Ferrero et al., 2004). Assessing IEJs at the single-cell level enables subtle nuances in inter-endothelial cells connectivity within the same monolayer to be highlighted. Nonetheless, measuring and classifying junctional signal by image analysis is technically challenging, and no currently available software can evaluate multiple images in a standardised way.

To overcome this obstacle, we automated analysis of CDH5 pattern and junction morphology using the ML capabilities of CPA and integrated this into ECPT. Previous studies on EC junction (Brezovjakova et al., 2019; Seebach et al., 2015; Wiseman et al., 2019) did not analyse the proportion of different junction pattern; our method classifies objects from whole junctions using an expert-trained ML algorithm (CPA with Fast Gentle Boosting, Appendix 1 at https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers). In line with the literature in the field (Seebach et al., 2016), we classified junctions in a scale from 0 to 5 (Fig. 5A). J0, J1 and J2 are discontinuous, highly jagged and jagged junctions, respectively. J3 and J4 are linear, with J4 having a continuous CDH5 signal distributed over a larger area than J3. J5 junctions are visually reminiscent of the reticular junctions previously described (Millán et al., 2010) but could also appear as transient structures (Kim and Cooper, 2018; Seebach et al., 2020). J0 junctions might result from mis-segmented cells; therefore in downstream analyses we considered cells with less than 20% J0 junction type.

Fig. 5.

Image-based EC junction analysis. (A) Examples of junctional structures imaged in EC monolayers. The junctions are classified as J0 (not a junction or highly discontinuous), J1 (highly jagged and discontinuous), J2 (jagged), J3 (linear), J4 (linear reinforced signal) and J5 (linear/reticular). (B) Percentage of cells in the different categories for across cell types under control (CTRL) conditions and cultured for 96 h. (C) Repartition of the cells in the different classes after 48 h or 96 h in culture. (D) Percentage of cells in the different classes in basal conditions or after VEGF treatment.. The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each statistical comparison indicated as annotations in individual plots. *P<0.05, **P<0.01, ***P<0.001, &P<0.001 against HPMECs, #P<0.001 against HUVECs (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 5.

Image-based EC junction analysis. (A) Examples of junctional structures imaged in EC monolayers. The junctions are classified as J0 (not a junction or highly discontinuous), J1 (highly jagged and discontinuous), J2 (jagged), J3 (linear), J4 (linear reinforced signal) and J5 (linear/reticular). (B) Percentage of cells in the different categories for across cell types under control (CTRL) conditions and cultured for 96 h. (C) Repartition of the cells in the different classes after 48 h or 96 h in culture. (D) Percentage of cells in the different classes in basal conditions or after VEGF treatment.. The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each statistical comparison indicated as annotations in individual plots. *P<0.05, **P<0.01, ***P<0.001, &P<0.001 against HPMECs, #P<0.001 against HUVECs (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 5B shows the average percentage of junctions of each type across the three EC lines examined. HAoECs had a prevalence of the J1, J2 and J3 types with scarce J4. HPMECs were characterised by a high proportion of J4 junctions, whereas HUVECs had a high proportion of the J2 type. Overall, HPMECs and HUVECs had opposing phenotypes, whereas HAoECs were more heterogeneous. Time in culture upon seeding has been shown to affect junction response to perturbation (Andriopoulou et al., 1999), therefore we set to evaluate this effect under our experimental conditions. Fig. 5C shows the average percentage of junctions of each type per cell per field comparing cells cultured for either 48 h or 96 h. The effect of longer culture in all cells was an increase in J4 and J5 types, with a corresponding decrease in J3. Furthermore, a decrease in the highly jagged J1 type with a proportional increase in the less jagged J2 type was also noted. Overall, we conclude that longer time in culture significantly and positively affects junction quality in line with previous reports (Andriopoulou et al., 1999; Bazzoni and Dejana, 2004). Finally, we set to evaluate the effect of VEGF treatment in our system. As shown in Fig. 5D, the main effect of VEGF treatment on confluent cells cultured for 96 h upon seeding was a sharp increase in J3 with a corresponding decrease in all the other junction types.

Altogether, these data demonstrate that the EC lines are very heterogeneous in respect to IEJs, and that many cells have a composition of different junctions, which might also reflect differential connectivity with different neighbours. This junctional heterogeneity has been previously described and can now be detected and quantified by our ECPT. This suggests that IEJ architecture can be regulated locally as all cells were cultured under identical experimental conditions and all individual ECs within monolayers exposed to the same environment. Overall, we showed that image analysis and supervised ML can be used to characterise EC junctions and highlight differences at the population level. Our workflow was able to efficiently distinguish changes in junction morphology after different time in culture or VEGF treatment and the junctional status of single cells can be used to study intra-population heterogeneity.

Analysis of NOTCH activation

NOTCH signalling is a key modulator of EC development and function, but assessment of signal heterogeneity in EC monolayers and its potential role in regulating these functions is currently lacking.

To assess heterogeneity in NOTCH signalling in our experiments we measured the intensities of intra-nuclear NOTCH1 (nN1) and intra-nuclear HES1 (nHES1) by ECPT. Fig. 6A illustrates the density distribution of the two signal intensities for all cells in our dataset at baseline, demonstrating bulk differences between EC types. Analysis of mean signal intensity by microscopy field did not highlight striking differences across EC types and treatment. However, inspection of density distributions highlighted differences, and suggests that these might originate in different repartitions of cells across signal intensities. To investigate this aspect, we binned signal intensities as illustrated by the banding in Fig. 6A and calculated summary statistics by microscopic field for these bins. The results of these analyses for nN1 and nHES1 are shown in Fig. 6B,C. Analysis of nN1 highlighted differences in proportion of cells with either intermediate (28) or high (210) signal. At baseline, HAoECs had the lowest percentage of cells in the intermediate intensity category and the highest percentage of cells in the high category. The opposite was true for HUVECs, whereas HPMECs displayed a similar percentage in the two categories and were intermediate between the other EC types. VEGF treatment set the moderate cells at intermediate levels for all cell types, while it increased the percentage of high cells in HUVECs.

Fig. 6.

Analysis of NOTCH activation. (A) Scaled density distributions of nN1 and nHES1 for HAoECs, HPMECs and HUVECs (coloured traces corresponding to legend). Distances between distributions are shown as annotations in plots [Kolmogorov–Smirnov (KS) D, P<0.001 for all reported D]. Green/transparent banding on plots indicates boundaries of selected intensity bins. AU, arbitrary units; HA, HAoECs; HP, HPMECs; HU, HUVECs. (B) Percentage of cells per image (microscopy field) pertaining to different intensity bins (nN1) and EC types. (C) Percentage of cells per image (microscopy field) pertaining to different intensity bins (nHES1) and EC types. (D) Intensity of nuclear HES signal (single cells) across phases of cell cycle, cell types (HAoEC, HPMEC, HUVEC) and treatment (CTRL, VEGF). The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Points outside this range are indicated as outlier dots. n of observations for each statistical comparison indicated as annotations under or in individual plots. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 6.

Analysis of NOTCH activation. (A) Scaled density distributions of nN1 and nHES1 for HAoECs, HPMECs and HUVECs (coloured traces corresponding to legend). Distances between distributions are shown as annotations in plots [Kolmogorov–Smirnov (KS) D, P<0.001 for all reported D]. Green/transparent banding on plots indicates boundaries of selected intensity bins. AU, arbitrary units; HA, HAoECs; HP, HPMECs; HU, HUVECs. (B) Percentage of cells per image (microscopy field) pertaining to different intensity bins (nN1) and EC types. (C) Percentage of cells per image (microscopy field) pertaining to different intensity bins (nHES1) and EC types. (D) Intensity of nuclear HES signal (single cells) across phases of cell cycle, cell types (HAoEC, HPMEC, HUVEC) and treatment (CTRL, VEGF). The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Points outside this range are indicated as outlier dots. n of observations for each statistical comparison indicated as annotations under or in individual plots. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 6C shows the analysis for nHES1 where cells were binned according to a log10 scale (bins referred to as low, intermediate, high and very high hereafter). HAoECs and HPMECs had both a lower average percentage of cells in the low intensity category and a greater percentage of cells in the high intensity category (13.7% and 56.9% in HAoECs, and 13.5% and 51% in HPMECs, respectively) in comparison to HUVECs (22.3% and 46.1%). VEGF treatment induced a marked decrease in the average percentage of low intensity cells and a corresponding increase in high intensity cells in HAoECs (6.8% and 61.5%) and had the opposite effect in HUVECs (30.8% and 44.7%).

We then evaluated the effect of differential time in culture by comparing nHES1 distribution in cells cultured for 48 h or 96 h.  Fig. S1A displays density distributions of nHES1 intensity, which demonstrated marked changes in HAoECs and HUVECs and less difference in HPMECs. Fig. S1B shows that the overall effect of longer culture was a marked reduction of cells in the low intensity bin and a corresponding increase in the high bin in HAoECs and HUVECs, whereas HPMECs were mostly unaffected.

Altogether, these results are in line with baseline gene expression levels of HES1 in the three EC types (Fig. S2) and with previous gene expression studies (Chi et al., 2003) where arterial ECs demonstrated higher levels of NOTCH signalling than venous ECs.

We also demonstrated previously unappreciated high levels of NOTCH signalling in HPMECs, which were higher than HUVECs, as also confirmed by gene expression analysis (Fig. S2).

Importantly, our analysis highlights that: (1) NOTCH signalling is not homogeneous in ECs within monolayers, (2) all monolayers examined have similar ranges of nN1 and nHES1 signal intensities, and (3) bulk differences (e.g. qRT-PCR, Fig. S2) result from different proportions of low, intermediate or high signalling cells.

To estimate functional consequences of these observations, we evaluated correlation between downstream NOTCH signalling and the cell cycle. We found that the single-cell levels of nHES1 were increased in cells engaged in cell cycle (∼2-fold increase between G0/G1 and S and a further ∼2-fold increase between S and G2/M) and decreased in dividing cells (∼2-fold among LM and G0/G1) (Fig. 6D). These observations were invariant across cell types, treatment and time in culture. It might be expected that cycling ECs would have low NOTCH signalling as it has been previously shown that NOTCH signalling causes cell cycle arrest and inhibition of proliferation (Fang et al., 2017; Luo et al., 2021). In fact, all cell types in all phases of cell cycle had a significant proportion of cells with low nHES1 (Fig. 6D, shown as outlier dots in boxplots) and mitotic LM cells had an average or low level of nHES1.

We conclude from all data on nHES1 that, as expected and reported before (Fang et al., 2017), confluent ECs have relatively high basal levels of NOTCH signalling. However, we also observed that, within the same population, some cells have low levels of nHES1. To interpret previous reports and our current data, we hypothesize that sustained NOTCH signalling in confluent EC monolayers acts as a molecular break to limit cell cycle progression. Further, we postulate that this exists as an escape mechanism by which cells can have low local levels of nHES1 allowing cell cycle finalisation despite sustained signalling at population level. Many previous studies have highlighted the role of NOTCH in mediating tissue patterning by mechanisms of lateral inhibition and in certain context lateral induction. We therefore set out to estimate the involvement of these processes using our ECPT.

Spatial autocorrelation analysis of NOTCH signalling

To evaluate the role of lateral inhibition and lateral induction mechanisms, we used spatial autocorrelation analysis (Moran's analysis) to assess the distribution of nN1 and nHES1 signals in cell neighbourhoods and across entire monolayers. We first ran the population level analysis (global Moran's, as detailed in Fig. S3 and Materials and Methods) across all images in our dataset using either nN1 or nHES1 signals as inputs and length of junction object as the weighting parameter since estimated strength of neighbours' interaction depends on extent of cell–cell contact.

Fig. 7A shows that, on average, 25.7–40.4% of microscopic fields across cell types and VEGF treatment had a statistically significant positive Moran's index (pGMi). Furthermore, 1.9–4.7% of microscopic fields had a statistically significant negative Moran's Index (nGMi), whereas 57.1–71.7% had randomly distributed cell intensities for nHES1. HUVECs had the highest percentage of microscopic fields with a pGMI.

Fig. 7.

Spatial autocorrelation analysis. (A) Box plots representing percent of fields with either significant (P<0.05) positive or negative global Moran's index (pGMi, nGMi) or random cell distributions (P>0.05 in both pGMI and nGMi) by cell type and treatment. (B) Representative maps of cells with significant pMGi and relative Local Moran's analysis. Green coloured maps refer to intensity of either nuclear NOTCH1 (log10) or HES1 (log2) signal (white=low intensity, dark green high intensity). Red coloured maps refer to P-value of local Moran's index LMi (white P<0.05, dark red P=1). Outlines in intensity maps indicated neighbourhood of selected cells with significant LMi associated with either positive (black) or negative (red) LMi. (C) Boxplot (by EC type) representing percentages of dividing cells per images grouped in bins with low <2 or higher (2–10) numbers of cells with negative and significant LMi per field. The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each statistical comparison indicated as annotations under or in individual plots. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Fig. 7.

Spatial autocorrelation analysis. (A) Box plots representing percent of fields with either significant (P<0.05) positive or negative global Moran's index (pGMi, nGMi) or random cell distributions (P>0.05 in both pGMI and nGMi) by cell type and treatment. (B) Representative maps of cells with significant pMGi and relative Local Moran's analysis. Green coloured maps refer to intensity of either nuclear NOTCH1 (log10) or HES1 (log2) signal (white=low intensity, dark green high intensity). Red coloured maps refer to P-value of local Moran's index LMi (white P<0.05, dark red P=1). Outlines in intensity maps indicated neighbourhood of selected cells with significant LMi associated with either positive (black) or negative (red) LMi. (C) Boxplot (by EC type) representing percentages of dividing cells per images grouped in bins with low <2 or higher (2–10) numbers of cells with negative and significant LMi per field. The box represents the 25–75th percentiles, and the median (line) and mean (circle) is indicated. The whiskers extend to the maximum or minimum values within a 1.5× interquartile range (IQR) distance from the box hinge. Outliers not represented. n of observations for each statistical comparison indicated as annotations under or in individual plots. *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001 (one-way ANOVA followed by Tukey's post-hoc HSD test).

Thus, the first key finding upon global Moran's analysis is that most microscopic fields analysed (irrespective of cell type and treatment) had randomly distributed cell intensities, highlighting the presence of intermediate cell states at the local population level. However, a sizeable fraction of the populations analysed shown evidence of clustered or sparse cell distributions, suggesting that lateral inhibition and lateral induction are effectively active in the system. The average positive and negative global Moran's indexes for all images analysed was ∼0.4 and ∼−0.3. These are intermediate values (extremes are 1 and −1 as illustrated in Fig. S3), which suggest that clustering or sparse distributions were detected in cell neighbourhoods smaller than the microscopic field analysed. Overall, global Moran's analysis confirms that the observed bulk population distributions of nHES1 (Fig. 6A) are reflective of local arrangement after considering spatial relationships, and suggests that clustering or sparseness can emerge in the wider population. To confirm and extend these observations, we performed local Moran's analysis, which considers the immediate cell neighbourhood.

Fig. 7B illustrates one example of such analysis for both NOTCH (upper panels) and HES1 (lower panels). The intensity maps associated with nN1 and nHES1 highlight local heterogeneity in signalling and the P-value maps indicate which cells in the map have statistically significant local Moran's index (LMi, either positive or negative).

Comparisons of intensity and P-value maps highlight clusters of cells with either similar or dissimilar signal values, which were predicted by global Moran's analysis and confirms that these clusters are smaller than the microscopic field (5–20 cells). 84% of all images had at least one cell with a statistically significant negative LMI (nHES1) and 50% of all images had at least one cell with a statistically significant positive LMI (nHES1). The frequency of nLMi were 1.6±1.2, 1.4±1.3 and 2.8±1.8 nLMI cells/per field (mean±s.d.) in HAoECs and HPMECs and HUVECs. Average frequencies of pLMi were 0.6±0.8, 0.6±1 and 1±1.4 cells/field for HAoECs, HPMECs and HUVECs, and VEGF treatment did not induce substantial changes in nLMi or pLMi frequencies.

These observations confirm those of the global analysis and show that the distribution of nHES1 within EC monolayers is largely random and independent of cell type or exposure to VEGF treatment. However, monolayers are interspersed with neighbourhoods of cells that are engaged in lateral induction or inhibition processes as exemplified in Fig. 7B.

If cell proliferation is correlated with these phenomena as hypothesized above, we expect that cells able to progress through the cell cycle would be laterally inhibited to express lower levels of HES1 and escape cell cycle arrest. If that were the case, we would also expect to find more cells progressing through the cell cycle in areas where lateral inhibition is prevalent. To test this hypothesis within the limits of our experimental setup, we compared microscopy fields with below – or above – average number of nLMi cells. Fig. 7C shows the results of this analysis, confirming that microscopy fields with above average numbers of nLMI cells also had higher numbers of cells in cell cycle (Fig. 7C). As expected, no difference was found when considering abundance of pLMi cells.

Overall, we conclude that local lateral inhibition mechanisms can inhibit HES1 expression in a few cells within an EC population and thus increase the probability of individual cells progressing through the cell cycle.

Finally, as we had observed that longer times in culture quantitatively affected the distribution of nHES1 intensities (Fig. S1A,B), we sought to evaluate qualitative effects in Lmi. Fig. S1C shows comparisons between counts of images containing 0 to 8 cells with negative Lmi (nLMi) and demonstrates that longer time in culture reduced the proportion of images with higher number of nLMi cells (without change in pLMi cell counts). Overall, longer time in culture reduced the number of cells engaged in lateral inhibition without subverting the overall distributions of cells with scattered presence of either small homogeneous cell clusters or laterally inhibited cells.

ECs exert an outstanding variety of specialised functions (Augustin and Koh, 2017; Rafii et al., 2016), which are reflected in EC phenotypical heterogeneity. EC diversity arises during development and is preserved through homeostasis. The activation of differential genetic programmes interplays with microenvironmental factors (Adams and Alitalo, 2007; Chi et al., 2003). Interestingly heterogeneity within the same vascular bed, where microenvironmental factors are likely only undergoing minor fluctuations, has been observed in vivo (McCarron et al., 2019) and deemed crucial to understand variable functions within the same EC type.

These mechanisms require similar cells to preserve discrete, diverse and adaptable phenotypes, and the time scale at which these changes happen can be smaller than that of gene expression (Chapman et al., 2016; Stepanova et al., 2021). Thus, to justify rapid phenotype switches it is necessary to involve inherently dynamic processes, such as junctional plasticity. Direct evaluation of dynamic changes in cell connectivity still poses considerable technical hurdles, and despite great advances in live-microscopy, comprehensive cell phenotyping of single cells at the population level still relies heavily on end-point experiments. Although observation of fixed cells cannot provide information on the exact dynamics of the system, extensive profiling of single cells can inform about cell dynamic behaviours and help to develop hypotheses, which can then be evaluated experimentally or with computational methods (Altschuler and Wu, 2010).

To gain a better understanding of EC heterogeneity, connectivity and emerging dynamic behaviours (McCarron et al., 2019, 2017), we developed ECPT and characterised three EC lines under standard conditions or supra-physiological levels of VEGF (like those observed in cancer). In line with previous studies (Chi et al., 2003), our data demonstrate that human aortic, umbilical vein and pulmonary microvascular ECs present distinct phenotypes (Figs 17), which are maintained at steady state and are independent of microenvironmental stimuli. HAoECs, HPMECs and HUVECs displayed different intrinsic proliferative potential and demonstrated a diffuse ability to proliferate, in contrast to the idea that ECs undergo terminal differentiation and lose proliferative capacity (Yoder, 2018). Our image-based analysis of cell cycle based on DNA intensity highlighted interesting differences between the different EC lines in subconfluent conditions, with HUVECs and HAoECs proliferating faster than HPMECs, and an overall small increase in the dividing cell number after VEGF treatment. In confluent conditions, all cell types settled to a lower proliferation rate, which was similar for all EC types.

It is well established that cytoskeleton arrangement correlates with differential EC phenotypes in different contexts. ECs exposed to flow display stress fibres aligned with direction of flow and migrating ECs display similar stress fibres along the direction of migration. Analysis of the EC cytoskeleton (Fig. 4) demonstrated distinctive profiles in line with proliferation statuses in HAoECs, HPMECs and HUVECs, with microvascular HPMECs having the lowest numbers of stress fibres. At the same time, all cells demonstrated intra-population heterogeneity with intermixing of cells with or without stress fibres in individual monolayers (Fig. 4). This suggests a high degree of cellular rearrangement as indicated by time-lapse microscopy experiments showing that live ECs in monolayers constantly rearrange their shape and connections (Kim and Cooper, 2018).

Differential EC functions lead to variance in IEJs, which in turn regulate traffic of molecules and solutes in and out of capillaries and inhibit coagulation by preventing exposure of the underlying subintimal layer in arteries and veins. Analysis of IEJ in three different EC populations highlighted a high degree of junctional heterogeneity. In line with the data discussed above and existing literature, we show that HAoECs (arterial) had more linear junctions than HUVECs (venous), which instead tended to form highly dynamic jagged junctions. Importantly, HPMECs had the most linear/stabilised junctions, whereas HAoECs had more heterogeneous jagged/linear junctions and stress fibres (Fig. 4). This is reminiscent of the unique role the pulmonary and arterial vessels play in regulating solute exchange or adapting to high flow, respectively. This is in line with recent sc-RNA sequencing data demonstrating that >40% of murine HAoEC express mesenchymal genes (Lukowski et al., 2019) and are therefore expected to display morphological and junctional heterogeneity accompanied by stress fibres.

In agreement with results reported before (Bazzoni and Dejana, 2004), longer culture times and VEGF treatment in our experiments induced characteristic changes in cell morphology and promoted more continuous junction types across all EC types.

Distinctive features among different EC populations seem to be hard coded within cell gene expression programme, as cells from individual donors had vastly overlapping phenotypes and the features were maintained in all lines upon several passages in culture.

If the decision to enter the cell cycle or to maintain stable junctions were dependent on clear-cut differential gene expression, we would expect to find a limited collection of homogeneous phenotypes. However, we find a continuous range of EC phenotypes within the same EC monolayer. To interpret EC intra-population heterogeneity, we focussed on NOTCH signalling for its well-established role as a coordinator of the opposing functions of EC proliferation (Fang et al., 2017; Luo et al., 2021) and junctional complex stabilisation (Bentley et al., 2014).

Furthermore, the NOTCH pathway is one of the main drivers of endothelial cell heterogeneity and is linked to vascular maturation, arteriovenous specification and angiogenic sprouting (Fish and Wythe, 2015; Hellström et al., 2007; Potente et al., 2011; Torres-Vázquez et al., 2003). Although NOTCH signalling dynamics have been studied in HUVECs and retina, it is still unclear how it affects EC heterogeneity in different organs. Nuclear NOTCH1 (N1-ICD; where ICD is intracellular domain, generically NICD) acts as a transcription co-factor to induce the expression of target genes, such as the HES/HEY family of transcription factors, and is therefore active in the nucleus. In line with previous results on HEY1 and HEY2 expression in HUVECs (Aranguren et al., 2013), we hypothesised that ECs within monolayers would have a homogeneous level of NOTCH signalling that was either low or absent in venous ECs and high in arterial ECs. One study has reported appreciable baseline NICD levels in HPMECs by western blotting (Zong et al., 2018) but not in comparison with other EC lines. Recent single-cell RNA sequencing data of human lungs has highlighted heterogeneity among arterial, venous, and microvascular ECs, which displayed intermediate phenotypes (Kalucka et al., 2020). Interestingly, activation of NOTCH downstream target genes also seems heterogeneous across and within EC populations (Travaglini et al., 2020). Fluctuations in NOTCH signalling have been hypothesised in EC monolayers and demonstrated in EC during sprouting angiogenesis where NOTCH acts as a bistable switch (Ubezio et al., 2016). It has been demonstrated that leading tip cells induce NOTCH signalling in trailing stalk cells via Dll4–NOTCH1 leading to lateral inhibition of the tip cell phenotype (Ubezio et al., 2016), whereas ECs within a stabilised monolayer cannot acquire a tip cell phenotype and all ECs receive, in principle, similar NOTCH stimulation from neighbours.

By measuring NOTCH1 and HES1 levels in individual ECs, we found that different ECs within their monolayer had different levels of NOTCH activation, which was also reflected in different bulk gene expression measures of NOTCH target genes (HES1). However, the clearest differences across cell types and treatments were found in relative proportions of cells with differential signal within the same EC population, suggesting that individual cells can acquire differential phenotypes, which can be modulated by VEGF treatment (Fig. 6). This degree of heterogeneity in NOTCH signalling in neighbouring cells also strongly suggests that the NOTCH phenotype is regulated dynamically in EC monolayers, with lateral inhibition and lateral induction being two potential candidate mechanisms. To investigate this hypothesis, we performed spatial autocorrelation analysis (Fig. 7), which demonstrated a high degree of heterogeneity of NOTCH and HES1 in ECs within the same monolayer. If either lateral inhibition or lateral induction were the prevailing mechanism in regulating the NOTCH phenotype in ECs, we would expect sparse or uniform spatial distribution of NOTCH and HES1 activation. However, we found a high degree of randomness in the spatial distribution of these signals. Alternatively, the two mechanisms might act in concert to produce qualitatively different cell distributions, which seems plausible as all EC types analysed expressed both Dll4 and Jagged1 ligands (Fig. S2), which have been shown to be involved in lateral inhibition and lateral induction, respectively (Pedrosa et al., 2015). When we performed local spatial autocorrelation analysis, we confirmed that both mechanisms seem to occur together and our results are in line with previous mathematical models (Boareto et al., 2016) showing that concurrent lateral inhibition and induction can generate intermediate cell states. Our current data demonstrate that intermediate phenotypes in both nN1 and nHES1 are common in all ECs analysed, irrespective of treatment, and that lateral induction or lateral inhibition patterns emerge locally. This complex spatial distribution does not seem compatible with stabilised phenotypes and rather suggests a scenario where cell phenotype is dynamically regulated and can change over time to exert differential functions while maintaining a stable balance across the wider population.

Novel spatialised mathematical models of NOTCH signalling could be used in future work to assess whether NOTCH is sufficient to generate the heterogeneity we found in our experiments or whether further layers of regulation need to be accounted for.

To evaluate the functional consequences of our findings, we asked whether spatial patterning accompanied by differential signalling could affect other parameters in our dataset. We did not find any correlation between the extent or spatial distribution of nN1 and nHES and junction status, possibly because junctions are remodelled at a very fast pace in cultured ECs (Kim and Cooper, 2018) and thus our experimental context using fixed cells might fail to resolve processes at this timescale. When we considered cell proliferation, we found that HES1 was, on average, higher in cells progressing through the cell cycle but several cells had low levels irrespective of their cell cycle status. When we compared the abundance of laterally inhibited cells with dividing cells at the population level, we consistently found that populations containing more of these cells were also proliferating at a slightly higher pace. Overall, our data suggest that the decision to initiate and progress through cell cycle in continuous monolayers with high basal levels of NOTCH signalling is regulated by the extent (and possibly duration) of lateral inhibition in the local cell neighbourhood. Moreover, we show that the formation of patches of cells with similar NOTCH signalling is a common finding in EC monolayers in vitro. It will be interesting to investigate whether this is reflected in vivo and the functional consequences of this phenomenon considering previous work has identified patches of cells with differential Ca2+ signalling and density of endothelial M3 muscarinic acetylcholine receptors (AchRM3s) in vivo (McCarron et al., 2019, 2017). However, in vivo investigations will still pose significant technical challenges to high throughput analyses such as those presented here.

It will be also interesting in future work to evaluate whether the patterns observed in our experiments are stable or remodelled over time, and the timescales of these changes. The ECPT presented in this work is a step forward in comparison to available platforms (see Table 1) and represents a bridge between experimental and computational setups, where experimental ECPT data can be used to guide development of spatialised mathematical models, which in turn can be used to guide further experimentations. The experimental setup introduced in the present work does not extend to resolving highly dynamic processes; however, ECPT paves the way towards quantifying more advanced experiments, such as measurements of gene transcription by fluorescence in situ hybridization (FISH) or live imaging using fluorescent reporters.

Overall, the ECPT allows for the evaluation of cell-intrinsic mechanisms of monolayer maintenance and plasticity, excluding variability caused by microenvironmental factors, in line with previous observations on blood vessel heterogeneity in vivo originating from EC rather than perivascular cells (Chavkin and Hirschi, 2020). It will be important to elucidate in future studies how the crosstalk with perivascular cells and microenvironment can affect emerging endothelial behaviours.

We envisage that coupling our ECPT workflow with live imaging setups, computational modelling and single cells transcriptomic analysis will open the way to a much deeper understanding of emerging dynamic endothelial behaviours and thus help to develop novel more effective therapies for regenerative medicine, prevention of cardiovascular diseases and the treatment of cancer.

Cell culture

HAoECs, HPMECs and HUVECs (all PromoCell) were plated on 10 μg/ml fibronectin (from human plasma, Promocell)-coated flasks, grown in EGMV2 medium (Promocell) in the absence of antibiotics, detached with accutase (Thermo Fisher Scientific, Waltham, MA), and used by passage 5. We analysed two distinct donors for each cell type, which were chosen excluding diseases affecting the vasculature (e.g. diabetes). Donor's age (HAoECs, HPMECs) was between 50 and 63 years. For experiments, 4×104 ECs per well were seeded in fibronectin-coated 96-well plates (µClear, Greiner) and cultured for 48 or 96 h under basal (EGMV2, Promocell) or activated (EGMV2+50 ng/ml VEGFA, Peprotech, London, UK) conditions in triplicate paralleling conditions described previously (Andriopoulou et al., 1999). The EC formed confluent monolayers at microscopic inspection [phase contrast, 10×–20× original magnification (OM; as indicated on objective)] at the time of immunostaining and image acquisition.

Immunostaining

Cells were fixed with 2% paraformaldehyde in phosphate-buffered saline (PBS) for 10 min at room temperature. Cells were blocked 1 h with PBS supplemented with 1% fetal bovine serum (FBS) and permeabilised with 0.1% Triton X-100. Cells were then incubated for 1 h at room temperature with primary antibodies against CDH5 (VE-cadherin; Novusbio NB600-1409, 1 μg/ml final), NOTCH1 (Abcam, ab194122, Alexa Fluor 647-conjugated, 1 μg/ml final) and Hes1 (Abcam, ab119776, 1 μg/ml final). Plates were washed and incubated 1 h with 1 μg/ml secondary Alexa Fluor 488-conjugated and Alexa Fluor 555-conjugated antibody (Thermo Fisher Scientific), Hoechst 33342 (1 μg/ml, Sigma) and Phalloidin-Atto 647N (Sigma). All antibodies used in immunostaining were previously validated for specificity by western blot analysis (Fig. S4).

Image acquisition

We obtained images from slides with an Operetta CLS system (PerkinElmer, Waltham, MA) equipped with a 40× water immersion lens (NA 1.1). In each well, three areas were acquired. Each area is composed of nine microscopic fields at 40× OM (Fig. S1). We standardised acquisition parameters (light-emitting diode power, exposure time) throughout different experiments and used HUVECs as a standard for calibration in all experiments. We analysed an image database containing 28,000 images (7000 fields in four fluorescence channels) extracted from nine independent experiments conducted on EC lines from two different donors each. Two intraexperiment replicates were conducted for each experiment.

Endothelial cell profiling tool

We used a combination of machine learning-aided image segmentation (ImageJ) (Schindelin et al., 2012) and an image-based cell profiling tool (CellProfiler) (Carpenter et al., 2006) to extract the phenotype of single EC in monolayers. Our workflow enables the measurement of EC morphology (area, perimeter, shape descriptors and cell neighbours), NOTCH1, HES1, CDH5 and DNA intensities, and the characterisation of inter-endothelial junctions (IEJs). In the present study, we chose to analyse only selected features with recognised functions in EC biology. Image texture features were measured and only used during the training of ML algorithms (Caicedo et al., 2017; Jones et al., 2008), which in turn were used to classify junction morphology and LM cells. ECPT scripts and methods including FIJI/ImageJ macros for image importing and pre-processing are detailed in Appendix 1 at https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers.

We imported and cleaned the results into R studio excluding artefacts and mis-segmented cell objects (extreme values in cell area or signal intensity, NAs in measurements). We then calculated continuous and categorical (cell cycle) parameters. Following guidelines suggested by Caicedo et al. (2017), we pre-processed the database to exclude mis-segmented cells and to normalise the measurements prior to dimensionality reduction or single factor analysis. We reformatted and tidied the database and calculated summary statistics using packages from the Tidyverse packages collection (https://www.tidyverse.org; Dplyr, Tidyr, Tibble, Forcats, Purr and Ggplot2).

We conducted PCA analysis on ECPT parameters (Cell morphology, intensity, neighbourhood and junctional status,  Table S1) using the prcomp function (R Stats package) generating 12 PCA components. All plots in figures are generated using the Ggplot2 R package.

We created a Shiny application (https://CRAN.R-project.org/package=shiny) for interactive selection and visualisation of data of interest. All files containing the code required to reproduce all the plots in the paper and to run the Shiny application using our dataset are available on https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers.

Western blotting

Cells were scraped in the presence of RIPA buffer (Sigma-Aldrich) containing protease (Millipore, UK) and phosphatase inhibitors (Sigma-Aldrich, UK), left on ice for 15 min, and centrifuged at 13,000 g for 5 min in a refrigerated microfuge. Supernatants were assessed for total protein using the BCA protein quantitation kit (Thermo Fisher Scientific). 15 µg of protein were separated on NuPAGE 4–12% Bis-Tris gels (Invitrogen) before being transferred to nitrocellulose membranes (GE Healthcare, Amersham, UK). After probing with primary and secondary antibodies, membranes were developed using Clarity™ Western ECL Substrate (Bio-Rad, UK) and read using a ChemiDoc system (Bio-Rad). Antibodies for NOTCH1 (1:500; Abcam, ab194122), HES1 (1:1000; Abcam, ab119776), VE-cadherin (1:500; Novus bio NB600-1409) and β-tubulin (1:1000; Cell Signaling Technology) were used. Goat anti-mouse-IgG and anti-rabbit-IgG horseradish peroxidase (HRP)-conjugated antibodies were from Dako (Agilent).

RNA extraction and qRT-PCR

Total RNA was extracted and purified using the Monarch total RNA miniprep kit according to the manufacturer's instructions. The resulting RNA was quantified using Nanodrop (ISOGEN Life Science). For quantitative real-time PCR (qRT-PCR), 1 μg of RNA was used for reverse transcription using the iScript cDNA synthesis kit (Bio-Rad). The gene expression analysis was carried out using the SsoAdvancedTM Universal SYBR Green Supermix (Bio-Rad) and analysed by means of a Stratagene Mx3000P (Agilent Technologies) in real time, primers used are listed in Table S2.

Statistical analysis

To compare multiple groups, we used one-way ANOVA followed by Tukey's HSD post-hoc test. We considered P<0.05 (*) statistically significant and P<0.01 (**), P<0.001 (***) and P<0.0001(****) highly significant.

To evaluate linear correlation between continuous variables we calculated Pearson's correlation coefficient I and considered P<0.001 as highly statistically significant.

Comparisons between intensity distributions were performed by two-sided Kolmogorov–Smirnov (KS) test as implemented in the ‘stats’ R package. We considered comparisons with P<0.01 to be statistically significant and reported corresponding distance index (D).

Statistical significance in global and local Moran's analyses is computed by random permutations as implemented in the adespatial R package we used 999 permutations in each test. We considered P<0.05 in Gmi or Lmi to be statistically significant (Dray, 2011).

We wish to thank Dr V. La Ferla for creating the FIJI importer macro including GUI. A substantial proportion of these methods have built from previous work funded by the Wellcome Trust and the UK Medical Research Council (MRC) through the Human Induced Pluripotent Stem Cell Initiative (WT098503). D.D. also gratefully acknowledges funding from the Department of Health via the National Institute for Health Research comprehensive Biomedical Research Centre award to Guy's & St. Thomas' National Health Service Foundation Trust in partnership with King's College London and King's College Hospital NHS Foundation Trust.

Author contributions

Conceptualization: F.C., L.V.; Methodology: F.C., D.D., L.V.; Software: J.H., E.R., M. Branco, R.S., L.V.; Validation: F.C., J.H., J.L.C., L.V.; Formal analysis: L.V.; Investigation: F.C., A.P., J.L.C., M. Battilocchi, E.E.; Resources: D.D.; Data curation: E.R., L.V.; Writing - original draft: F.C., E.R., J.L.C., D.D., L.V.; Writing - review & editing: F.C., J.H., E.R., J.L.C., D.D., L.V.; Visualization: F.C., M. Branco, R.S., E.E., L.V.; Supervision: D.D., L.V.; Project administration: L.V.

Funding

This work is supported by an internal King's College London Dental Institute seed fund awarded to L.V. with D.D. as collaborator. F.C. is supported by a studentship from King's College London, Faculty of Dentistry Oral and Craniofacial Sciences (FoDOCS).

Data availability

All raw and elaborated data are available at https://github.com/exr98/HCA-uncovers-metastable-phenotypes-in-hEC-monolayers. To facilitate data exploration and re-use we have developed a dedicated data browser using the Shiny app environment. Original images are available upon request to the authors.

Abu Taha
,
A.
,
Taha
,
M.
,
Seebach
,
J.
and
Schnittler
,
H.-J.
(
2014
).
ARP2/3-mediated junction-associated lamellipodia control VE-cadherin-based cell junction dynamics and maintain monolayer integrity
.
Mol. Biol. Cell
25
,
245
-
256
.
Acar
,
M.
,
Mettetal
,
J. T.
and
van Oudenaarden
,
A.
(
2008
).
Stochastic switching as a survival strategy in fluctuating environments
.
Nat. Genet
40
,
471
-
475
.
Adams
,
R. H.
and
Alitalo
,
K.
(
2007
).
Molecular regulation of angiogenesis and lymphangiogenesis
.
Nat. Rev. Mol. Cell Biol.
8
,
464
-
478
.
Altschuler
,
S. J.
and
Wu
,
L. F.
(
2010
).
Cellular heterogeneity: do differences make a difference?
Cell
141
,
559
-
563
.
Andriopoulou
,
P.
,
Navarro
,
P.
,
Zanetti
,
A.
,
Lampugnani
,
M. G.
,
Dejana
,
E.
(
1999
).
Histamine induces tyrosine phosphorylation of endothelial cell-to-cell adherens junctions
.
Arterioscler. Thromb. Vasc. Biol
19
,
2286
-
2297
.
Aranguren
,
X. L.
,
Beerens
,
M.
,
Coppiello
,
G.
,
Wiese
,
C.
,
Vandersmissen
,
I.
,
Lo Nigro
,
A.
,
Verfaillie
,
C. M.
,
Gessler
,
M.
and
Luttun
,
A.
(
2013
).
COUP-TFII orchestrates venous and lymphatic endothelial identity by homo- or hetero-dimerisation with PROX1
.
J. Cell Sci.
126
,
1164
-
1175
.
Augustin
,
H. G.
and
Koh
,
G. Y.
(
2017
).
Organotypic vasculature: From descriptive heterogeneity to functional pathophysiology
.
Science
357
,
eaal2379
.
Bazzoni
,
G.
and
Dejana
,
E.
(
2004
).
Endothelial cell-to-cell junctions: molecular organization and role in vascular homeostasis
.
Physiol. Rev.
84
,
869
-
901
.
Bentley
,
K.
,
Franco
,
C. A.
,
Philippides
,
A.
,
Blanco
,
R.
,
Dierkes
,
M.
,
Gebala
,
V.
,
Stanchi
,
F.
,
Jones
,
M.
,
Aspalter
,
I. M.
,
Cagna
,
G.
et al. 
(
2014
).
The role of differential VE-cadherin dynamics in cell rearrangement during angiogenesis
.
Nat. Cell Biol
16
,
309
-
321
.
Boareto
,
M.
,
Jolly
,
M. K.
,
Goldman
,
A.
,
Pietilä
,
M.
,
Mani
,
S. A.
,
Sengupta
,
S.
,
Ben-Jacob
,
E.
,
Levine
,
H.
and
Onuchic
,
J. N.
(
2016
).
Notch-Jagged signalling can give rise to clusters of cells exhibiting a hybrid epithelial/mesenchymal phenotype
.
J. R. Soc. Interface
13
,
20151106
.
Brezovjakova
,
H.
,
Tomlinson
,
C.
,
Mohd Naim
,
N.
,
Swiatlowska
,
P.
,
Erasmus
,
J. C.
,
Huveneers
,
S.
,
Gorelik
,
J.
,
Bruche
,
S.
and
Braga
,
V. M.
(
2019
).
Junction Mapper is a novel computer vision tool to decipher cell–cell contact phenotypes
.
eLife
8
,
e45413
.
Caicedo
,
J. C.
,
Cooper
,
S.
,
Heigwer
,
F.
,
Warchal
,
S.
,
Qiu
,
P.
,
Molnar
,
C.
,
Vasilevich
,
A. S.
,
Barry
,
J. D.
,
Bansal
,
H. S.
,
Kraus
,
O.
et al. 
(
2017
).
Data-analysis strategies for image-based cell profiling
.
Nat. Methods
14
,
849
-
863
.
Carpenter
,
A. E.
,
Jones
,
T. R.
,
Lamprecht
,
M. R.
,
Clarke
,
C.
,
Kang
,
I. H.
,
Friman
,
O.
,
Guertin
,
D. A.
,
Chang
,
J. H.
,
Lindquist
,
R. A.
,
Moffat
,
J.
et al. 
(
2006
).
CellProfiler: image analysis software for identifying and quantifying cell phenotypes
.
Genome Biol.
7
,
R100
.
Chapman
,
G.
,
Major
,
J. A.
,
Iyer
,
K.
,
James
,
A. C.
,
Pursglove
,
S. E.
,
Moreau
,
J. L. M.
and
Dunwoodie
,
S. L.
(
2016
).
Notch1 endocytosis is induced by ligand and is required for signal transduction
.
Biochim. Biophys. Acta
1863
,
166
-
177
.
Chavkin
,
N. W.
and
Hirschi
,
K. K.
(
2020
).
Single cell analysis in vascular biology
.
Front. Cardiovasc. Med.
7
,
42
.
Chi
,
J.-T.
,
Chang
,
H. Y.
,
Haraldsen
,
G.
,
Jahnsen
,
F. L.
,
Troyanskaya
,
O. G.
,
Chang
,
D. S.
,
Wang
,
Z.
,
Rockson
,
S. G.
,
van de Rijn
,
M.
,
Botstein
,
D.
et al. 
(
2003
).
Endothelial cell diversity revealed by global expression profiling
.
Proc. Natl. Acad. Sci. U. S. A
100
,
10623
-
10628
.
Dray
,
S.
(
2011
).
A New Perspective about Moran's Coefficient: Spatial Autocorrelation as a Linear Regression Problem. Moran系数的新视角:空间自相关视为线性回归问题: Spatial Autocorrelation as a Linear Regression Problem
.
Geogr. Anal
43
,
127
-
141
.
Dueck
,
H.
,
Eberwine
,
J.
and
Kim
,
J.
(
2016
).
Variation is function: Are single cell differences functionally important?: Testing the hypothesis that single cell variation is required for aggregate function
.
BioEssays News Rev. Mol. Cell. Dev. Biol.
38
,
172
-
180
.
Eling
,
N.
,
Morgan
,
M. D.
and
Marioni
,
J. C.
(
2019
).
Challenges in measuring and understanding biological noise
.
Nat. Rev. Genet
20
,
536
-
548
.
Fang
,
J. S.
,
Coon
,
B. G.
,
Gillis
,
N.
,
Chen
,
Z.
,
Qiu
,
J.
,
Chittenden
,
T. W.
,
Burt
,
J. M.
,
Schwartz
,
M. A.
and
Hirschi
,
K. K.
(
2017
).
Shear-induced Notch-Cx37-p27 axis arrests endothelial cell cycle to enable arterial specification
.
Nat. Commun.
8
,
2149
.
Fernández-Martín
,
L.
,
Marcos-Ramiro
,
B.
,
Bigarella
,
C. L.
,
Graupera
,
M.
,
Cain
,
R. J.
,
Reglero-Real
,
N.
,
Jiménez
,
A.
,
Cernuda-Morollón
,
E.
,
Correas
,
I.
,
Cox
,
S.
et al. 
(
2012
).
Crosstalk between reticular adherens junctions and platelet endothelial cell adhesion molecule-1 regulates endothelial barrier function
.
Arterioscler. Thromb. Vasc. Biol
32
,
e90
-
102
.
Ferrero
,
E.
,
Scabini
,
S.
,
Magni
,
E.
,
Foglieni
,
C.
,
Belloni
,
D.
,
Colombo
,
B.
,
Curnis
,
F.
,
Villa
,
A.
,
Ferrero
,
M. E.
and
Corti
,
A.
(
2004
).
Chromogranin A protects vessels against tumor necrosis factor α–induced vascular leakage
.
FASEB J.
18
,
554
-
556
.
Fish
,
J. E.
and
Wythe
,
J. D.
(
2015
).
The molecular regulation of arteriovenous specification and maintenance
.
Dev. Dyn. Off. Publ. Am. Assoc. Anat
244
,
391
-
409
.
Ginzberg
,
M. B.
,
Chang
,
N.
,
D'Souza
,
H.
,
Patel
,
N.
,
Kafri
,
R.
and
Kirschner
,
M. W.
(
2018
).
Cell size sensing in animal cells coordinates anabolic growth rates and cell cycle progression to maintain cell size uniformity
.
eLife
7
,
e26957
.
Hellström
,
M.
,
Phng
,
L.-K.
,
Hofmann
,
J. J.
,
Wallgard
,
E.
,
Coultas
,
L.
,
Lindblom
,
P.
,
Alva
,
J.
,
Nilsson
,
A.-K.
,
Karlsson
,
L.
,
Gaiano
,
N.
et al. 
(
2007
).
Dll4 signalling through Notch1 regulates formation of tip cells during angiogenesis
.
Nature
445
,
776
-
780
.
Itkin
,
T.
,
Gur-Cohen
,
S.
,
Spencer
,
J. A.
,
Schajnovitz
,
A.
,
Ramasamy
,
S. K.
,
Kusumbe
,
A. P.
,
Ledergor
,
G.
,
Jung
,
Y.
,
Milo
,
I.
,
Poulos
,
M. G.
et al. 
(
2016
).
Distinct bone marrow blood vessels differentially regulate haematopoiesis
.
Nature
532
,
323
-
328
.
Jones
,
T. R.
,
Carpenter
,
A. E.
,
Lamprecht
,
M. R.
,
Moffat
,
J.
,
Silver
,
S. J.
,
Grenier
,
J. K.
,
Castoreno
,
A. B.
,
Eggert
,
U. S.
,
Root
,
D. E.
,
Golland
,
P.
et al. 
(
2009
).
Scoring diverse cellular morphologies in image-based screens with iterative feedback and machine learning
.
Proc. Natl. Acad. Sci. USA
106
,
1826
-
1831
.
Jones
,
T. R.
,
Kang
,
I. H.
,
Wheeler
,
D. B.
,
Lindquist
,
R. A.
,
Papallo
,
A.
,
Sabatini
,
D. M.
,
Golland
,
P.
and
Carpenter
,
A. E.
(
2008
).
CellProfiler Analyst: data exploration and analysis software for complex image-based screens
.
BMC Bioinformatics
9
,
482
.
Kalucka
,
J.
,
de Rooij
,
L. P. M. H.
,
Goveia
,
J.
,
Rohlenova
,
K.
,
Dumas
,
S. J.
,
Meta
,
E.
,
Conchinha
,
N. V.
,
Taverna
,
F.
,
Teuwen
,
L.-A.
,
Veys
,
K.
et al. 
(
2020
).
Single-cell transcriptome atlas of murine endothelial cells
.
Cell
180
,
764
-
779.e20
.
Kim
,
J.
and
Cooper
,
J. A.
(
2018
).
Septins regulate junctional integrity of endothelial monolayers
.
Mol. Biol. Cell
29
,
1693
-
1703
.
Kusumbe
,
A. P.
,
Ramasamy
,
S. K.
and
Adams
,
R. H.
(
2014
).
Coupling of angiogenesis and osteogenesis by a specific vessel subtype in bone
.
Nature
507
,
323
-
328
.
Lampugnani
,
M. G.
,
Resnati
,
M.
,
Raiteri
,
M.
,
Pigott
,
R.
,
Pisacane
,
A.
,
Houen
,
G.
,
Ruco
,
L. P.
and
Dejana
,
E.
(
1992
).
A novel endothelial-specific membrane protein is a marker of cell-cell contacts
.
J. Cell Biol.
118
,
1511
-
1522
.
Lampugnani
,
M. G.
,
Corada
,
M.
,
Caveda
,
L.
,
Breviario
,
F.
,
Ayalon
,
O.
,
Geiger
,
B.
and
Dejana
,
E.
(
1995
).
The molecular organization of endothelial cell to cell junctions: differential association of plakoglobin, beta-catenin, and alpha-catenin with vascular endothelial cadherin (VE-cadherin)
.
J. Cell Biol.
129
,
203
-
217
.
Lord
,
S. J.
,
Velle
,
K. B.
,
Mullins
,
R. D.
and
Fritz-Laylin
,
L. K.
(
2020
).
SuperPlots: Communicating reproducibility and variability in cell biology
.
J. Cell Biol.
219
,
e202001064
.
Lukowski
,
S. W.
,
Patel
,
J.
,
Andersen
,
S. B.
,
Sim
,
S.-L.
,
Wong
,
H. Y.
,
Tay
,
J.
,
Winkler
,
I.
,
Powell
,
J. E.
and
Khosrotehrani
,
K.
(
2019
).
Single-cell transcriptional profiling of aortic endothelium identifies a hierarchy from endovascular progenitors to differentiated cells
.
Cell Rep
27
,
2748
-
2758
.
e3
.
Luo
,
W.
,
Garcia-Gonzalez
,
I.
,
Fernández-Chacón
,
M.
,
Casquero-Garcia
,
V.
,
Sanchez-Muñoz
,
M. S.
,
Mühleder
,
S.
,
Garcia-Ortega
,
L.
,
Andrade
,
J.
,
Potente
,
M.
and
Benedito
,
R.
(
2021
).
Arterialization requires the timely suppression of cell growth
.
Nature
589
,
437
-
441
.
MacLean
,
A. L.
,
Hong
,
T.
and
Nie
,
Q.
(
2018
).
Exploring intermediate cell states through the lens of single cells
.
Curr. Opin. Syst. Biol.
9
,
32
-
41
.
McCarron
,
J. G.
,
Lee
,
M. D.
and
Wilson
,
C.
(
2017
).
The endothelium solves problems that endothelial cells do not know exist
.
Trends Pharmacol. Sci.
38
,
322
-
338
.
McCarron
,
J. G.
,
Wilson
,
C.
,
Heathcote
,
H. R.
,
Zhang
,
X.
,
Buckley
,
C.
and
Lee
,
M. D.
(
2019
).
Heterogeneity and emergent behaviour in the vascular endothelium
.
Curr. Opin. Pharmacol
45
,
23
-
32
.
Millán
,
J.
,
Cain
,
R. J.
,
Reglero-Real
,
N.
,
Bigarella
,
C.
,
Marcos-Ramiro
,
B.
,
Fernández-Martín
,
L.
,
Correas
,
I.
and
Ridley
,
A. J.
(
2010
).
Adherens junctions connect stress fibres between adjacent endothelial cells
.
BMC Biol.
8
,
11
.
Park-Windhol
,
C.
and
D'Amore
,
P. A.
(
2016
).
Disorders of vascular permeability
.
Annu. Rev. Pathol
11
,
251
-
281
.
Pedrosa
,
A. R.
,
Trindade
,
A.
,
Fernandes
,
A. C.
,
Carvalho
,
C.
,
Gigante
,
J.
,
Tavares
,
A. T.
,
Diéguez-Hurtado
,
R.
,
Yagita
,
H.
,
Adams
,
R. H.
and
Duarte
,
A.
(
2015
).
Endothelial Jagged1 antagonizes Dll4 regulation of endothelial branching and promotes vascular maturation downstream of Dll4/Notch1
.
Arterioscler. Thromb. Vasc. Biol
.
35
,
1134
-
1146
.
Potente
,
M.
,
Gerhardt
,
H.
and
Carmeliet
,
P.
(
2011
).
Basic and therapeutic aspects of angiogenesis
.
Cell
146
,
873
-
887
.
Rafii
,
S.
,
Butler
,
J. M.
and
Ding
,
B.-S.
(
2016
).
Angiocrine functions of organ-specific endothelial cells
.
Nature
529
,
316
-
325
.
Roukos
,
V.
,
Pegoraro
,
G.
,
Voss
,
T. C.
and
Misteli
,
T.
(
2015
).
Cell cycle staging of individual cells by fluorescence microscopy
.
Nat. Protoc
10
,
334
-
348
.
Schindelin
,
J.
,
Arganda-Carreras
,
I.
,
Frise
,
E.
,
Kaynig
,
V.
,
Longair
,
M.
,
Pietzsch
,
T.
,
Preibisch
,
S.
,
Rueden
,
C.
,
Saalfeld
,
S.
,
Schmid
,
B.
et al. 
(
2012
).
Fiji: an open-source platform for biological-image analysis
.
Nat. Methods
9
,
676
-
682
.
Seebach
,
J.
,
Cao
,
J.
and
Schnittler
,
H. J.
(
2016
).
Quantitative dynamics of VE-cadherin at endothelial cell junctions at a glance: basic requirements and current concepts
.
Discov. Craiova Rom.
4
,
e63
.
Seebach
,
J.
,
Klusmeier
,
N.
and
Schnittler
,
H.
(
2020
).
Autoregulatory “multitasking” at endothelial cell junctions by junction-associated intermittent lamellipodia controls barrier properties
.
Front. Physiol.
11
,
586921
.
Seebach
,
J.
,
Taha
,
A. A.
,
Lenk
,
J.
,
Lindemann
,
N.
,
Jiang
,
X.
,
Brinkmann
,
K.
,
Bogdan
,
S.
and
Schnittler
,
H.-J.
(
2015
).
The CellBorderTracker, a novel tool to quantitatively analyze spatiotemporal endothelial junction dynamics at the subcellular level
.
Histochem. Cell Biol.
144
,
517
-
532
.
Siu
,
D. M. D.
,
Lee
,
K. C. M.
,
Lo
,
M. C. K.
,
Stassen
,
S. V.
,
Wang
,
M.
,
Zhang
,
I. Z. Q.
,
So
,
H. K. H.
,
Chan
,
G. C. F.
,
Cheah
,
K. S. E.
,
Wong
,
K. K. Y.
et al. 
(
2020
).
Deep-learning-assisted biophysical imaging cytometry at massive throughput delineates cell population heterogeneity
.
Lab. Chip
20
,
3696
-
3708
.
Stepanova
,
D.
,
Byrne
,
H. M.
,
Maini
,
P. K.
and
Alarcón
,
T.
(
2021
).
A multiscale model of complex endothelial cell dynamics in early angiogenesis
.
PLoS Comput. Biol.
17
,
e1008055
.
Tikhonova
,
A. N.
,
Dolgalev
,
I.
,
Hu
,
H.
,
Sivaraj
,
K. K.
,
Hoxha
,
E.
,
Cuesta-Domínguez
,
Á.
,
Pinho
,
S.
,
Akhmetzyanova
,
I.
,
Gao
,
J.
,
Witkowski
,
M.
et al. 
(
2019
).
The bone marrow microenvironment at single-cell resolution
.
Nature
569
,
222
-
228
.
Torres-Vázquez
,
J.
,
Kamei
,
M.
and
Weinstein
,
B. M.
(
2003
).
Molecular distinction between arteries and veins
.
Cell Tissue Res.
314
,
43
-
59
.
Travaglini
,
K. J.
,
Nabhan
,
A. N.
,
Penland
,
L.
,
Sinha
,
R.
,
Gillich
,
A.
,
Sit
,
R. V.
,
Chang
,
S.
,
Conley
,
S. D.
,
Mori
,
Y.
,
Seita
,
J.
et al. 
(
2020
).
A molecular cell atlas of the human lung from single-cell RNA sequencing
.
Nature
587
,
619
-
625
.
Ubezio
,
B.
,
Blanco
,
R. A.
,
Geudens
,
I.
,
Stanchi
,
F.
,
Mathivet
,
T.
,
Jones
,
M. L.
,
Ragab
,
A.
,
Bentley
,
K.
and
Gerhardt
,
H.
(
2016
).
Synchronization of endothelial Dll4-Notch dynamics switch blood vessels from branching to expansion
.
eLife
5
,
e12167
.
Veschini
,
L.
,
Belloni
,
D.
,
Foglieni
,
C.
,
Cangi
,
M. G.
,
Ferrarini
,
M.
,
Caligaris-Cappio
,
F.
and
Ferrero
,
E.
(
2007
).
Hypoxia-inducible transcription factor–1 alpha determines sensitivity of endothelial cells to the proteosome inhibitor bortezomib
.
Blood
109
,
2565
-
2570
.
Veschini
,
L.
,
Crippa
,
L.
,
Dondossola
,
E.
,
Doglioni
,
C.
,
Corti
,
A.
and
Ferrero
,
E.
(
2011
).
The vasostatin–1 fragment of chromogranin A preserves a quiescent phenotype in hypoxia–driven endothelial cells and regulates tumor neovascularization
.
FASEB J.
25
,
3906
-
3914
.
Vestweber
,
D.
,
Winderlich
,
M.
,
Cagna
,
G.
and
Nottebaum
,
A. F.
(
2009
).
Cell adhesion dynamics at endothelial junctions: VE-cadherin as a major player
.
Trends Cell Biol.
19
,
8
-
15
.
Wiseman
,
E.
,
Zamuner
,
A.
,
Tang
,
Z.
,
Rogers
,
J.
,
Munir
,
S.
,
Di Silvio
,
L.
,
Danovi
,
D.
and
Veschini
,
L.
(
2019
).
Integrated multiparametric high-content profiling of endothelial cells
.
SLAS Discov. Adv. Life Sci. R D
24
,
264
-
273
.
Yoder
,
M. C.
(
2018
).
Endothelial stem and progenitor cells (stem cells): (2017 Grover Conference Series)
.
Pulm. Circ
8
,
204589321774395
.
Zong
,
D.
,
Li
,
J.
,
Cai
,
S.
,
He
,
S.
,
Liu
,
Q.
,
Jiang
,
J.
,
Chen
,
S.
,
Long
,
Y.
,
Chen
,
Y.
,
Chen
,
P.
et al. 
(
2018
).
Notch1 regulates endothelial apoptosis via the ERK pathway in chronic obstructive pulmonary disease
.
Am. J. Physiol. Cell Physiol
315
,
C330
-
C340
.

Competing interests

D.D. is an employee of King's College London and an employee of bit.bio. D.D. declares no other affiliations with or involvement in any organization or entity with any financial or non-financial interest in the subject matter or materials discussed in the manuscript. All other authors declare no competing or financial interests.

Supplementary information