Data driven and biophysical insights into the regulation of trafficking vesicles by extracellular matrix stiffness.

Biomechanical signals from remodeled extracellular matrix (ECM) promote tumor progression. Here, we show that cell-matrix and cell-cell communication may be inherently linked and tuned through mechanisms of mechanosensitive biogenesis of trafficking vesicles. Pan-cancer analysis of cancer cells' mechanical properties (focusing primarily on cell stiffness) on substrates of varied stiffness and composition elucidated a heterogeneous cellular response to mechanical stimuli. Through machine learning, we identified a fingerprint of cytoskeleton-related proteins that accurately characterize cell stiffness in different ECM conditions. Expression of their respective genes correlates with patient prognosis across different tumor types. The levels of selected cytoskeleton proteins indicated that cortical tension mirrors the increase (or decrease) in cell stiffness with a change in ECM stiffness. A mechanistic biophysical model shows that the tendency for curvature generation by curvature-inducing proteins has an ultrasensitive dependence on cortical tension. This study thus highlights the effect of ECM stiffness, mediated by cortical tension, in modulating vesicle biogenesis.

Introduction

Cellular functioning, normal or pathological, is orchestrated by interactions with other cells, extracellular matrix (ECM), and other components such as extracellular vesicles in the microenvironment (Bissell and Hines, 2011). Besides providing structural support to tissues and organs, the ECM influences basic cellular processes ranging from differentiation and apoptosis to migration and proliferation (Hynes, 2009). A microenvironment more conducive to tumorigenesis and metastasis is created by tumors leveraging dysregulated ECM remodeling mechanisms (Bonnans et al., 2014; Winkler et al., 2020). One of the critical mechanical alterations is the stiffening of ECM in tumors (Masuzaki et al., 2007; Levental et al., 2009; Laklai et al., 2016). Cellular mechanosensing via transmembrane receptors (such as integrins) of a stiff matrix enhances oncogenic signaling pathways promoting growth, survival, and invasion leading to tumor progression (Kai et al., 2016; Miroshnikova et al., 2017; Northey et al., 2017). A stiff ECM also promotes myeloid cell polarization and alters T cell function, broadening its effects to regulate tumor immunity (Springer and Fischbach, 2016; Meli et al., 2019; Chirivì et al., 2021).

Cells respond to mechanical stimuli (such as a change in substrate stiffness) by altering their mechanics such as cell stiffness, motility, adhesion, and contractility (Janmey et al., 2020). The mechanical properties of cells have been shown to be modulated by the organization of the cytoskeleton (Fletcher and Mullins, 2010). Studies have also demonstrated mechanosensing-based remodeling of the cytoskeleton (Solon et al., 2007; Trichet et al., 2012; Gupta et al., 2015). Although dynamic reorganization of the cytoskeleton is essential in the cell cycle, morphogenesis, and cell motility, these processes often go awry in cancer. Deregulated cytoskeleton dynamics have been implicated in tumor cell migration, invasion, and metastasis (Yamaguchi and Condeelis, 2007; Hall, 2009; Sun et al., 2015).

Communication between a cell and its surroundings is mediated through biochemical and physical interactions of the plasma membrane with the microenvironment constituents. Endocytosis, involving the internalization of the plasma membrane (lipids and proteins), extracellular material, and subsequent cargo trafficking, contributes to several cellular processes such as cell adhesion and migration, signal transduction, and cell polarity (Grant and Donaldson, 2009; Doherty and McMahon, 2009). Endosomal cargo can either be recycled to the membrane or culminate in lysosomes for degradation or secreted extracellularly through the release of exosomes. Studies have reported aberrations in endocytic pathways of cancer cells (Mellman and Yarden, 2013). For instance, enhanced integrin recycling promotes invasive migration of cancer cells (Caswell and Norman, 2008; Muller et al., 2009), mechanisms that undermine growth factor receptors (such as EGFR and MET) degradation via lysosomes enhance proliferation (Mosesson et al., 2008), and deregulated E-cadherin trafficking disrupts cell polarity and induces epithelial-to-mesenchymal transition (Paredes et al., 2012). The regulatory machinery of endocytic trafficking being “hijacked” in tumors is further highlighted by the expression of chemokine receptors on tumor cells being crucial in metastasis (Zlotnik et al., 2011; Marchese, 2014; Chow and Luster, 2014).

Among various extracellular vesicles in the microenvironment, exosomes are nano-sized lipid vesicles produced through the endocytic pathway and secreted via exocytosis (Edgar, 2016; Gurung et al., 2021). Exosomes are efficient mediators of cell-to-cell communication, with the exosomal cargo of parent cell-specific bioactive molecules affecting the function and behavior of recipient cells (Gurung et al., 2021). Tumor-derived exosomes (TEXs) have been implicated in cancer progression and metastasis (Whiteside, 2016), with studies showing their role in promoting tumor cell invasion (Hoshino et al., 2013) and migration (Sung et al., 2015) as well as helping establish pre-metastatic niches at distant tissue sites (Peinado et al., 2017). TEXs are also a vital component of the ability of cancer cells to evade the immune system (Shinohara et al., 2017; Chen et al., 2018; Kurywchak et al., 2018).

From clathrin-mediated endocytosis at the plasma membrane to the biogenesis of exosomes as intraluminal vesicles through inwards budding in late endosomes, membrane curvature remodeling involving highly curved structures (such as vesicles and tubules) is an essential component of trafficking pathways. Among several mechanisms postulated for membrane remodeling (McMahon and Boucrot, 2015), a class of peripheral membrane proteins capable of generating curvature plays a central role. The importance of epsin N-terminal homology (ENTH) domain in inducing curvature in clathrin-coated pits during endocytosis has been well established (Ford et al., 2002; Agrawal et al., 2010). Bin-Amphiphysin-Rvs (BAR) domains have also been reported to play a primary role in high membrane curvature generation and stabilization in endocytic machinery (Dawson et al., 2006). Endosomal Sorting Complex Required for Transport (ESCRT) III protein complexes are considered to be primary inducers of curvature necessary for vesiculation in ILV formation (Henne et al., 2013). Curvature generation by these proteins will depend on the membrane’s mechanical properties (such as tension and bending rigidity). Forces from the cytoskeleton also strongly influence membrane deformability by modulating the effective tension in the membrane through membrane-to-cortex attachments (Diz-Muñoz et al., 2013; Sens and Plastino, 2015). Cytoskeleton-related tension will, therefore, affect the curvature-generating capabilities of the curvature-inducing proteins.

The abnormally altered ECM stiffness observed in tumorigenesis and the dependence of cytoskeletal dynamics on matrix elasticity suggest a mechanosensitive mechanism employed by tumors for regulating trafficking pathways. We thus asked whether stiff ECM affects the membrane curvature generation needed in trafficking pathways through its effect on cortical tension. To answer this question, in this study, we analyzed a pan-cancer mechanobiology dataset to elucidate the dependence of cell mechanics of different cancer cells on ECM stiffness and composition. Next, we used information-theoretic and machine learning-based approaches to identify key cytoskeletal proteins mediating the effect of ECM stiffness on cell mechanics. Survival analysis was performed to ascertain the importance of selected proteins from a clinical perspective. We then employed a continuum mesoscale membrane model to delineate how membrane deformations (such as vesiculation and tubulation) orchestrated by curvature-inducing proteins depend on tension owing to cytoskeleton. This connection between cortical tension and vesicle biogenesis confers an intriguing paradigm of mechanosensitivity to cellular functions and further underscores the need for considering biophysical effects on signaling pathways, to gain a holistic understanding of how mechanical stimuli such as ECM stiffening aid in tumor initiation and progression.

Results

Characterization of cell stiffness in different cancer cells across various substrates

The effects of changes in the mechanical properties of tissues are increasingly implicated in the normal development of tissues and the neoplastic functions of cells in various disease states (Discher et al., 2005). One classic example of this process is the stiffening of liver tissue in the early stages of fibrotic disease, which is caused by increasing stiffness of the extracellular matrix (ECM) well before changes in cell morphology or ECM deposition is observed (Georges et al., 2007). The effect of substrate stiffness and composition, therefore, might be expected to affect many aspects of both normal and cancer cell function. The National Cancer Institute (NCI) Leidos mechanobiology dataset is compiled from high-throughput studies on the mechanical properties of cancer cells in various substrates mimicking the ECM. The dataset has cell mechanics-related data (such as cell stiffness and motility) for 20 cancer cell lines grown on seven different ECM substrates. Based on the cell stiffness (Young’s modulus, ) data in the Leidos dataset (Figure 1A), we characterize how cancer cell stiffness changes in different ECM substrates of varying stiffness for different cell types. For our analysis, we computed the kernel density estimates (KDE) for summarizing the distribution of Young’s modulus: across all cell types and all substrates (KDEall) and for each cell type across all substrates (Figures S1 and S2).

Figure 1

Cell stiffness characterization for cancer cell lines grown on different ECM substrates

(A) 20 cancer cell lines grown on seven different matrix conditions explored in cell stiffness measurements. (HA: Hyaluronic Acid).

(B) Cell stiffness (Pa) for each cell line averaged over respective cells grown on different ECM substrates (left). Ratio of cell counts in high vs. low stiffness categories where cutoff used for demarcation is the median of KDEall (right).

(C and D) COL: Collagen, FN: Fibronectin, HACOL: Collagen coated HA, and HAFN: Fibronectin coated HA.

(C) Average cell stiffness (Pa) for each cell line on each of the seven ECM substrates.

(D) Ratio of cell counts in high vs. low stiffness categories where cutoff used for demarcation is the median of the cell type-specific KDE of cell stiffness over all substrates.

(E) For each cell line, the average value of the ratio of the Young’s Moduli of the cells in moving from one substrate to another as a measure of the cell stiffness sensitivity to substrate change. (COL fold: 30kPa COL/500Pa COL, FN fold: 30kPa FN/500Pa FN, HYAL fold: HACOL/HAFN, Glass fold: Glass/500Pa COL, HAFN fold: HAFN/500Pa FN, HACOL fold: HACOL/500Pa COL).

First, we examined how cell stiffness varies across different cell types. For each cell type, the average value of across the seven different substrates was computed (Figure 1B). Although skin-derived cells did not show significant variation in stiffness across different cell lines, other tissue types had significant variation in cell stiffness among their respective cell lines, with breast cancer cell lines showing the most variation (3.5 kPa–11 kPa). We then set out to quantify the relative counts of low and high stiffness cells within each cell type using the median of KDEall to split the cells into low and high categories (Figure 1B). Some interesting results can be seen. A375, WM266-4 (skin cancer cell lines), SW480 (colon cancer cell line), and MCF-10A-JSB (breast cancer cell line) have a very high proportion of their cells with high cell stiffness across all ECM substrate conditions. The remaining skin, colon, and breast cell types (except SW620 and T-47D) also show a higher proportion of high stiffness cells. On the other hand, pancreatic, prostate, and lung cancer cell lines have a more significant proportion of low stiffness cells. Overall, the results suggest a variation in the distribution of low stiffness and high stiffness cancer cells within tissue types and across tissue types.

In order to gain a more granular understanding of how cell stiffness varies within a cell type on a particular ECM substrate, we computed the mean of the cell stiffness for a cell line on each of the seven ECM substrates (Figure 1C). An interesting picture emerges from Figure 1C, which shows a variation in cell stiffness within tissue-specific cell types in different substrate conditions and across tissue types. The count distribution of the cell type and substrate-specific ratio of high stiffness to low stiffness cells is shown in Figure 1D. Unlike before, splitting cells of a cell type on a particular ECM substrate into low and high categories is based on the median of that cell-type-specific KDE of cell stiffness over all substrates. This bucketing was conducted because we were interested in comparing the behavior of a cell line on a specific substrate with the average behavior of that cell line over all substrates. Some interesting observations stand out. Breast, colon, and skin cancer cell lines have a higher proportion of high stiffness cells to low stiffness cells in 30 kPa collagen as compared to 500 Pa collagen (a softer ECM). However, not all of the breast, colon, and skin cell lines show a similar trend on fibronectin substrates. Lung and pancreatic tissue cell types show significantly greater proportions of high stiffness cells on 500 Pa fibronectin as compared to 30 kPa fibronectin. This analysis highlights that different cell types have varying dependence of cell stiffness on ECM stiffness and composition.

The cell stiffness sensitivity to ECM substrate (which is essential in cell mechanosensing) was further explored by analyzing the relative change in cell stiffness when moving from one substrate to another. Six different substrate sensitivity metrics were considered, namely COL : 30kPa COL/500Pa COL, FN fold: 30kPa Fn/500Pa Fn, HYAL fold: HACOL/HAFN, Glass fold: Glass/500Pa COL, HAFN fold: HAFN/500PA FN, HACOL fold: HACOL/500PA COL. The average ratio of cell stiffness for each cell line was calculated in each of the six substrate sensitivity regimes (Figure 1E).

We see a wide variation in the relative change in cell stiffness as a function of substrate sensitivity. For instance, Panc-1 (pancreatic cancer cell line) and lung cancer cell lines show a decrease in cell stiffness with increasing ECM stiffness in both collagen and fibronectin matrices. This is in contrast to the breast and colon tissue cell types which primarily have higher cell stiffness when grown on more stiff substrate. However, for an increase in fibronectin matrix stiffness (FN fold) HCT116 (colon cancer cell line) and MCF10A-JSB (breast cancer cell line) show a decrease in cell stiffness. Overall, the results show the heterogeneity in how cells within and across tissue-specific cancer cell lines respond to mechanical cues such as increasing ECM stiffness.

Essential cytoskeleton proteins mediating effects of extracellular matrix stiffness on cell stiffness

Having characterized the cell stiffness of various cancer cells in various ECM substrates of varying stiffness and composition, we next asked two questions. (1) Which proteins seem to be implicated in the ability of cells to respond to mechanical cues in the context of modulating their cell stiffness? Moreover, (2) is there a defined expression signature of proteins that can capture the cell stiffness variations in response to changes in ECM mechanics? To gain critical insights related to these questions, we used the genomics and proteomics data available in the NCI Leidos mechanobiology dataset (in addition to cell mechanics-related data) for nine cell lines (A375, RWPE-1, 22RV1, SW620, SW480, T-47D, MDA-MB-231, U-87, and T98G) grown on seven different ECM substrates. The richness of the multi-omics data at single-cell resolution provides unique new opportunities to construct predictive data-driven models to relate genotype to cell states to phenotype across different microenvironment mechanics and chemistry.

In order to understand which proteins could be implicated in cell response to mechanosensing, mRNA diff-seq analysis of the Leidos genomics dataset to analyze gene expression profiles of cancer cells mechanosensing of different substrate microenvironments was performed, which yielded ∼800 genes, of which ∼180 were reflected in the proteomics dataset. Gene ontology (GO) analysis of these 180 genes showed a preponderance of several gene categories: signaling, extracellular region, ECM organization, actin-binding, cell adhesion, and cell junction. We focused on the role of cytoskeletal-related proteins in modulating cell stiffness in response to ECM mechanics.

The cytoskeleton comprises hundreds of diverse signaling and scaffolding proteins that modulate cell shape, migration, and cellular mechanosensing (Fletcher and Mullins, 2010). A manually curated list of cytoskeletal-related proteins was compiled. These proteins include the major components of the cytoskeleton itself: microfilaments (actin), intermediate filaments (vimentin, keratin, and lamin), and microtubules (tubulin). We surmised that changes to the quantities of these proteins should correlate strongly with significant morphological changes to the cell and its stiffness. Accessory proteins, including the actin cross-linkers filamin alpha-actinin, nesprin, and plectin, were included in the list. Finally, the list included proteins responsible for cell adhesion and force transduction, namely, integrin adhesion receptors, components talin and vinculin, and myosin motor protein. In order to couple the protein expression levels of these proteins to cell stiffness, we applied an information-theoretic measure of mutual information (I normalized to lie between 0 and 1). Proteins from the curated list whose I with cell stiffness ≥0.8 were selected, producing a resultant list of 18 proteins (Figure 2A).

Figure 2

Cytoskeleton proteins correlated with and predictive of cell stiffness for cancer cell lines grown on different ECM substrates

(A) 18 cytoskeletal-related proteins whose mutual information with cell stiffness ≥0.8.

(B) PCA heatmap of first five principal components (PC1-5) of protein expression levels in low stiffness and high stiffness cells. Pairwise angles between the respective PC1-5 of low and high stiffness cells.

(C) Shapley values of the selected 18 cytoskeletal-related proteins used as features in neural network to predict cell stiffness category (with balanced accuracy of 98% on the test set). Histogram showing the distribution of protein expressions of the top five proteins in low and high cell stiffness categories.

Cells were split into low and high cell stiffness regimes with demarcation based on the median of the KDE of cell stiffness across all nine cell types and all substrates in order to streamline further analysis. We applied principal component analysis (PCA) to the protein count of the selected 18 proteins in the low and high cell stiffness datasets and focused on the first five principal components (PC1-5), which accounted for approximately 90% of the variation in the data (Figure 2B). There is a clear distinction between low stiffness and high stiffness regime cells, as observed from the heatmaps. We computed the pairwise angles between PC1-5 of the low stiffness and high stiffness cells to quantify their difference, which clearly showed that the protein expression signature significantly differed (Figure 2B). Based on these insights, we explored the predictive power of these 18 proteins in predicting cell stiffness (low/high stiffness). A neural network predicted with balanced accuracy (BACC) of 98% on a test dataset if a cell was likely to be belonging to the low stiffness or high stiffness category.

A model explainability metric based on treating the features (protein counts of the proteins in Figure 2A) as players in a cooperative game and applying a game-theoretic measure called the Shapley value (Lundberg and Lee, 2017) was applied to the neural network (Figure 2C). This gives a measure of how much each feature contributes to the prediction. The top five contributors to the prediction of cell stiffness in our neural network are alpha-actinin-4 (ACTN4), myosin-9 (MYH9), plectin (PLEC), talin-1 (TLN1), and vinculin (VCL). The histogram of protein counts of these proteins shows that their protein expression signatures are generally higher in high stiffness cells than in low stiffness cells (Figure 2C). We also trained a neural network having the same architecture as before but input features consisting only of Shapley-based top five proteins which gave a high balanced accuracy of 85% on the test set.

This section gives an intriguing insight into how the highly diverse behavior among different cancer cell lines in terms of their cell stiffness on varied ECM substrates (Figure 1) can be captured through a low-dimensional fingerprint involving certain essential cytoskeletal proteins. Altogether, the results suggest that the 18 selected cytoskeleton-related proteins are implicated in modulating cell stiffness in response to ECM mechanical cues. Furthermore, considering that the selected proteins contain several actin cross-linking proteins (alpha-actinin and filamins), cytoskeletal linker protein (plectin) and myosin-II subunit (myosin-9) that have good correlative and predictive power with respect to cell stiffness of different cell types on different ECM substrates, the analysis indicates that the underlying cortical stiffness and tension significantly modulate the cell stiffness.

Clinical significance of the selected cytoskeleton genes

To further clarify the significance of the selected cytoskeleton genes, we performed a survival analysis using the clinical patient data for different cancer types from the TGCA PanCancer Atlas dataset on the cBioPortal platform (http://cbioportal.org) (Cerami et al., 2012; Gao et al., 2013). The patient dataset for each tumor type was divided into three categories based on the normalized mRNA expression () of a gene. In brief, the -score of a gene refers to how many standard deviations away is the gene expression for a patient as compared to the mean gene expression in the reference population. The reference population used was all profiled tumor samples of a particular cancer type. For each gene, the three categories are then defined as (i) samples having low expression (-)), (ii) unaltered samples (-), and (iii) samples having high expression (). Following this, the clinical data of the patients in each of the three groups were used for Kaplan-Meier (KM) analysis, wherein we considered the overall survival and the progression-free survival curves and compared unaltered and low/high categories.

Table 1 summarizes the results from the survival analysis for different cancer types in terms of median months of overall survival and progression-free survival. Note that, for both overall and progression-free survival, the difference between the KM curves of unaltered and low/high samples was considered statistically significant for logrank test p-value of less than 0.05. Additionally, only those KM plots that had enough samples (and events) for computing the 95% confidence interval for median months were taken into account. An interesting spread of median months of overall survival and progression-free survival was observed. Skin cutaneous melanoma patients with higher expression of some of the selected cytoskeleton genes had lower median survival. For ovarian serous cystadenocarcinoma and lung adenocarcinoma, there are also cases where low expression of a particular gene leads to lower median survival. The effect of a particular gene varied across tissue types. For instance, low LMNA expression patients had lower and higher median overall survival for ovary and brain cancers, respectively. Overall, the survival analysis suggests that the expression of the essential cytoskeletal-related genes identified in this study could have a significant effect on the clinical outcome of a patient across different tumor types.

Table 1

Median months (for overall and progression-free survival) from survival analysis performed using cBioPortal on patient-derived clinical data from TCGA PanCancer Atlas dataset

Tissue Type	Ovary			Skin			Brain			Lung
	Low	Unaltered	High	Low	Unaltered	High	Low	Unaltered	High	Low	Unaltered	High
ACTN4		44.84a,c	29.03a,c
MYH9					37.38a,d	19.50a,d
PLEC	12.49a,c	18.90a,c						14.73a,c	10.82a,c		53.33b,c	29.75b,c
											42.51b,c	23.80b,c
TLN1					94.98a,c	58.52a,c
VCL											46.72a,c	37.71a,c
ACTN1	16.83a,c	18.51a,c
LMNA	34.03a,c	45.01a,c			94.98a,c	61.28a,c	19.82a,c	13.61a,c
FLNB		44.32a,c	41.06a,c		93.01b,d	53.19b,d
	14.43a,d	18.90a,c,d	16.83a,c
LMNB2					40.47b,c	22.32b,c
ACTN3		44.84a,c	33.67a,c							44.61a,c	53.65a,c
FLNC											54.34b,c	21.57b,c
KRT18											53.33b,c	29.46b,c
ACTG1					94.98a,c	55.59a,c
GSN								14.73a,c	7.36a,c
								7.59a,c	5.79a,c
FLNA					94.98a,c	55.59a,c

Low and high correspond to samples having normalized mRNA expression ( −) and , respectively. Unaltered refers to the profiled tumor samples having gene expression within “threshold” standard deviations of the mean expression of all profiled tumor samples of a particular cancer type. Median months for progression-free survival are shown in italic. See also Figures S5–S11 for the KM plots. Ovary: ovarian serous cystadenocarcinoma, Skin: skin cutaneous melanoma, Brain: glioblastoma multiforma, Lung: lung adenocarcinoma. See Table S1 for details about the TCGA cancer datasets. Note that only cases that had logrank test p-value < 0.05, and also had sufficient events to establish a 95% confidence interval for median months have been reported in the table.

logrank test p-value < 0.05.

logrank test p-value < 0.01 for statistically significant difference between survival curves of unaltered and low/high samples.

z-score threshold = 1.0.

z-score threshold = 1.5.

Ultrasensitive response of membrane deformation by curvature-inducing proteins to changes in cortical tension

Curvature-inducing proteins are part of a specialized class of cytosolic proteins that can sense, induce and stabilize high curvature local morphologies (like tubes and buds) of membranes (Shibata et al., 2009; Kozlov et al., 2014). One of the key features in the curvature generation mechanism observed for these proteins is the spatial localization of functional units via oligomerization and via membrane-mediated forces for inducing sufficient curvature followed by self-assembly of these functional units on the membrane to generate macroscopic high curvature morphologies such as vesiculation and tubulation (Yu and Schulten, 2013; Simunovic et al., 2013; Zhao et al., 2013; Tourdot et al., 2014; Daum et al., 2016). A previous study (Tourdot et al., 2015) proposed a micellization-like model for membrane tubulation, wherein tubes formed by self-assembly of curvature-inducing proteins are analogous to micelle aggregates formed by self-assembly of surfactant molecules. Here, we show that this self-assembly mechanism is crucial to explain the tension dependence of vesiculation.

In cells, the deformation of a patch of membrane is also dependent on the cortical actin tension owing to the associated cortex (Diz-Muñoz et al., 2013; Granger et al., 2014). To address how the curvature-generating capabilities of curvature-inducing proteins is affected by cortical tension, we used dynamically triangulated Monte Carlo (DTMC) simulations to explore the morphological conformational space of planar to highly curved membranes (Ramakrishnan et al., 2014; Tourdot et al., 2014, 2015). We utilized the DTMC framework to simulate a membrane patch and model the surface deformations induced by curvature-inducing proteins, which are represented as curvature fields (each field represents 10 protein units). The dependence of cortical tension is incorporated by simulating membrane patches with different excess areas (A/A, A is the curvilinear area of the membrane patch, and is the projected area) where increasing excess area corresponds to decreasing cortical tension (Ramakrishnan et al., 2018).

As our focus in this study is highly curved structures produced by morphological transitions (such as observed in endocytic events), we use the universality that both vesiculation and tubulation can be described in terms of the self-assembly of curvature-inducing proteins to generate a morphological change in the membrane (Tourdot et al., 2014, 2015). Using the DTMC framework, we analyzed how the membrane morphology changes with increasing excess area for a given total number of protein-induced curvature fields (n) (Figure 3A). The extent of curvature change can be seen in terms of how many protein fields are in the tubular region (n) for the equilibrium membrane morphologies. The equation for describing the simulation data (n vs A/A) was obtained by generalizing the micelle model proposed by Tourdot et al. (2015) as described in STAR Methods (also see Figure S3) and is given by:where a, b, and c are constants for a given n. The parameter n can be viewed as a Hill-type coefficient indicating cooperativity between protein fields to induce tubulation. For a fixed value of n, Equation 1 was fitted to the data (n vs A/A) obtained from simulations for different n values. The heatmap in Figure 3B summarizes the mean square errors (MSE) for different values of across various nₚ used in simulations. The MSE values for were much higher than those for n ≤ 7. Among the values of 2–7, n = 4 had lower MSE values compared to n = {2,3} for n ≤ 16 and n = {5,6,7} for n > 16. Furthermore, n = 4 is consistent with the number of protein fields per tube being approximately four post onset of tubulation in (Tourdot et al., 2015). Therefore, Equation 1 with n = 4 was chosen as the representative of n as a function of A/A. Figure 3C shows the ultrasensitive dependence of on excess area (A/A) for different values of n.

Figure 3

Effect of cortical tension on curvature generation by curvature-inducing proteins

(A) Snapshot illustrating the change in membrane curvature owing to curvature-inducing proteins (here, modeled as curvature fields).

(B) Mean square errors (MSE) between analytical (non-linear fit of Equation 1) and simulation data for different values of degree n in Equation 1 across various values of the total number of protein fields (n) on the simulated membrane patch.

(C) Number of protein fields in the tubular region (n) as a function of excess area (A/A, proxy for cortical tension) for different values of n on the simulated membrane patch (Eqn: non-linear fits of Equation 1 with n = 4, Sim: simulation data).

The number of proteins on tubular region (n) plateaus with increasing excess area indicating that the extent of morphological change in a membrane patch induced by curvature-generating proteins eventually saturates with lowering of cortical tension. Parameter a′ in Equation 1 defines the saturation value of n and was found to increased linearly with n (Figure S4). Parameters ′b′ and ′c′ in Equation 1 characterize how fast the saturation is achieved, and their values monotonically increased and decreased, respectively, as a function of n (Figure S4). This highlights the observation that with more curvature-inducing proteins on the membrane (higher n), protrusion formation can happen at higher cortical tension (Figure 3C).

It is important to note the micellization-based model used as the basis for Equation 1 which underscores an interesting proposition wherein the ability of curvature-inducing proteins to self-assemble for the formation of membrane protrusions is modulated by cortical tension. Hence, the Hill function-like dependence of curvature generation on the excess area (proxy for cortical tension) suggests that lower cortical tension increases the tendency for curvature-inducing proteins to congregate and thereby, promotes morphological changes such as tubulation and vesiculation.

Discussion

The role of a stiff ECM in promoting tumor progression has been well established (Levental et al., 2009; Kai et al., 2016; Miroshnikova et al., 2017). To investigate the response of cancer cells to mechanical cues such as ECM stiffening, we used the cell stiffness data from the pan-cancer NCI Leidos Mechanobiology Dataset. Our analysis suggests that cancer cells of the same cell type respond to mechanical stimuli (as quantified by ECM substrate stiffness) in varied ways and that while mechanosensing mechanisms are functions of cell types, even within a given cell type, there is a clear distribution of emergent behavior. This spread is consistent with reported heterogeneity in the stiffness of cells both within and between breast cancer cell lines (Shen et al., 2020) and substrate-dependent sensitivity of cell stiffness for melanoma cell lines (Lekka et al., 2012).

Using information theory and machine learning-based analysis of the Leidos proteomics dataset, we identified a list of 18 cytoskeletal-related proteins which might be critical potentiators of the regulation of cell stiffness in response to mechano-stimuli such as ECM stiffness. The survival analysis further highlighted how the expression of respective genes of the selected cytoskeleton proteins correlates with the prognosis of patients across different tumor types. Additionally, there are numerous independent measurements that support the importance of the proteins identified here. For example, a specific C-terminal actin-binding site of talin is required for coupling the membrane to the cytoskeleton to regulate membrane tension (Schulte et al., 2016), and a separate site on talin mediates integrin-dependent regulation of membrane trafficking (Margadant et al., 2012). Mechanical unfolding of filamin is important in the tension-induced autophagy pathway of myocytes and immune cells (Ulbricht et al., 2013), and filamin is required for cytoskeletal responses to tension applied at the cell membrane that leads to altered Ca2+ influx (Glogauer et al., 1998). Recruitment of alpha actinin by myosin one- generated tension at the cell membrane is important for strengthening cell-cell cohesion in epithelial cell monolayers (Kannan and Tang, 2018). Vinculin is required for the membrane tension dependence of the dynamin-independent CLIC/GEEC (CG) endocytic pathway (Thottacherry et al., 2018).

Among the top three proteins based on Shapley analysis, the higher expression of the crucial actin-actin and actin-integrin cross-linking protein (ACTN4) and the versatile cytoskeletal linker protein (PLEC) point toward an increased cortical stiffness (Wiche, 1998; Ehrlicher et al., 2015), while higher expression of the myosin-II subunit (MYH9) is indicative of increased cortical actin tension (Clark et al., 2007). This finding suggests that an increase (or decrease) in cell stiffness with increasing ECM stiffness could be seen as an increase (or decrease) in cortical stiffness and tension. In our analysis of cell stiffness data, the increase in cell stiffness observed for several cancer cell lines with increasing ECM stiffness can, therefore, be explained by the upregulation of RhoA/ROCK pathway leading to increased actomyosin tension (Paszek et al., 2005; Rath and Olson, 2012). In contrast, some cancer cell lines showed a reduction in cell stiffness with increasing ECM stiffness, which could be potentiated by MYC (a proto-oncogene (Venkateswaran and Conacci-Sorrell, 2017)) through a feedback loop that negatively regulates specific RhoA downstream signaling branches connected to the induction of stress fibers, focal adhesions and actomyosin contractility (Sauzeau et al., 2010).

Altered components of trafficking pathways (Grant and Donaldson, 2009) have been implicated in promoting tumor invasion and proliferation of cancer cells (Mellman and Yarden, 2013). We asked if membrane curvature generation (critical in endocytic pathways) is regulated by ECM stiffness through its effect on cortical tension. Based on purely biophysical considerations, our results suggest an ultrasensitive response with respect to cortical tension in membrane deformation induced by curvature-inducing proteins (Figure 4). The Hill-type ultrasensitivity suggests that the effect of cortical tension on membrane deformation is captured by the change in the propensity of curvature-inducing proteins to self-assemble on the membrane for generating high curvature structures.

Figure 4

Schematic of how ECM stiffness affects the formation of curved structures, essential in endocytic pathways, by curvature-inducing proteins

Here, the generation of highly curved structures (such as vesicles and tubules) as a function of excess area A/A (proxy for cortical tension) has been quantified using the results from our biophysical model (Figure 3C) as the number of protein fields in membrane protrusions n) normalized by the total number of protein fields (n). Fold change in the curvature generation (i.e. ratio of the normalized metric) with increase (or decrease) of cortical tension when going from low ECM stiffness to high ECM stiffness is shown for different amounts of curvature-inducing proteins on the membrane. MVE: multivesicular endosomes, ECM: Extracellular Matrix. Created with Biorender.com.

Our model predicts that a decrease in cortical tension could lead to a higher propensity for curvature generation. However, increasing cortical tension does not necessarily mean a decrease in curvature induced by proteins, as the cortical tension increase could be within the saturation region of hill-type dependence (Figure 4). This dependency suggests that the influence of ECM stiffness on membrane deformation depends on which region (saturated or sensitive) of the Hill-type function and by how much the cortical tension changes, which would be dependent on the cell type. For instance, an increase in cortical tension with stiffening ECM could be within the saturation region and thereby, not significantly affect the generation of highly curved structures.

Recent evidence indicates that tension experienced by membranes also modulates the recruitment dynamics of curvature-inducing proteins with upregulation in binding to membranes at low tension (Kai et al., 2021; Mercier et al., 2020; Tsujita et al., 2021). This tension-dependent recruitment of curvature-inducing proteins combined with our model’s prediction of curvature generation by recruited proteins being sensitive to cortical tension suggests an interesting mechanosensitive regulatory axis that could be exploited to alter/hijack intracellular trafficking mechanisms by cancer cells to gain fitness advantage.

The plasticity of cancer cells along the epithelial-to-mesenchymal transition (EMT) axis has been shown to be regulated by biomechanical properties of the microenvironment such as matrix stiffness. Nuclear localization of EMT transcription factor TWIST1 is promoted by the mechanosensitive EPHA2/LYN protein complex, wherein EPHA2 signaling promoted by high matrix stiffness activates LYN kinase which then phosphorylates cytoplasmic TWIST1 (Wei et al., 2015; Fattet et al., 2020). Using a computational framework, Deng and co-workers demonstrated a positive feedback loop involving EMT transcription factor ZEB1, collagen cross-linking enzyme LOXL2 and ECM stiffness (Deng et al., 2021). EMT-induced morphology has been reported to have lower plasma membrane tension as compared to parental cells, leading to increased assembly of BAR proteins on the membrane for cell motility (Tsujita et al., 2021). Cells can also exploit endocytic machinery for advancing changes (such as cell polarity) to gain mesenchymal phenotype, and/or carry out EMT signaling (Corallino et al., 2015). Based on these studies, coupled with our insights about the effect of ECM stiffness on cell stiffness and vesicle trafficking, it is tempting to speculate about the intricate dual role of ECM stiffness in EMT-based invasive phenotype, involving biophysical and biochemical axes.

While we have focused on how ECM stiffness might regulate membrane curvature deformation through cortical tension, it is essential to consider the signaling pathways involved in cargo trafficking being altered in cancer. The importance of mechano-based mechanisms being employed by tumors in the deregulation of trafficking is emphasized by the recent study by Patwardhan and coworkers (Patwardhan et al., 2021) who showed ECM stiffness-dependent regulation of exosome secretion through the YAP-TAZ pathway in breast cancer for MDA-MB-231 and MCF-7 cell lines with higher exosome release on stiff matrix.

The Rab family of small GTPases is pivotal in several signaling steps, ranging from transport to fusion of membrane-bound organelles and vesicles, of trafficking pathways (Grant and Donaldson, 2009). Numerous studies have reported the dysregulation of Rab-mediated vesicle trafficking associated with tumor progression (Tzeng and Wang, 2016; Lanzetti and Fiore, 2017). A vital signaling component of trafficking involves the sorting of various cargo. Although in normal cellular behavior endosomal sorting is tightly regulated, growing evidence suggests that the altered sorting pathways could be key contributors to cancerous behavior. Case in point, aberrant functioning of prominent protein families involved in sorting, such as sorting nexins (SNXs; orchestrates endosomal sorting (Cullen, 2008)) and Cbl (mediates targeting of receptor tyrosine kinases for lysosomal degradation (Swaminathan and Tsygankov, 2006)), has been noted to have oncogenic roles (Liyasova et al., 2015; Hanley and Cooper, 2021). For tethering of secretory vesicles to the cellular membrane in order to fuse and subsequently release their content into the extracellular milieu, cells employ the evolutionary conserved exocyst complex consisting of 8 (EXOC) subunits (Heider and Munson, 2012). Components of the exocyst complex also have roles in a diverse range of cellular functions (such as cell migration, nanotube formation, EMT, and ciliogenesis), and they have been implicated in the carcinogenesis of different types of tumors (Tanaka et al., 2017; Zago et al., 2019; Zaman et al., 2021; Saha et al., 2022). An important class of trafficking vesicles is autophagosomes, involved in autophagy that plays a critical role in cancer cell survival under conditions of stress such as nutrient limitation (Amaravadi et al., 2016). Autophagy-related (ATG) proteins, regulatory protein family that controls the dynamic membrane events in autophagosome biogenesis (Mizushima et al., 2011), have emerged as promising targets for cancer therapy with several inhibitory drugs in clinical and preclinical trials (Chen et al., 2014; Mulcahy Levy and Thornburn, 2020). These examples underpin the importance of altered signaling components of trafficking pathways in conferring an array of functional advantages to cancer cells.

To illustrate how these key signature proteins (Rab, SNX, EXOC, and ATG) involved in regulating vesicle trafficking might be altered in the patients of cancer types considered in Figure 1, we used the cBioPortal (http://cbioportal.org) (Cerami et al., 2012; Gao et al., 2013) to query for and analyze the mutations and copy number variations of the signature proteins in the extensive patient-derived cancer genomics data from The Cancer Genome Atlas (TCGA). Table 2 summarizes the results from this genomic analysis. We can see that a significant proportion of patients with cancer across different cancer types showcase genotypic variations in the protein families critical in trafficking pathways, further highlighting the importance of the signaling (regulatory) portion of intra- and inter-cellular cargo trafficking machinery in tumorigenesis.

Table 2

Genomic analysis performed using cBioPortal on patient-derived cancer genomic data from TGCA PanCancer Atlas

Tissue Type	Samples Profiled		Rab		SNX		EXOC		ATG
			max	≥4%	max	≥4%	max	≥4%	max	≥4%
Skin	440	Mut%	5	Rab27B	5.7	SNX29, 13, 19, 31, 9	5	EXOC1, 7, 4, 6B	10	ATG2A, 2B, ULK1, 2
	367	CNA%	3.8		2.7		4.6	EXOC3	3
Prostate	494	Mut%	0.6		0.8		0.6		0.8
	489	CNA%	4.3	Rab2A, 3C	8.6	SNX14, 31, 3, 16	2.2		8.6	ATG5
Ovary	523	Mut%	0.6		1		1		0.8
	572	CNA%	10.8	Rab8A, 19, 3D, 6A, 30, 25, 2A, 4A, 38, 3A, 13, 22A, 39B, 44	12.4	SNX31, 27	9.1	EXOC3, 4, 2, 8	7.9	ATG4D, 9B, 16L2, 3
Lung	566	Mut%	1.2		2.5		3		3.7
	511	CNA%	9.4	Rab13, 25, 22A, 29	10	SNX27, 6, 16, 31	12.9	EXOC3	2
Colon	534	Mut%	2.6		4.9	SNX13	3.7		6.2	ATGA2A, 2B, 9B
	592	CNA%	7.6	Rab22A	7.6	SNX21	0.8		1.5
Breast	1066	Mut%	0.7		0.7		1		1
	1070	CNA%	9.2	Rab4A, 29, 13, 25, 2A, 22A, 37, 26, 34	11.1	SNX31, 27, 16	9.4	EXOC8	4	ATG16L2
Brain	397	Mut%	0.5		1.3		1.3		2
	575	CNA%	1.4		1.7		5.2	EXOC1	1.6

Mut% and CNA% correspond to the percentage of samples with one or more mutations and copy number alterations respectively. Samples profiled is the number of patient samples having the genomic data needed for that particular analysis. Here, max is the maximum value of Mut% (CNA%) for a protein family and ≥ 4% refers to the proteins having Mut% (CNA%) greater than 4%. See Table S1 for details about the TCGA datasets used for each of the tissue types. Rab: Rab GTPases, SNX: sorting nexins, EXOC: exocyst complex components, and ATG: autophagy-related proteins.

Future work will be focused on integrating signaling pathways triggered by mechanical cues from ECM with the biophysical model to understand better how vesicular biogenesis and trafficking are (de)regulated in tumors by mechanobiology-based mechanisms. This would extend the mechanosensitive paradigm for cellular functioning presented in this study wherein tumors exploit the connection between cortical tension and biogenesis of tubulo-vesicular membranous compartments along with altered signaling components in trafficking pathways to survive and thrive.

Limitations of the study

Single-cell proteomics data was not available in the machine learning model, and average protein spectral counts were used instead of single-cell data. An additional limitation is that all the data in NCI Leidos Mechanobiology dataset is derived from subconfluent cells on a flat substrate, and therefore additional factors owing to cell-cell contacts and three dimensionality are missing. In our mesoscale membrane model, we have focused on surface deformation owing to curvature-inducing proteins. However, several other mechanisms also contribute to curvature generation. For instance, actin polymerization and motor protein forces can directly contribute to highly curved structure generation processes in endocytic pathways, such as the elongation of tubules in recycling endosomes and contribution to vesicle secession at the plasma membrane (Anitei and Hoflack, 2011). Cell membrane adhesion, which can be seen from a biophysical perspective as pinning of membranes (such as membrane-cytoskeleton linker proteins), has also been implicated in promoting curvature generation (Kutty Kandy and Radhakrishnan, 2019). Another factor to consider is the glycocalyx-induced membrane remodeling (Kutty Kandy and Radhakrishnan, 2022). We have also not undertaken in vivo experimental studies for the signature proteins identified by our analysis. Given the prior supporting data discussed in this study, future work focusing on in vivo experiments (such as live-cell interventions and patient-derived xenograft models) could provide further biological relevance to the observations from our investigation and identify potential druggable targets.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ravi Radhakrishnan.

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Cell growth

All cell lines were from ATCC (Manassas, VA) and grown according to ATCC standard operating procedures (SOPs), which are provided in the Synapse site along with the source data (see key resources table). Media and reagents used for cell culture were either provided by ATCC or purchased from an outside supplier as per ATCC SOP instructions. Thawed cells were taken to be primary cultures (p0). Cells were passaged until passage 3 (p3) and plated on polyacrylamide or HA gels at single cell density (25,000 cells/gel). Experiments were performed 24 hours after cell plating.

Polyacrylamide gel fabrication

Polyacrylamide gels were created by using a combination of acrylamide and bisacrylamide from Bio-Rad Laboratories (Hercules, CA), along with water and NHS (N-hydroxysuccinimide ester) from Sigma (St. Louis, MO) dissolved in toluene (Fisher Scientific, Waltham, MA). Briefly, NHS was dissolved in toluene to create a saturated solution. The solution was then spun for 5 minutes on a tabletop centrifuge to remove any undissolved NHS. The NHS-toluene solution was then added to the water/polyacrylamide/bisacrylamide mixture and vortexed for ten seconds. To make the 30kPa polyacrylamide gels, a final concentration of 13.5% acrylamide was used, while 500 Pa polyacrylamide gels had a final concentration of 3.5%. Both gels had a bisacrylamide concentration of 0.33%. Acrylamide and bisacrylamide were from BioRad Laboratories (Hercules, CA). After vortexing the solution, we centrifuged at 1000 rpm in a tabletop centrifuge to separate the toluene from the polyacrylamide solution. We then removed the polyacrylamide solution from under the toluene layer, being careful not to disturb the toluene. We then aliquoted out the polyacrylamide. To begin polymerization, we then added to a final concentration of 0.3% and 0.06% of Temed and APS (ammonium persulfate solution) respectively. Our solution was immediately pipetted onto functionalized coverslips and a siliconized coverslip was added on top. Gels were allowed to polymerize for 20 minutes.

Top coverslips were removed, and gels were rinsed twice for 5 minutes each in PBS. To crosslink protein to the surface of the gel, the gels were incubated in a 0.1 mg/mL solution of collagen (BD Biosciences, San Jose, CA) or fibronectin (EMD Millipore, Billerica, MA) in HEPES pH 8.0 either overnight at 4C or at RT for 4 hours. Gels were then rinsed and kept in PBS for storage. One additional set of 500 Pa Collagen gels were made along with the other PAA gels but containing fluorescent beads and were sent to our collaborators for volume and traction force measurements.

Hyaluronic acid gel fabrication

HA gels were created by polymerizing thiol-modified hyaluronic acid with Extra Link (PEGDA) (Ascendance Biotechnology, Alameda, CA) according to the manufacturer’s instructions. A 1mL vial of lyophilized HA was reconstituted with 875ul of degassed, deionized water for 30 minutes at 37C. After 30 minutes, 125ul of a 1 mg/mL solution of either collagen or fibronectin was added to the vial. Gels were polymerized at a ratio of 1-part Extra Link to 4 parts HA/protein solution. The desired volume of this mixture was pipetted on a glutaraldehyde-functionalized glass cover slip and then a siliconized coverslip was placed on top. Gels were allowed to polymerize for twenty minutes before removing the top coverslip and storing in PBS.

Method details

NCI Leidos mechanobiology dataset

The National Cancer Institute (NCI) commissioned the Leidos study to understand the relative importance of substrate stiffness on the structure, motility, and proliferation rate of cancer cells in substrates with defined stiffness and chemical composition. 25 different human cancer cell lines derived from human tumors from the breast, colon, brain, ovary, pancreas, prostate, and skin were studied along with five immortalized but non-tumorigenic cells from the same set of tissues. 20 of the cell lines had a complete data universe comprising the measurements of 9 physical quantities (cell stiffness, motility, adherent area, circularity, aspect ratio, proliferation, cell volume, nuclear volume, and contractility) characterizing the cell phenotype, of which 9 cell lines (A375, RWPE-1, 22RV1, SW620, SW480, T-47D, MDA-MB-231, U-87, and T98G) had measurements of genetic variants, gene and transcript expression values, and proteomics data. These cells were grown in 7 different ECM substrates (of varying stiffness); namely, Glass, 500Pa Fibronectin (FN)-coated polyacrylamide gels, 30kPa FN-coated polyacrylamide gels, 500Pa Collagen (COL) I-coated polyacrylamide gels, 30kPa Collagen I-coated polyacrylamide gels, 500Pa FN-coated hyaluronan (HAFN) gels, 500Pa Collagen I-coated hyaluronan (HACOL) gels. The dataset was created to provide a corpus of well-characterized molecular and cellular phenotypes measured on the same cell lines in the same substrate conditions that would provide a rich source of information about making inferences between these phenotypes and their corresponding genotypes. We focus on the analysis of the cell stiffness dataset as characterized by the Young’s Modulus (E) measuring cell stiffness (Figure 1A) and the proteomics data specifying the intensity-based absolute quantitation (iBAQ) spectral counts of the protein expression of the cell types in various substrates.

AFM, RNA-seq, proteomics

Detailed description of the experimental procedure followed for obtaining Young’s modulus via atomic force microscopy (AFM), RNA-seq data in Leidos genomic dataset, and proteomics data used in this study are provided on Synapse along with the dataset. DOIs are listed in the key resources table.

KDEs of cell stiffness data

We used SciPy 1.6.2 (Virtanen et al., 2020), an open-source Python library for scientific and technical computing, to compute the kernel density estimates (KDE, estimation of the probability density function) of the Young’s modulus E using Gaussian kernels with bandwidth estimator as Scott’s Rule. KDEs were calculated: across all cell types and substrates, and across all substrates for each cell type (Figure S1). In order to obtain a more discrete view of the data, we introduced a binarization scheme using the median of the kernel density estimate (either across all cell types and all substrates or across all substrates for each cell type) to split the Young’s modulus data into low or high categories based on if the value is below or above the median respectively. To quantify the distribution of cell count in the low E and high E regime, we computed the ratio of the count of cells in high E over low E.

Mutual information calculations

The coupling of the phenotype (cell stiffness in our analysis) to proteomics was achieved through the use of an information-theoretic measure, mutual information (I). Mutual Information (I) is a statistical measure of two random variables X, Y which describes the amount of information contained in one random variable relative to the other. In this study, we used normalized mutual information such that 0≤I≤1 . Therefore, a normalized mutual information I with a value close to 1 indicates a high informational dependency between the two variables X, Y. We used this framework to select proteins that have high mutual information with the cell phenotype. We assume that the protein spectral count and the cell stiffness (E) are random variables generated by some biologically driven stochastic process. We are interested in finding the amount of information that the protein spectral count contains about the cell phenotype variable E. The MIT licensed library of information-theoretic methods pyitlib 0.2.2 was used to calculate the normalized mutual information.

Neural network, PCA and shapley analysis

The proteins that are selected through the mutual information-driven process are used as features in a neural network to predict if a cell belongs to low E or high E regime (where low or high is determined if below the median or above the median respectively). Note that the machine learning model for proteomics is applied to bulk level data while the phenotype data is at single cell resolution. The neural network was implemented in TensorFlow 2.5.0 (Abadi et al., 2016) using the Keras API with one single hidden layer. For training, ‘binary_crossentropy’ was used as the loss function, and Adam optimizer was used with 0.001 as the learning rate. A 70:30 split of the dataset was done to obtain a training set and test set. K-fold cross-validation (with k = 5) was used to determine the optimal number of neurons (hyperparameter) in the hidden layer, which was found to be 5 neurons. Scikit-learn’s (Pedregosa et al., 2011) balanced accuracy metric was used to estimate the accuracy of the neural network prediction.

Principal component analysis (PCA) of the protein spectral count features was performed using Scikit-learn’s PCA library. The Shapley game-theoretic metric to calculate the feature importance of the neural network protein spectral count features was implemented using the SHAP python package (Lundberg and Lee, 2017).

Continuum membrane model

The core methodology followed for performing the simulations is the same as that reported in Tourdot et al. (2015). In brief, the membrane patch is modeled as a thin elastic sheet and discretized into vertices, each of characteristic size a0, interlinked by L links that form T triangles. The thermodynamic behavior of the membrane is well captured by the elastic energy functional given by:which is a discretized form of the Canham-Helfrich Hamiltonian (Helfrich 1973); where the material properties and σ represent the bending rigidity and the bare surface tension, respectively. C1, and C2, are the curvatures at vertex v along the two principal axes, computed as in Ramakrishnan et al. (2010), and A is the corresponding surface area. H0, is a spontaneous curvature field at vertex v which captures the curvature-inducing interactions between the protein and membrane and is given bywhere r and R denote the position of vertex v and i protein field, respectively. This approach of treating the effect of the curvature-inducing proteins as a curvature field in the continuum field formulation has been utilized in prior studies (Agrawal et al., 2010; Liu et al., 2012; Zhao et al., 2013). See Ramakrishnan et al., 2014 for a comprehensive review on the methodological aspects of the membrane model used in this work.

The results presented in this study are for a membrane patch with N = 900, and σ = 0. A typical value of for the systems we model is a0 = 10nm; this value ensures that the and used in this study are primed to model curvature-inducing proteins such as BAR and ENTH domains, as justified in previous works (Tourdot et al., 2014, 2015). Note that a detailed sensitivity analysis for the mesoscale model has been presented in an earlier study (Tourdot et al., 2015) for the parameters of the protein fields (C0 and ). Monte Carlo techniques evolve the membrane morphology and diffusion of protein fields on the membrane patch. The membrane is allowed to equilibrate for 3 million MC steps. The number of protein fields in the tubular regions (n) are calculated using a clustering algorithm with a mean curvature (H) cutoff of 0.5 , i.e. membrane regions with tubulation are characterized as (see Tourdot et al., 2015 for detailed discussion). All reported data is averaged over four independent ensembles.

The non-linear fits for Equation 1 to the simulation data at different values of n (total number of curvature-inducing protein fields on the membrane) were obtained using the curve_fit function of python library SciPy 1.6.2. Parameters used for curve_fit were ‘lm’ (Levenberg-Marquardt algorithm) as the method for optimization, [1, 0, 0] as the initial condition for [‘a’, ‘b’, ‘c’] in Equation 1, and 8000 as maximum function evaluations.

Excess area dependence from micelle model

In our previous work (Tourdot et al., 2015), we showed the analogy between membrane tubulation and micellization, and determined the relation between the number of protein fields in the planar and tubular regions aswhere n is total number of proteins fields and is the number of protein fields in the planar region. N is the absolute number of protein units within each tubule given as N = nζ, where n is the average number of protein fields per tubule (which from simulation was found to be ) and ζ refers to the number of protein units per curvature-inducing protein field. Critical aggregation number (n refers to the critical number of protein fields that promote membrane tubulation (see Tourdot et al. (2015) for detailed discussion). Taking n to be number of tubes each containing n protein fields, the number of protein fields in the tubular region (n is then given by n = nn. Figure S3 shows the excellent agreement between the micelle model in Equation 4 and results from our simulation for different excess areas. Substituting in Equation 4, we get

In asymptotic limit of , the above expression can be simplified to

The value for n for different excess area (A/Ap was obtained from simulation data, and a non-linear fit for n as a function of A/A was evaluated (see Figure S3D) which is given by

Substituting the above expression for n in Equation 8, we get

Equation 10 suggests that the relation between n and A/A can be written as

cBioPortal: genomic analysis of patient data

The cancer genomic data exploration platform of cBioPortal (Cerami et al., 2012; Gao et al., 2013) was used to obtain measures such as percentage of mutations, percentage of copy number variations for protein families of Rab GTPases, sorting nexins (SNX), exocyst complex component (EXOC) and autophagy-related proteins (ATG) for different tumor types. For copy number variations, both amplifications and deep deletions were considered. Table S1 lists the tissue type and the corresponding TGCA PanCancer Atlas data used for the analysis.

cBioPortal: survival analysis

Survival analysis was performed using the cBioPortal (Cerami et al., 2012; Gao et al., 2013) web server for all the tumor types listed in Table S2. For each of the selected cytoskeleton genes, the control and subject groups were defined based on the normalized mRNA expression ( -score) of the gene. Here, -score is defined aswhere r is the raw expression value, μ and σ are the mean and standard deviation of the reference population (taken as all profiled tumor samples of a particular cancer type), respectively. Note that we used log transformed data for computing the -score. For each gene, two subject groups were defined, (i) samples having low expression ( -)), and (ii) samples having high expression ( -score ≥ threshold). Control group was defined as samples satisfying -. The Onco Query Language (OQL) was used to define the control and subject groups on the cBioPortal. See the OQL documentation on cBioPortal for more details. The queries used in our analysis were ‘: EXP <= -threshold’ and ‘: EXP >= threshold’ for defining the subject groups. The samples not satisfying these conditions are separately grouped as unaltered (control group) by cBioPortal. Subsequently, the Kaplan-Meier plots for overall and progression free survival, produced on cBioPortal using the clinical data of the patients in control and subject groups, were analyzed. We examined values of 1.0 and 1.5 (higher values had too few events in subject groups and thereby, made their KM curves much less reliable).

Quantification and statistical analysis

Log-rank test was used to compare the survival curves of two groups on the cBioPortal web server, and p-value < 0.05 was considered a significant difference.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
NCI Leidos mechanobiology dataset (PS-ON Cell Line Characterization)	Synapse	Synapse data: https://doi.org/10.7303/syn7248578
Code for analysis of the mechanobiology dataset	This paper	https://github.com/jnukpezah/iScience_StiffnessCode
Experimental models: Cell lines
Human: SK-MEL-2	ATCC	HTB-68;RRID:CVCL_0069
Human: A375	ATCC	CRL-1619;RRID:CVCL_0132
Human: WM266-4	ATCC	CRL-1676;RRID:CVCL_2765
Human: RWPE-1	ATCC	CRL-11609;RRID:CVCL_3791
Human: 22Rv1	ATCC	CRL-2505;RRID:CVCL_1045
Human: Panc-1	ATCC	CRL-1469;RRID:CVCL_0480
Human: SK-OV-3	ATCC	HTB-77;RRID:CVCL_0532
Human: NCI-H2126	ATCC	CCL-256;RRID:CVCL_1532
Human: NL20	ATCC	CRL-2503;RRID:CVCL_3756
Human: HCT116	ATCC	CCL-247;RRID:CVCL_0291
Human: HT-29	ATCC	HTB-38;RRID:CVCL_0320
Human: SW620	ATCC	CCL-227;RRID:CVCL_0547
Human: SW480	ATCC	CCL-228;RRID:CVCL_0546
Human: T-47D	ATCC	HTB-133;RRID:CVCL_0553
Human: MCF7	ATCC	HTB-22;RRID:CVCL_0031
Human: MDA-MB-231	ATCC	CRM-HTB-26RRID:CVCL_0062
Human: hTERT-HME1	ATCC	CRL-4010;RRID:CVCL_3383
Human: MCF10A-JSB	ATCC	CRL-10317;RRID:CVCL_VH36
Human: U-87	ATCC	HTB-14;RRID:CVCL_0022
Human: T98G	ATCC	CRL-1690;RRID:CVCL_0556
Software and algorithms
Python 3.9.12	Python Software Foundation	https://www.python.org/
SciPy v1.6.2	Virtanen et al., 2020	https://scipy.org/
pyitlib v0.2.2	MIT licensed library	https://pafoster.github.io/pyitlib/
TensorFlow v2.5.0	Abadi et al., 2016	https://www.tensorflow.org/
scikit-learn v1.0.2	Pedregosa et al., 2011	https://scikit-learn.org/stable/
SHAP python package	Lundberg and Lee, 2017	https://shap.readthedocs.io/en/latest/index.html
Source code for DTMC membrane model	Ramakrishnan et al., 2014	https://doi.org/10.1016/j.physrep.2014.05.001
Other
ATCC SOPs for cell growth	Synapse	https://www.synapse.org/#!Synapse:syn7248578/wiki/410632
Experimental procedure for obtaining Young’s modulus using atomic force microscopy	Synapse	https://doi.org/10.7303/syn7248585
Experimental procedure for RNA-seq data in Leidos genomic dataset	Synapse	https://doi.org/10.7303/syn11510898
Experimental procedure for obtaining proteomics data in Leidos dataset	Synapse	https://doi.org/10.7303/syn9697791
cBioPortal : web application for analysis of TCGA database	Cerami et al., 2012; Gao et al., 2013	http://cbioportal.org