Population Heterogeneity in the Epithelial to Mesenchymal Transition Is Controlled by NFAT and Phosphorylated Sp1

Abstract

Epithelial to mesenchymal transition (EMT) is an essential differentiation program during tissue morphogenesis and remodeling. EMT is induced by soluble transforming growth factor [Beta] (TGF-[Beta]) family members, and restricted by vascular endothelial growth factor family members. While many downstream molecular regulators of EMT have been identified, these have been largely evaluated individually without considering potential crosstalk. In this study, we created an ensemble of dynamic mathematical models describing TGF-[Beta] induced EMT to better understand the operational hierarchy of this complex molecular program. We used ordinary differential equations (ODEs) to describe the transcriptional and post-translational regulatory events driving EMT. Model parameters were estimated from multiple data sets using multiobjective optimization, in combination with cross-validation. TGF-[Beta] exposure drove the model population toward a mesenchymal phenotype, while an epithelial phenotype was enhanced following vascular endothelial growth factor A (VEGF-A) exposure. Simulations predicted that the transcription factors phosphorylated SP1 and NFAT were master regulators promoting or inhibiting EMT, respectively. Surprisingly, simulations also predicted that a cellular population could exhibit phenotypic heterogeneity (characterized by a significant fraction of the population with both high epithelial and mesenchymal marker expression) if treated simultaneously with TGF-[Beta] and VEGF-A. We tested this prediction experimentally in both MCF10A and DLD1 cells and found that upwards of 45% of the cellular population acquired this hybrid state in the presence of both TGF-[Beta] and VEGF-A. We experimentally validated the predicted NFAT/Sp1 signaling axis for each phenotype response. Lastly, we found that cells in the hybrid state had significantly different functional behavior when compared to VEGF-A or TGF-[Beta] treatment alone. Together, these results establish a predictive mechanistic model of EMT susceptibility, and potentially reveal a novel signaling axis which regulates carcinoma progression through an EMT versus tubulogenesis response.

About the Authors:

Russell Gould

Contributed equally to this work with: Russell Gould, David M. Bassen

Affiliation: Nancy E. and Peter C. Meinig School of Biomedical Engineering, Cornell University, Ithaca, New York, United States of America

David M. Bassen

Contributed equally to this work with: Russell Gould, David M. Bassen

Affiliation: Nancy E. and Peter C. Meinig School of Biomedical Engineering, Cornell University, Ithaca, New York, United States of America

Anirikh Chakrabarti

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York, United States of America

Jeffrey D. Varner

Affiliation: Robert Frederick Smith School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York, United States of America

Jonathan Butcher

* E-mail: [email protected]

Affiliation: Nancy E. and Peter C. Meinig School of Biomedical Engineering, Cornell University, Ithaca, New York, United States of AmericaAbstract

Epithelial to mesenchymal transition (EMT) is an essential differentiation program during tissue morphogenesis and remodeling. EMT is induced by soluble transforming growth factor [Beta] (TGF-[Beta]) family members, and restricted by vascular endothelial growth factor family members. While many downstream molecular regulators of EMT have been identified, these have been largely evaluated individually without considering potential crosstalk. In this study, we created an ensemble of dynamic mathematical models describing TGF-[Beta] induced EMT to better understand the operational hierarchy of this complex molecular program. We used ordinary differential equations (ODEs) to describe the transcriptional and post-translational regulatory events driving EMT. Model parameters were estimated from multiple data sets using multiobjective optimization, in combination with cross-validation. TGF-[Beta] exposure drove the model population toward a mesenchymal phenotype, while an epithelial phenotype was enhanced following vascular endothelial growth factor A (VEGF-A) exposure. Simulations predicted that the transcription factors phosphorylated SP1 and NFAT were master regulators promoting or inhibiting EMT, respectively. Surprisingly, simulations also predicted that a cellular population could exhibit phenotypic heterogeneity (characterized by a significant fraction of the population with both high epithelial and mesenchymal marker expression) if treated simultaneously with TGF-[Beta] and VEGF-A. We tested this prediction experimentally in both MCF10A and DLD1 cells and found that upwards of 45% of the cellular population acquired this hybrid state in the presence of both TGF-[Beta] and VEGF-A. We experimentally validated the predicted NFAT/Sp1 signaling axis for each phenotype response. Lastly, we found that cells in the hybrid state had significantly different functional behavior when compared to VEGF-A or TGF-[Beta] treatment alone. Together, these results establish a predictive mechanistic model of EMT susceptibility, and potentially reveal a novel signaling axis which regulates carcinoma progression through an EMT versus tubulogenesis response.

Author Summary

Tissue formation and remodeling requires a complex and dynamic balance of interactions between epithelial cells, which reside on the surface, and mesenchymal cells that reside in the tissue interior. During embryonic development, wound healing, and cancer, epithelial cells transform into a mesenchymal cell to form new types of tissues. It is important to understand this process so that it can be controlled to generate beneficial effects and limit pathological differentiation. Much research over the past 20 years has identified many different molecular species that are relevant, but these have mainly been studied one at a time. In this study, we developed and implemented a novel computational strategy to interrogate the key players in this transformation process to identify which are the major bottlenecks. We determined that NFATc1 and pSP1 are essential for promoting epithelial or mesenchymal differentiation, respectively. We then predicted the existence of a partially transformed cell that exhibits both epithelial and mesenchymal characteristics. We found this partial cell type develops a network of invasive but stunted vascular structures that may be a unique cell target for understanding cancer progression and angiogenesis.

Citation: Gould R, Bassen DM, Chakrabarti A, Varner JD, Butcher J (2016) Population Heterogeneity in the Epithelial to Mesenchymal Transition Is Controlled by NFAT and Phosphorylated Sp1. PLoS Comput Biol 12(12): e1005251. doi:10.1371/journal.pcbi.1005251

Editor: Feilim Mac Gabhann, Johns Hopkins University, UNITED STATES

Received: January 14, 2015; Accepted: November 17, 2016; Published: December 27, 2016

Copyright: © 2016 Gould et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting Information file S1_Text.pdf. The code and files for recreating the model are uploaded to GitHub (https://github.com/varnerlab/Reduced_EMT_Model). If more information is required, authors may be contacted at [email protected] or [email protected].

Funding: The project described was supported by Award Number #U54CA143876 from the National Cancer Institute to JDV, National Science Foundation CBET-0955172 to JB, National Science Foundation Graduate research fellowship to DMB, and National Institute of Health HL110328 to JB. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Cancer Institute or the National Institutes of Health. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The epithelial to mesenchymal transition (EMT) is a broadly participating, evolutionarily conserved differentiation program essential for tissue morphogenesis, remodeling and pathological processes such as cancer [1]. During EMT polarized, tightly adhered epithelial cell monolayers are transformed into non-interacting motile mesenchymal cells that simultaneously degrade and synthesize extracellular matrix (ECM) components and invade into the underlying tissue space [2]. EMT is the fundamental initiator of developmental processes such as embryonic gastrulation and valvulogenesis [3]. Transforming growth factor [Beta] (TGF-[Beta]) family members are important inducers of both developmental and pathological EMT [4, 5]. Decades of research has focused on identifying molecular regulators of EMT, but almost all on a single gene and in a nearly binary yes/no level of qualitative understanding. Medici and coworkers identified a core signaling program by which TGF-[Beta] isoforms induce EMT across a variety of cell lines [6, 7]. This program involves carefully orchestrated rounds of gene expression driven by the Smad and Snail families of transcription factors as well as other key factors such as lymphoid enhancer-binding factor 1 (LEF-1), nuclear factor of activated T-cells, cytoplasmic 1 (NFATc1), and specificity protein 1 (Sp1). Coregulators such as [Beta]-catenin, NF-KB, and the ErbB family of receptor tyrosine kinases however also participate in EMT regulation, but the degree of each's influence is difficult to ascertain in isolation [8-11]. EMT also exhibits complex temporal dynamics that are often intractable in gain/loss of function studies. Elucidating the master regulatory architecture controlling EMT therefore requires inclusion of these complex overlapping and non-binary behaviors.

Systems biology and mathematical modeling are essential tools for understanding complex developmental programs like EMT [12]. Previous computational models of TGF-[Beta] induced differentiation focused on single biological factors or EMT in single cells. For example, Chung et al., constructed a model of TGF-[Beta] receptor activation and Smad signaling using ordinary differential equations and mass-action kinetics. Their model suggested that a reduction of functional TGF-[Beta] receptors in cancer cells may lead to an attenuated Smad2 signal [13]. Similarly, Vilar et al. suggested that specific changes in receptor trafficking patterns could lead to phenotypes that favor tumorigenesis [14]. Coarse grained modeling approaches have also been applied to EMT; Steinway et al. used discrete dynamic modeling to study developmental EMT and known dysregulation in invasive hepatocellular carcinoma [15]. Although these models provided insight into the role of receptor dynamics, EMT induction involves many other components, including competing second messengers and interconnected transcriptional regulatory loops. Integrating these additional scales of molecular signaling while maintaining the capacity for robust prediction requires a new and expanded computational and experimental strategy. Data-driven systems approaches [16] or logical model formulations [17] are emerging paradigms that constrain model complexity through the incorporation of training and validation data. These are interesting techniques because the data informs model structure (which can be expanded as more data becomes available). Alternatively, Bailey proposed more than a decade ago that a qualitative understanding of a complex biological system should not require complete definition of its structural and parametric content [18]. Shortly thereafter, Sethna and coworkers showed that complex model behavior is often controlled by only a few parameter combinations, a characteristic seemingly universal to multi-parameter models referred to as "sloppiness" [19]. Thus, reasonable model predictions are often possible with only limited parameter information. Taking advantage of this property, we developed sloppy techniques for parameter identification using ensembles of deterministic models [20]. Furthermore, we proposed that the sloppy behavior of biological networks may also be seen as a source of cell-to-cell [21] or even patient-to-patient heterogeneity [22]. Bayesian parameter identification techniques have also been used to explore cell-to-cell heterogeneity [23, 24], where a population of cells could be viewed as a dynamic ensemble of context-specific biochemical networks [25].

In this study, we developed a family of mathematical models describing the induction of EMT by TGF-[Beta] isoforms in the presence and absence of vascular endothelial growth factor A (VEGF-A). We integrated a simple rule-based description of activity and gene expression regulation with traditional ordinary differential equation (ODE) modeling to describe an EMT interaction network containing 97 gene, protein or mRNA components interconnected through 169 interactions. This integration allows the description of complex regulatory interactions in the absence of specific mechanistic information, it also allowed to build a predictive yet compact model. A family of model parameters was estimated using 41 molecular data sets generated in DLD1 colon carcinoma, MDCKII and A375 melanoma cells using the Pareto optimal ensemble technique (JuPOETs) multiobjective optimization algorithm. JuPOETs generated an ensemble of approximately 1400 models for analysis. Analysis of the model population suggested that both MCF10A and DLD1 cells could exhibit phenotypic heterogeneity if treated simultaneously with TGF-[Beta]1/2 and VEGF-A. This heterogeneity was characterized by a significant fraction of the population being in a "hybrid state" having both high E-cadherin and high Vimentin expression. We tested these predictions using qRT-PCR and flow-cytometry studies in a variety of experimental conditions. Validation studies confirmed that upwards of 45% of the cellular population could be put into the hybrid state in the presence of both TGF-[Beta]1/2 and VEGF-A. Moreover, this response depended upon both activation of Sp1 by MAPK and NFATc1 transcriptional activity consistent with the predicted molecular signaling. Lastly, the hybrid populations of both DLD1 and MCF10A cells exhibited different functional behavior than those from either TGF-[Beta] or VEGF-A treatment. The extent of ductal branch formation significantly increased with MCF10A cells in the hybrid phenotype, compared with cells treated with VEGF-A alone. Together, these results establish a predictive mechanistic model of EMT susceptibility, and reveal a novel signaling axis, which possibly regulates carcinoma progression through an EMT versus tubulogenesis response.

Results

The model population captured key features of TGF-[Beta] induced EMT

The EMT model architecture, based upon curated molecular connectivity, described the expression of 23 genes following exposure to TGF-[Beta] isoforms and VEGF-A (Fig 1). The EMT model contained 74 molecular species interconnected by 169 interactions. Model equations were formulated as ordinary differential equations (ODEs) augmented with rule-based descriptions of activity and gene expression regulation. ODEs are common tools to model biochemical pathways [26-28]. However, while ODE models can simulate complex intracellular behavior, they require estimates for model parameters which are often difficult to obtain. The EMT model had 251 unknown model parameters, 169 kinetic constants 38 control constants and 44 saturation constants. In addition to constants appearing in the signaling kinetics, characteristic gene expression and translation parameters were estimated from literature, and then corrected for each protein and transcript (four correction factors per gene product, see materials and methods). As expected, the unknown parameters were not uniquely identifiable given the training data [29]. Thus, instead of identifying a single best fit (but uncertain) model, we estimated a sloppy population of models (each consistent with the training data) by simultaneously minimizing the difference between model simulations and 41 molecular data sets using the Pareto Optimal Ensemble Technique (JuPOETs). The training data were generated in DLD1 colon carcinoma, MDCKII, and A375 melanoma cells following exposure to TGF-[Beta] isoforms [7]. We organized these data sets into 11 objective functions which were simultaneously minimized by JuPOETs. Additionally, we used data generated in this study (S1 Fig), and 12 molecular data sets generated in HK-2 cells following VEGF-A exposure to train VEGF-A responsive model processes [30]. To guard against overfitting, we augmented the multiobjective optimization with leave-one-out cross validation to independently estimate both the training and prediction error for each objective. Thus, we generated 11 different model ensembles. Lastly, we compared model predictions with independent data sets not used during training (both at the molecular and model population levels) to evaluate the predictive power of the parameter ensemble.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 1. Model connectivity recreates the core architecture during EMT.

The EMT network contains 97 nodes (proteins, mRNA, and genes) interconnected by 169 interactions. Central to EMT induction, activation of the MAPK cascade occurs through TGF-[Beta]1/2 binding which activates the AP-1/Sp1 transcriptional axis. AP-1/Sp1 drives an autocrine response of TGF-[Beta]3, which activates the Smad cascade, leading to phenotypic change. Conversely, VEGF-A binding can stabilize an epithelial phenotype through NFAT activation. Downstream activation of [Beta]-catenin signaling due to E-cadherin loss and GSK3 inactivation of [Beta]-catenin confinement is critical to the complete activation of the EMT program. The complete list of molecular interactions that comprise the model is given in the supplement.

http://dx.doi.org/10.1371/journal.pcbi.1005251.g001

JuPOETs generated a population of probable signaling models which captured the multiple phases of EMT induction (Fig 2). JuPOETs sampled well over 104 probable models during each stage of the cross-validation using global random sampling. From this analysis, N ≈ 1400 models were selected for further analysis. The selected models all had the same possible molecular connectivity, but different values for model parameters. Transcription and translation rates, as well as mRNA and protein degradation terms, were set using physical values from the literature [31], and allowed to vary by a scaling factor, see methods. Model selection was based upon Pareto rank, the prediction and training error across all objectives. The model population recapitulated key signaling events following TGF-[Beta] exposure. We subdivided the response to TGF-[Beta] exposure into two phases. First, TGF-[Beta]1/2 signaling initiated a program which downregulated E-cadeherin expression in a MAPK dependent manner while simultaneously upregulating TGF-[Beta]3 expression. Second, TGF-[Beta]3 secretion initiated an autocrine feedback which upregulated the expression of mesenchymal markers such as Vimentin and key upstream transcription factors such as LEF-1 in a SMAD dependent manner. TGF-[Beta]3 expression was also able to sustain [Beta]-catenin release by inhibiting its sequestration by the APC complex through PI3K mediated GSK3, which was captured by the model (S2 Fig). Each phase involved the hierarchal expression and/or post-translational modification of several key transcription factors. During the first phase, stimulation with TGF-[Beta]1/2 (10 a.u.) activated both the SMAD and MAPK pathways. MAPK activation resulted in the phosphorylation of the transcription factor activator protein 1 (AP-1), which in-turn upregulated the expression of Snail, a well established transcriptional repressor (Fig 2A). Snail expression was MAPK-dependent; the MEK inhibitor U0126 blocked AP-1 activation and Snail expression following TGF-[Beta]1/2 exposure (Fig 2A, Lane 3). Overexpression of either Snail or Slug upregulated TGF-[Beta]3 expression (Fig 2C) while simultaneously downregulating E-cadeherin expression (Fig 2F). During the second phase, TGF-[Beta]3 secretion and the subsequent autocrine signaling resulted in the upregulation of mesenchymal marker expression. The TGF-[Beta]3 induced gene expression program involves a complex hierarchy of transcriptional and post-translational regulatory events. Absence of E-cadherin indirectly promoted TGF-[Beta]3 expression through the [Beta]-catenin/TCF4 complex following Snail or Slug expression (Fig 2C, Lane 2 or 3). Conversely, over-expression of E-cadherin inhibited the TGF-[Beta]3 autocrine production by sequestering cytosolic [Beta]-catenin, thereby blocking EMT (Fig 2C, Lane 4 or 5). TGF-[Beta]3 signaled through the Smad pathway to regulate LEF-1 expression and downstream target EMT genes (Fig 2G). TGF-[Beta]3 (10 a.u.) in combination with downstream inhibitors (DN-Smad4 and DN-LEF-1) completely inhibited Vimentin expression, while elevating E-cadherin expression (Fig 2H and 2I).

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 2. Training and validation simulations. The population of EMT models qualitatively captured TGF-[Beta]-induced EMT signaling.

(A-I) The population was generated using JuPOETs and trained using 11 different objective functions (41 data sets) taken from Medici et al. [7]. The model captured the simulated experiments for all cases when compared to randomized controls. (J-L) The model populations were also compared against untrained temporal data to measure the effectiveness as a pure prediction. The western blot data was reproduced from Medici et al. [7]. The intensity of each band was estimated using the ImageJ program. These blot intensities were then used as the training data for the EMT parameter estimation studies.

http://dx.doi.org/10.1371/journal.pcbi.1005251.g002

The predictive power of the ensemble was tested using cross validation and by comparing simulations with data not used for model training. In whole, all of our training objectives were statistically significant (at a 95% confidence interval) compared to a randomized parameter family (N = 100) generated from a random starting point. Conversely, we predicted all of the training objectives, at a 95% confidence interval compared to randomized parameters (Wicoxon non-parametric test). The model also captured the temporal gene expression responses of E-cadherin, pSmad2, and LEF-1 (not used for model training) to within one-standard deviation (up to the 48 hr time-point) (Fig 2J-2L). Taken together, the model captured the key signaling events revealed by Medici et al. [7] that drive the phenotypic conversion. A listing of objective function values resulting from training, cross validation and the random parameter control is given in the supplement (S3 Fig).

TGF-[Beta]1/2 and VEGF-A exposure promotes phenotype heterogeneity through NFATc and phosphorylated Sp1

While we captured the central tendency of many of the molecular features of EMT induction following TGF-[Beta]1/2 exposure, an often neglected but important emergent feature of developmental and pathological programs is population heterogeneity [32]. We (and others) have previously hypothesized that deterministic model ensembles can simulate population behavior, at least at a course grained level [21]. We tested this hypothesis by analyzing the response of the population of EMT models to extracellular cues and then comparing this response to flow cytometry studies. We quantified the phenotypic response of the individual members of the ensemble to TGF-[Beta]1/2 stimulation for two downstream phenotypic markers, Vimentin (mesenchymal) and E-cadherin (epithelial) following the addition of TGF-[Beta]1/2 alone (Fig 3), and/or VEGF-A in combination with NFATc inhibitors (Fig 3).

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 3. Simulated VEGF-A and TGF-[Beta]1/2 exposure promoted phenotype heterogeneity.

Simulated response to TGF-[Beta]1/2 and VEGF-A exposure with and without axis specific inhibitors. Vimentin and E-cadherin abundances (in nM) were used to quantify the shift in population at 48 hrs. (A-C) VEGF-A (50 a.u.) treatment resulted in a population with enhanced epithelial properties. This was contrary to the addition of TGF-[Beta]2 (10 a.u.), which shifted the population towards a mesenchymal phenotype. Interestingly, the combined effects of TGF-[Beta]2 and VEGFA was found to increase both ecadherin and vimentin levels, creating a heterogeneous population (black arrow), which can also be seen in a minority of untreated cells (gray arrow). (D-F) To isolate the effect of NFAT, we inhibited NFAT de-phosphorylation in combination with VEGFA. This negated the increase in ecadherin expression and shifted the population towards a mesenchymal phenotype. Likewise, combining NFAT inhibition with TGF-[Beta] mitigated all VEGF enhanced ecadherin expression. Lastly, combination of TGF-[Beta]2, VEGFA, and NFAT inhibition nearly mitigated all effects of VEGFA, shifting the heterogeneous population towards a mesenchymal phenotype. In whole, high levels of phosphorylated-Sp1 correlated with vimentin expression, while NFAT was responsible for maintaining E-cadherin expression in the presence of other factors, although neither were mutually exclusive.

http://dx.doi.org/10.1371/journal.pcbi.1005251.g003

We identified model subpopulations that exhibited different behaviors following exposure to TGF-[Beta]1/2 (Fig 3B). Analysis of the molecular signatures of these subpopulations suggested the abundance, localization and state of the Sp1, AP-1 and NFATc transcription factors controlled population heterogeneity. The majority of models (>80%) responded to treatment, moving away from the untreated population (Fig 3A-3F, gray). These models showed the classically expected behavior, a switch from an epithelial to mesenchymal phenotype following TGF-[Beta]1/2 exposure. Some models resembled untreated cells; they had elevated phosphorylated Sp1, relative to non-induced cells, which decreased E-cadherin expression through Slug-mediated inhibition, which in turn increased Vimentin expression through TGF-[Beta]3 autocrine signaling and the liberation of [Beta]-catenin. However, the most biologically interesting behavior was exhibited by cells achieving a hybrid phenotype, most notable in a dual treatment condition (Fig 3C, black arrow), but also present in a small percentage of untreated cells (Fig 3B, gray arrow). Models with this hybrid phenotype had elevated Sp1 and NFAT transcriptional activity, resulting in simultaneously increased Vimentin and E-cadherin expression (Fig 4A). However, these conclusions are likely sensitive to the components that we have included in the model, and could also involve proteins that we have not considered.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 4. Analysis of underlying signaling responses.

(A) We examined the distribution of NFATc1 and AP1/SP1 in cells containing the hybrid phenotype (VEGF-A + TGF-[Beta]2 case), showing the potential for cells to express both SP1 and NFATc1 in a non exclusive manner. (B) In the absence of YREG1, most of the population failed to consistently to retain a stable epithelial state. (C) We identified a novel regulator of LEF1 called YREG1 that allows Snail/Slug to emulate an inducer by repressing YREG1, which was required to stabilize the untreated population. YREG1 overexpression revealed an enhanced epithelial phenotype, while some inherently transformed cells moved towards a hybrid phenotype.

http://dx.doi.org/10.1371/journal.pcbi.1005251.g004

To better understand the hybrid phenotype, we simulated the response of the model population to TGF-[Beta]1/2 and VEGF-A treatment with and without NFATc inhibitors (Fig 3). As expected, stimulation with VEGF-A (50 a.u.) maintained an epithelial population (Fig 3A), while TGF-[Beta]1/2 (10 a.u.) exposure shifted the population from an epithelial to a mesenchymal phenotype (Fig 3B). On the other hand, combined stimulation with TGF-[Beta]1/2 (10 a.u.) and VEGF-A (50 a.u.) increased both E-cadherin and Vimentin expression, resulting in a hybrid phenotype with both epithelial and mesenchymal characteristics (Fig 3C). Vimentin expression was correlated with high levels of nuclear phosphorylated Sp1, following TGF-[Beta]1/2 exposure. Conversely, elevated E-cadherin expression depended upon the activity of NFAT transcription factors downstream of VEGF-A stimulation. To further isolate the role of NFAT on this hybrid state, we simulated the inhibition of NFAT transcriptional activity across all conditions (all else being equal). NFAT inhibition in combination with VEGF-A or TGF-[Beta]1/2 treatments blocked increased E-cadherin expression in the case of VEGF-A (Fig 3D), but did not influence TGF-[Beta]1/2 signaling (Fig 3E). Lastly, NFATc inhibition in combination with simultaneous TGF-[Beta]1/2 and VEGF-A exposure repressed nearly all E-cadherin expression, shifting nearly the entire population towards a mesenchymal phenotype (Fig 3F). Taken together, high levels of nuclear localized phosphorylated Sp1 correlated with Vimentin expression, while NFATc transcriptional activity was critical for maintaining E-cadherin expression in the presence of competing signals.

Identification of a novel LEF-1 regulator

During model identification, we found that consistent TGF-[Beta] induced EMT from a stable epithelial cell population required an additional regulatory protein. This protein, which we called hypothetical regulator 1 (YREG1), was required to mediate between SNAIL/SLUG transcriptional activity and the upregulation of LEF-1 expression following TGF-[Beta]1/2 exposure. SNAIL/SLUG are well known transcriptional repressors [33-35], although there are a few studies which suggest that at least SNAIL can also act as a transcriptional activator [36]. In the model, we assumed the expression of SNAIL/SLUG was likely regulated by AP1/SP1 [37]. Thus, upon receiving direct SNAIL/SLUG and TGF-[Beta]3 signals, the model predicted enhanced SNAIL/SLUG expression, consistent with experimental observations. TGF-[Beta]1/2 stimulation also induces LEF-1 expression. However, literature evidence suggested that LEF-1 expression was not strongly dependent upon AP1/SP1 activity [38]. Thus, either SNAIL/SLUG are acting as inducers (contrary to substantial biochemical evidence) or, they are repressing the expression of an intermediate repressor. Given the biochemical evidence supporting SNAIL/SLUG as repressors, we created the hypothetical YREG1 repressor whose expression is downregulated by SNAIL/SLUG. The literature data therefore suggested that YREG1 had two transcriptional targets, LEF-1 and TGF-[Beta]3. By adding this regulator, our simulations became consistent with training and literature data. Medici et al. suggested that feedback between [Beta]-catenin and LEF-1 was likely, although this feedback had yet to be identified [7]. Low levels of YREG1 expression were present in all simulations to regulate the formation of the [Beta]-catenin-LEF-1 complex. To test the effect of YREG1 on the epithelial population, we conducted over-expression and knockdown simulations on untreated cells (Fig 4B and 4C). In the absence of YREG1, the population of models failed to consistently to retain a stable epithelial state (Fig 4B). Conversely, YREG1 amplification revealed an enhanced epithelial phenotype, while some inherently transformed cells moved towards a hybrid phenotype (Fig 4C). Elevated YREG1 repressed LEF-1 and TGF-[Beta]3 expression, thereby not allowing free [Beta]-catenin to form the [Beta]-catenin-LEF-1 complex, or TGF-[Beta]3 induced SMAD activation. Taken together, low YREG1 expression was required for the maintenance of a stable epithelial phenotype that was simultaneously inducible across TGF-[Beta]1/2, TGF-[Beta]3 and SNAIl/SLUG transfection.

Combined TGF-[Beta]2 and VEGF-A exposure drives heterogeneity in MCF10A and DLD1 cells

The EMT model simulations suggested the transcriptional activity of NFATc and Sp1 could be independently tuned to generate a hybrid cell population with both epithelial and mesenchymal characteristics. To test this hypothesis, we exposed either quiescent epithelial (MCFA10, Fig 5) or transformed epithelial cells (DLD1, S4 Fig) to combinations of TGF-[Beta]1/2 and/or VEGF-A. As expected, TGF-[Beta]1/2 treatment (10ng/ml) increased Slug and Vimentin expression, while repressing E-cadherin expression both at the transcript and protein levels in MCF10A (Fig 5A and 5B) and DLD1 cells (S5C Fig). Both MCF10A (Fig 5C) and DLD1 cells (S4E and S4G Fig) transitioned from quiescent cobblestone morphology to spread spindle shapes, consistent with EMT. As predicted, we found increased nuclear localization of phosphorylated Sp1 following TGF-[Beta]1/2 stimulation in both MCF10A (Fig 5B and 5C) and DLD1 cells (S4E and S4F Fig). Consistent with model predictions, VEGF-A (50ng/ml) treatment increased the abundance of NFATc1 and E-cadherin at both the transcript and protein level in both MCF10A (Fig 5A) and DLD1 cells (S4 Fig). We also found that NFATc1 nuclear localization significantly increased in both MCF10 (Fig 5B and 5C) and DLD1 (S4C and S4E Fig) cells treated with VEGF-A. Interestingly, combining VEGF-A (50ng/ml) with TGF-[Beta]1/2 (10ng/ml) resulted in significantly elevated expression of both E-cadherin and Vimentin at the transcript and protein levels in both MCF10A (Fig 5A and 5B) and DLD1 cells (S4 and S5 Figs). NFATc1 expression increased, while Sp1 expression was similar to the TGF-[Beta]1/2 case alone (Fig 5A and 5B, S4D and S4E Fig), supporting their independent regulation. The expression of Slug, and Vimentin significantly increased, while E-cadherin levels were increased in MCF10A cells (Fig 5A) and maintained at control levels in DLD1 cells (S4D Fig). As predicted, nuclear co-localization of both NFATc1 and phosphorylated Sp1 were apparent in MCF10A (Fig 5B and 5C) and DLD1 (S4E and S4F Fig) cells treated with both ligands. Taken together, combined VEGF-A and TGF-[Beta]1/2 treatment elicited a hybrid phenotype expressing both mesenchymal and epithelial characteristics in both MCF10A and DLD1 cells. This phenotype was driven by the transcriptional activity of two key transcription factors, Sp1 and NFATc, which could be modulated independently by TGF-[Beta]1/2 and VEGF-A exposure.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 5. Simultaneous TGF-[Beta]1/2 and VEGF-A treatment induced phenotype heterogeneity and is dependent upon NFAT activity in-vitro.

(A) In MCF10A, treatment with (10ng/ml) TGF-[Beta]2 increased Slug and vimentin, while ecadherin expression was inhibited at both the gene and protein level at 48 hrs. Conversely, VEGFA alone increased both NFATc1 and ecadherin gene expression. Simultaneous TGF-[Beta]2 (10ng/ml) and VEGFA (50ng/ml) treatment increased Slug, NFATc1, and vimentin expression, while also increasing ecadherin levels via qPCR. (B-C) Immunofluorescence confirmed these results and nuclear co-localization of both phospho-Sp1 and NFAT were found dependent upon TGF-[Beta]2 and VEGFA, respectively. (D) To isolate the effect of NFAT, treatment of VEGFA (50ng/ml) and VIVIT (10MM) reduced ecadherin expression at 48hrs (control-dashed line). Similarly, combined TGF-[Beta]2, VEGFA and VIVIT treatment increased Slug and vimentin expression, while inhibiting ecadherin levels via qPCR. (E) These findings were confirmed via immunofluorescence as the VIVIT peptide inhibited ecadherin and nuclear localization of NFATc1 in all three cases. (F) Quantitative flow cytometry also confirmed this trend. Similar experiments in DLD1 followed a similar trend (supplement). Magnification, 40x. Scale bars: 50Mm. C = Control, T = TGF-[Beta]2, V = VEGFA, VI = NFAT inhibitor (VIVIT). Asterisks signify statistical differences from each other according to a one-way ANOVA with Tukey's post hoc (p[precedes]0.05).

http://dx.doi.org/10.1371/journal.pcbi.1005251.g005

Our phenotypic analysis predicted that NFATc transcriptional activity was critical to maintaining E-cadherin expression in the presence of both VEGF-A and TGF-[Beta]1/2. We experimentally tested this hypothesis by exposing both MCF10A (Fig 5E and 5F) and DLD1 cells (S5 Fig) to combinations of VEGF-A and TGF-[Beta]1/2 in the presence or absence of VIVIT, a soluble peptide inhibitor of NFATc transcriptional activity [39]. Treatment with VEGF-A (50ng/ml) and VIVIT (10MM) in MCF10A cells significantly reduced E-cadherin expression compared to VEGF-A alone (Fig 5D and 5E). Co-treatment with VIVIT and TGF-[Beta]1/2 did not enhance EMT capacity of MCF10A cells above that of TGF-[Beta]1/2 alone (Fig 5A, 5B and 5E). Likewise, VIVIT in combination with both TGF-[Beta]1/2 and VEGF-A resulted in a loss of E-cadherin gene and protein expression, while Slug and Vimentin levels remained increased (Fig 5D and 5E). Quantitative flow cytometry confirmed these results in both MCF10A (Fig 5F) and DLD1 cells (S5C Fig). Both epithelial cell lines initially had high levels of E-cadherin expression, and low Vimentin abundance (Q1-99.5%), but both MCF10A and DLD1 cells shifted from an epithelial to mesenchymal phenotype (Q1-33.4%, Q4-42.8%) following TGF-[Beta]1/2 exposure. As expected, NFATc nuclear localization was repressed with VIVIT treatment regardless of ligand stimulation, while the abundance of nuclear phosphorylated Sp1 increased for both TGF-[Beta]1/2 and TGF-[Beta]1/2 + VIVIT conditions (Fig 5C and 5E). Combined TGF-[Beta]1/2 and VEGF-A increased both Vimentin and E-cadherin expression (Q1-42.1%, Q2-52.3%) compared to TGF-[Beta]1/2 alone. Together, these results demonstrate that NFATc and phosphorylated Sp1 are critical for regulating E-cadherin and Vimentin expression during phenotype heterogeneity in MCF10A and DLD1.

Ductal branching during acini formation is dependent upon phenotype heterogeneity in MCF10A and DLD1 cells

We finally employed established three-dimensional (3D) in vitro models of invasion, migration, compaction, and tubulogenesis [40] to determine the functional consequences of the hybrid phenotype (Fig 6). MCF10A and DLD1 cells were aggregated via hanging drop, placed on the surface of a collagen gel, and cultured for 72 hrs under various biochemical treatments. TGF-[Beta]1/2 stimulation significantly enhanced cell matrix invasion and matrix compaction, while in contrast VEGF-A stimulation promoted surface migration but no invasion or compaction (Fig 6B-6D). Interestingly, combined TGF-[Beta]1/2 and VEGF-A stimulation significantly increased cell migration potential above that of VEGF-A alone while maintaining 3D matrix compaction, though with decreased magnitude compared to TGF-[Beta]1/2 alone. Inhibition of NFATc transcriptional activity by VIVIT decreased migration following treatment with VEGF-A alone (Fig 6B). Co-treatment of VIVIT significantly decreased migration, while complementarily increasing invasion and compaction, when MCF10A cells were stimulated with both VEGF-A and TGF-[Beta]1/2 (Fig 6B-6D). The responses of DLD1 cells followed a similar trend to MCF10A, although the magnitudes of migration, invasion, and compaction were less. Cell circularity within 3D gels strongly and negatively correlated with both invasion and compaction regardless of treatment (Fig 6E). Circularity refers to the morphology of the cells. In general, a quiescent epithelial cells assumes a circular morphology in culture, while an active mesenchymal cell is highly elongated. The circularity index, a common means of quantifying cell morphology, relates cell area to perimeter. A perfect circle has a circularity index equal to 1.0, while a straight line has a circularity index equal to 0.0, see Butcher et al. [41]. TGF-[Beta]1/2 treatment alone resulted in irregular and spindle shaped morphology, while VEGF-A exposure promoted round quiescent cells (Fig 6A). Combined VEGF-A and TGF-[Beta]1/2 promoted morphology between these extremes. VIVIT mediated NFATc inhibition significantly reduced the circularity index, similar to TGF-[Beta]1/2 treatment (Fig 6F). VEGF-A treatment also induced the formation of tubular structures (acini), but the number of tubular branches relative to total acini was significantly increased upon combined TGF-[Beta]1/2 and VEGF-A. No tubular structures were identified within the DLD1 constructs during the 7 day tubulogenesis endpoints, supporting that MCF10A and DLD1 cells have some cell-type specific EMT sensitivity despite their underlying competency for acquiring a heterogeneous phenotype. This suggests that initial EMT sensitivity of a cell influences downstream functional response from TGF-b and VEGFA stimulation. Together, these results establish that VEGF-A and TGF-[Beta]1/2 ligand concentrations potentiate between acini and ductal branch formation in 3D culture, and are dependent upon NFATc activity.

Download:

*

PPT

PowerPoint slide

*

PNG

larger image

*

TIFF

original image

Fig 6. Ductal branching is dependent upon phenotype heterogeneity within MCF10A in 3-D culture.

MCF10A and DLD1 were formed into spheroids overnight and explanted to a collagen gel for 72 hrs. For compaction and tubular assays, cells were embedded into collagen gels for 72 hrs, and the extent of tubulogenesis was measured at 7 days. (A-D) Within MCF10A, TGF-[Beta]2 (10ng/ml) enhanced invasion and contractile properties while, VEGFA (50ng/ml) promoted increased migration. TGF-[Beta]2 with VEGFA significantly increased migration, while limiting with compaction. VIVIT (10MM) in combination with VEGFA and TGF-[Beta]2 decreased migration and compaction, while increasing invasion. (D) Likewise, cell morphology (circularity index) correlated with both invasion and compaction in MCF10A. (E-F) The size of tubular structures (acini) also increased significantly upon addition of VEGFA, while the number of ductal branches was most significant upon simultaneous TGF-[Beta]2 and VEGFA treatment (Red-Ecadherin, Green-Factin, Blue-Nuclear). DLD1 cells followed a similar trend, although the degree of migration, invasion, and compaction was less significant. In addition, no tubular structures were identified during the 7 day tubulogenesis endpoints. Scale bars: 500Mm, 1000Mm, 250Mm, and 80Mm, respectively. C = Control, T = TGF-[Beta]2, V = VEGFA, VI = NFAT inhibitor (VIVIT). Asterisks signify statistical differences from each other according to a one-way ANOVA with Tukey's post hoc (p[precedes]0.05). Boxes in the left-most panel identify regions identified by arrows that were then imaged in greater zoom in the panel immediately below. The box diagram was not repeated for arrows in the other panels for clarity, but the same method was applied.

http://dx.doi.org/10.1371/journal.pcbi.1005251.g006

Discussion

In this study, we developed a family of mathematical models describing the induction of EMT by TGF-[Beta] isoforms in the presence and absence of VEGF-A. The model, which contained 74 molecular species interconnected by 169 interactions, described the expression of 23 genes in response to growth factor stimulation. We estimated an ensemble of likely model parameters using the JuPOETs multiobjective optimization framework. The model population was trained and cross-validated to prescribe biological significance using 41 data sets generated in DLD1 colon carcinoma, MDCKII, and A375 melanoma cell lines [7]. Analysis of this population predicted possible phenotypic modes (and their associated signaling) that cells could exhibit when stimulated with TGF-[Beta] and/or VEGF-A. The most novel hypothesis generated from the analysis was that cells could operate in a hybrid state defined by both epithelial and mesenchymal traits when stimulated simultaneously with TGF-[Beta] and VEGF-A. We tested this hypothesis in MCF10A and DLD1 cells stimulated with combinations of TGF-[Beta] and VEGF-A. As expected, in the presence of TGF-[Beta] or VEGF-A alone, MCF10A and DLD1 cells were either mesenchymal or epithelial, respectively. However, with both TGF-[Beta] and VEGF-A, MCF10A and DLD1 cells exhibited a hybrid phenotype, having both epithelial and mesenchymal characteristics. Furthermore, we found that functional traits such as tubulogenesis and ductal branching were different for cells in this hybrid phenotype. Together, this study established a predictive model of EMT induction, determined that deterministic model ensembles could predict population heterogeneity, and proved the existence of a unique hybrid phenotype resulting from the simultaneous integration of extracellular growth factor signals.

Cells routinely process a multitude of signals simultaneously, especially when coordinating developmental or pathological programs. For example, oncogenic cells integrate both mechanical and chemical cues in their local microenvironment during tumorigenesis, including cytokines VEGF and TGF-[Beta] [42]. VEGF-A mediates pathological angiogenic remodeling of tumors [43], while TGF-[Beta] can elicit both protective and oncogenic responses [44, 45]. While much research has tested signaling pathways individually, far less is understood about combinatorial stimulation, such as with both VEGF-A and TGF-[Beta]. Both in vitro and in vivo studies have suggested that epithelial cells can exhibit heterogeneous phenotypes in addition to classically defined epithelial or mesenchymal states [46, 47]. For example, expression profiling in human epithelial cancer cell lines demonstrated a spectrum of phenotypes, including some that expressed both E-cadherin and Vimentin simultaneously [48, 49]. Zajchowski et al., speculated that these expression profiles were somehow important for maintaining epithelial properties, while simultaneously allowing other functional behavior such as proliferation and migration [50]. Whether and how heterogeneous phenotypes arise and participate in cancer progression, as well as their response to pharmacological inhibition are fundamental questions that should receive increased attention. In this study, we determined that a hybrid phenotype could be obtained through combined treatment with VEGF-A and TGF-[Beta], both common factors localized in the tumor microenvironment. Furthermore, our systematic simulation-experimentation strategy identified that the transcriptional activity of Sp1 and NFATc were the critical factors controlling this phenotypic heterogeneity. Several studies have highlighted the importance of NFATc as a key transcription factor involved in cell growth, survival, invasion, angiogenesis and cancer [51]. For example, proliferation and anchorage-independent growth of pancreatic tumor cells is dependent on calcineurin and NFATc1 activity, consistent with the high levels of nuclear NFATc1 found in pancreatic tumors [52]. Likewise, our results found that VEGF-A was a potent inducer of NFATc1 expression, which may be required for epithelial cell migration and tubulogenesis. Although specific NFATc isoforms were not distinguished in the model, our simulations suggested that NFATc transcriptional activity was capable of maintaining epithelial traits, even during TGF-[Beta] induced EMT. Experimentally, we found that E-cadherin expression was dependent upon NFATc dephosphorylation in response to simultaneous VEGF-A and TGF-[Beta]1/2 treatment. Thus, these results support the hypothesis that NFATc activity plays a critical role in maintaining cell-cell contacts, even during partial EMT.

Epithelial cells reproduce tissue-like organization when grown in a three-dimensional extracellular matrix (ECM) environment, and therefore are an attractive model to study morphogenic mechanisms. It is well established that MCF10A cells form structures that closely resemble acini (multi-lobed cluster of cells) in three-dimensional in vitro cultures [53]. It has been postulated that a cellular response reminiscent of partial EMT underlies this process, stimulating further branching and formation of acini [54]. Normally well controlled process such as tubulogenesis can be co-opted by cancer cells to break away from a primary lesion and invade through the surrounding stroma [55]. However, by retaining a transient hybrid EMT-like state, clusters of these tube-forming tumor cells can reform at a high rate after invasion, possibly explaining why invasive human carcinomas frequently appear to be cellular collections with varying degrees of gland-like differentiation [56]. In this study, we showed that our predicted hybrid phenotype generated by simultaneous treatment of epithelial cells with VEGF-A and TGF-[Beta] possessed altered migration and invasion, which enhanced tubular branching. A salient feature of this behavior, however, was the retention of cell-cell contacts that allowed cells to migrate without completely dissociating from their neighbors. Thus, our results support a mechanism in which hybrid cells can maintain some functional characteristics of epithelial cells such as cell-cell adhesion, which are normally lost in a fully differentiated mesenchymal state. The tumor microenvironment contains many soluble signals simultaneously, including VEGF and TGF-[Beta]. Thus, it is likely that some cancerous epithelial cells could exhibit hybrid EMT phenotypic states. This may explain why fibroblastoid morphology, a classical feature of EMT, is not commonly observed in human carcinomas [56]. This study focused on the combinatorial effects of two very different ligand families present together in the tumor environment. Additional modeling studies are required to unravel the global response of epithelial cells to the full spectrum of chemical, substrate, and mechanical cues. The simulation strategy presented here is readily adaptable to larger species sets, with the major advantage that experimentally testable hypotheses can be generated regarding how signals get integrated to produce global cellular response. Furthermore, by simulating multiple ensembles of parameter sets, subpopulations across a constellation of phenotypes can be created and mined for common and/or divergent signaling characteristics. This is a significant advantage over forced convergence to a single unique solution and thereby generating a potentially non-physiological homogeneous population.

The deterministic population of EMT models predicted heterogeneous behavior that was qualitatively consistent with experimental studies. There is a diversity of algorithmic approaches to estimate model parameters [57], as well as many strategies to integrate model identification with experimental design [58, 59]. However, despite these advances, the identification of models describing intracellular network behavior remains challenging. There are different schools of thought to deal with this challenge. One school has focused on model reduction. Data-driven approaches [16], boolean [60] or other logical model formulations [17, 61] are emerging paradigms that constrain model complexity by the availability of the training and validation data. Other techniques such as constraints based modeling, which is commonly used to model metabolic networks, have also been applied to model transcriptional networks, although primarily in lower eukaryotes and prokaryotes [62]. These techniques (and many others, see review [63]) are certainly exciting, with many interesting properties. Here, we used a traditional approach of mass action kinetics within an ordinary differential equation framework that also included transfer functions to simplify scenarios where reactions involving one species are controlled by several others (e.g. E-cadherin transcription). The identification problem for the EMT model was underdetermined (not uncommon for differential equation based models). However, a central criticism leveled by biologists is that model simplification is often done at the cost of biological reality, or done for reasons of computational expediency [64]. To avoid this criticism, we systematically identified an ensemble of likely models each consistent with the training data, instead of a single but uncertain best fit model. Previously, we (and others) have suggested that deterministic ensembles could model heterogeneous populations in situations where stochastic computation was not feasible [21]. Population heterogeneity using deterministic model families has previously been explored for bacterial growth in batch cultures [65]. In that case, distributions were generated because the model parameters varied over the ensemble, i.e., extrinsic noise led to population heterogeneity. In this study, parameters controlling physical interactions such as disassociation rates, or processes such as gene expression were distributed over the ensemble. Population heterogeneity can also arise from intrinsic thermal fluctuations, which are not captured by a deterministic population of models [66]. Thus, deterministic ensembles, provide a coarse-grained or extrinsic-only ability to simulate population diversity. Despite this limitation, our prediction of phenotypic heterogeneity (and the underlying signaling events responsible for the heterogeneity) was consistent with experimental observations. This suggested that deterministic ensembles could simulate disease or developmental processes in which heterogeneity plays an important role, without having to resort to stochastic simulation.

A common criticism of ODE modeling has been the poorly characterized effect of structural and parametric uncertainty. In this study, parametric uncertainty was addressed by developing an ensemble of probable models instead of a single best-fit but uncertain model using multiobjective optimization. While computationally complex, multiobjective optimization is an important tool to address qualitative conflicts in training data that arise from experimental error or cell line artifacts [67]. On the other hand, structural uncertainty is defined as uncertainty in the biological connectivity. The EMT model connectivity was assembled from an extensive literature review. However, several potentially important signaling mechanisms were not included. First, we identified a potential gap in biological knowledge surrounding the regulation of LEF-1 expression, that was filled by the addition of the hypothetical YREG1 transcriptional repressor. The LEF-1 transcription factor is expressed in tissues that undergo EMT during embryogenesis [68, 69], and has been suggested to promote an invasive phenotype in cancer cells [8, 70]. Low levels of YREG1 were important for stabilizing the interaction between LEF-1 and [Beta]-catenin, while elevated levels inhibited EMT by downregulating LEF-1 transcriptional activity. Recent evidence has established a complex role of Amino terminal Enhancer of Split (AES) and Groucho/TLE on suppressing LEF-1 activity. AES opposes LEF-1 transcriptional activation while Groucho/TLE binds with LEF-1 for a histone deacetylase repression. In addition, [Beta]-catenin directly displaces Groucho/TLE repressors from TCF/LEF-1 in Wnt-mediated transcription activation [71, 72]. Our model agrees with this newly discovered feedback system, as YREG1 regulates LEF-1 activity leading to EMT stabilization.

NF-KB may also play an essential role of in the epithelial transformation. NF-KB may influence Snail expression through the AKT pathway and directly stabilize Snail activity [73]. This is particularly important for integrating inflammation pathways, such as interleukin-6 (IL-6) and tumor necrosis factor-[alpha] (TNF-[alpha] ), which have been linked to EMT in pathological conditions [74]. Other pathways such as Notch have also been shown to act synergistically with TGF-[Beta] to express Slug in the developing embryo [75]. In addition, other regulatory proteins such as ZEB1 could be added to the next generation model. ZEB1 can activate the expression of mesenchymal genes such as collagens, smooth muscle actin and myosin, Vimentin while simultaneously contributing to the repression of epithelial genes [76]. Therefore, ZEB1 can act in parallel to SNAIL and is an alternative to the canonical signaling paradigm presented here. Hong et al. developed an ODE model of ZEB1 signaling in MCF10A cells which also demonstrated intermediate EMT states [77]. Thus, there may be alternative modes resulting in a hybrid phenotype that should be included in future model architecture. Lastly, while we have modeled classical protein signaling, we have not considered the role of regulatory RNAs on EMT. There is growing evidence that microRNAs (miRNAs) play a strong role in EMT, where several miRNAs, for example miR-21 and miR-31 are strongly associated with TGF-[Beta] exposure [78], and have suggested that the EMT process is driven by a miR-34/SNAIL and/or miR-200/ZEB feedback circuit [79-81]. Addressing missing structural components like these (and potentially others), could generate more insight into TGF-[Beta] signaling and its role in phenotypic transformation.

Materials and Methods

The model code and parameter ensemble is freely available under an MIT software license and can be downloaded from http://www.varnerlab.org.

Signaling network connectivity

The EMT model described the gene expression program resulting from TGF-[Beta] and VEGF-A signaling in a prototypical epithelial cell. The TGF-[Beta]-EMT network contained 97 nodes (proteins, mRNA or genes) and 169 interactions. The network connectivity was curated from more than 40 primary literature sources in combination with on-line databases [82, 83]. The model interactome was not specific to a single epithelial cell line. Rather, we assembled canonical pathways involved in TGF-[Beta] and VEGF-A signaling, defaulting to human connectivity when possible. Using a canonical architecture allowed us to explore general features of TGF-[Beta] induced EMT without cell line specific artifacts.

Our signaling network reconstruction was based on Medici et al. who identified the pathways though which MDCKII, DLD1 colon carcinoma, and A375 melanoma cells transition towards a mesenchymal phenotype [7]. Sequential activation of MAPK and Smad pathways were initiated upon addition of TGF-[Beta]1/2. Briefly, TGF-[Beta]2 signals through the RAS-RAF-MEK-ERK pathway to up-regulate Snail and Slug expression [6]. Snail, a known repressor of junctional proteins, inhibits the expression of E-cadherin [70]. This initial repression of E-cadherin leads to a release of [Beta]-catenin from the cell membrane. This release of [Beta]-catenin can then translocate to the nucleus and form transcriptional complexes with TCF-4 to drive TGF-[Beta]3 expression [7]. The PI3K to GSK3 pathway was included and acted as an activating mechanism of [Beta]-catenin signaling through TGF-[Beta]3 signaling [7]. GSK3 is known to act on [Beta]-catenin signaling through the ubiquitin-proteasome pathway [84, 85]. Thereby, further [Beta]-catenin release also resulted from by TGF-[Beta]3 signals to the cells interior by binding to type II receptors, which form heterodimers with type I receptors (ALK5) [86]. This activates the receptors serine/threonine kinase activity to phosphorylate and activate the receptor Smads 2/3 [87]. In the model, receptors are simplified and represented as either bound or unbound complexes with their ligands. Phosphorylated Smads 2/3 (pSmad2/3) form heterodimers and translocate to the nucleus. pSmads complexes up-regulate other transcription factors, such as LEF-1. The pSmad2/4-LEF-1 complex has been shown to directly repress the E-cadherin gene [88]. LEF-1 also binds with [Beta]-catenin to upregulate mesenchymal proteins such as fibronectin [89]. In the model, Smad signaling is consolidated into a single Smad species that can act in a co-dependent fashion with LEF1 to downregulate E-cadherin via a transfer function, eliminating the need for an explicity LEF-1, pSmad complex. The EMT gene expression program was initiated by the binding of TGF-[Beta] isoforms to TGF-[Beta] surface receptors, starting the downstream signaling program. Repression of E-cadherin expression is the central event in the transition from an epithelial to a mesenchymal phenotype [70]. However, this transition is not solely driven by transcriptional events. At the protein level, the repression of E-cadherin leads to a release of [Beta]-catenin from cell membrane. Cytosolic [Beta]-catenin then translocates to the nucleus and forms transcriptionally-active complexes with immunoglobulin transcription factor 2 (TCF-4) to drive TGF-[Beta]3 expression [7]. The PI3K to GSK3 pathway was included and acted as an activating mechanism of [Beta]-catenin signaling through TGF-[Beta]3 signaling [7]. GSK3 is known to act on [Beta]-catenin signaling through APC complex associated ubiquitin-proteasome pathway. The APC complex is represented in our model and serves as a second reservoir of [Beta]-catenin in untransformed cells whose sequestration is regulated by GSK3 [7, 84, 85]. Lastly, VEGF-A activation of NFATc1 takes place through calcineurin signaling leading to an enhancement of E-cadherin expression [90], as supported by our VEGF-A experimental data (S1 Fig).

Formulation, solution and analysis of the EMT model equations

The EMT model integrated kinetic ODE modeling with a rule-based description of signaling activity and gene expression regulation. Regulation of enzyme activity or gene expression processes was modeled by combining individual regulatory contributions, quantified by transfer-functions, into a control coefficients using an integration rule. These control coefficients then modified the kinetic rates that appeared in the model. This strategy is similar in spirit to the Constrained Fuzzy Logic (cFL) approach of Lauffenburger and coworkers which has been used to effectively model signal transduction pathways important in human health [17]. In our formulation, Hill-like transfer functions 0 ≤ f(x) ≤ 1 were used to calculate the influence of factor abundance upon enzyme activity or a gene expression process. In this context, factors can be individual metabolite or protein levels or some function of abundance, e.g., the product of metabolite or proteins levels. However, more generally, factors can also correspond to non-modeled influences, categorial variables or other abstract quantities. In the current study, we let factors correspond to the abundance of individual proteins, inhibitors or transcription factors, as well as products of transcription factors where it is necessary to represent a co-dependent regulatory activity. When a protein or gene was potentially sensitive to more than one regulatory input, logical integration rules were used to select which regulatory transfer function influenced enzyme activity at any given time. Thus, our EMT network model encoded complex signaling and regulatory features with a relatively small number of equations. This model formulation has been used previously to construct reduced order kinetic models of synthetic RNA circuits [91], cell free metabolic models [92] and reduced-order models of blood coagulation [93]. In this study, we extended this approach to describe gene expression processes in addition to signal transduction events.

EMT signaling events.

EMT signaling events were modeled using either saturation or mass-action kinetics within an ordinary differential equation (ODE) framework: (1) where denotes the number of signaling reactions and denotes the number of signaling proteins in the model. The quantity [tau]i denotes a time scale parameter for species i which captures un-modeled effects; in the current study [tau]i = 1 for all species. The quantity rj(x, [epsilon], k) denotes the rate of reaction j. Typically, reaction j is a non-linear function of biochemical and enzyme species abundance, as well as unknown model parameters k (). The quantity [sigma]ij denotes the stoichiometric coefficient for species i in reaction j. If [sigma]ij > 0, species i is produced by reaction j. Conversely, if [sigma]ij < 0, species i is consumed by reaction j, while [sigma]ij = 0 indicates species i is not connected with reaction j. Species balances were subject to the initial conditions x(to) = xo.

Rate processes were written as the product of a kinetic term () and a control term (vj) in the EMT model. The rate of enzyme catalyzed reactions was modeled using saturation kinetics: (2) where denotes the catalytic rate constant for reaction j, [epsilon]i denotes the abundance of the enzyme catalyzing reaction j, and Kjs denotes the saturation constant for species s in reaction j. On the other hand, mass action kinetics were used to model protein-protein binding interactions within the network: (3) where denotes the maximum rate for reaction j, [sigma]sj denotes the stoichiometric coefficient for species s in reaction j, and s ∈ mj denotes the set of reactants for reaction j. Reversible binding was decomposed into two irreversible steps.

The control terms 0 ≤ vj ≤ 1 depended upon the combination of factors which influenced rate process j. For each rate, we used a rule-based approach to select from competing control factors. If rate j was influenced by 1, ..., m factors, we modeled this relationship as where 0 ≤ fij([dot]) ≤ 1 denotes a regulatory transfer function quantifying the influence of factor i on rate j. The function is an integration rule which maps the output of regulatory transfer functions into a control variable. In this study, we used and hill transfer functions [92, 93]. If a process had no modifying factors, vj = 1.

EMT gene expression processes.

The EMT model described both signal transduction and gene expression events following the addition of TGF-[Beta] and VEGF-A. For each gene of the we considered, we modeled both the resulting mRNA (mj) and protein (pj): (4) (5) where . The terms rT,j and rX,j denote the specific rates of transcription, and translation while the terms [theta]

m,j and [theta]

p,j denote first-order degradation constants for mRNA and protein, respectively. Lastly, M denotes the specific growth rate, and [Lambda]j denotes the constitutive rate of gene expression for gene j. The specific transcription rate was modeled as the product of a kinetic term and a control term uj which described how the abundance of transcription factors, or other regulators influenced the expression of gene j. The kinetic rate of transcription was modeled as: (6) where the maximum gene expression rate was defined as the product of a characteristic transcription rate constant (kT) and the abundance of RNA polymerase, . The parameter [alpha] j denotes the gene specific correction to the characteristic transcription rate (estimated in this study). Similar to the signaling processes, the gene expression control term 0 ≤ uj ≤ 1 depended upon the combination of factors which influenced rate process j. For each rate, we used a rule-based approach to select from competing control factors. If the expression of gene j was influenced by 1, ..., m factors, we modeled this relationship as where 0 ≤ fij([dot]) ≤ 1 denotes a regulatory transfer function quantifying the influence of factor i on the expression of gene j. The function is an integration rule which maps the output of regulatory transfer functions into a control variable. In this study, we used and hill transfer functions [92, 93]. If a gene expression process has no modifying factors, uj = 1. The degradation rate constants were defined as [theta]

m,j = [delta]jkd,m and [theta]

p,j = [gamma]jkd,p where kd,m and kd,p denote characteristic degradation constants for mRNA and protein, respectively and [delta]j and [gamma]j denote transcript and protein specific correction factors to these constants.

Lastly, the specific translation rate was modeled as: (7) where denotes a characteristic maximum translation rate estimated from literature, [Beta]j denotes the transcript specific correction the characteristic translation rate, and KX denotes a translation saturation constant. The characteristic maximum translation rate was defined as the product of a characteristic translation rate constant (kX) and the abundance of ribosomes (RIBO), .

In this study, we estimated kT, kX, kd,m, kd,p, RNAP and RIBO directly from literature using estimates of transcription and translation rates, the half-life of a typical mRNA and protein and a typical value for the copies per cell of RNA polymerase and ribosomes. The saturation constants KX and KT were adjusted so that gene expression and translation resulted in gene products on a biologically realistic concentration scale. Next, we calculated the concentration for gene j by assuming, on average, that a cell had two copies of each gene at any given time. Thus, the bulk of our gene expression parameters were based directly upon literature values, and were not adjusted during model identification. The values used for the characteristic transcription/translation parameters, degradation constants and macromolecular copy number are given in the supplement along with the specific formulas required to calculate all derived constants.

The signaling and gene expression model equations were implemented in Julia and solved using the CVODE routine of the Sundials package [94, 95]. The model code and parameter ensemble is freely available under an MIT software license and can be downloaded from http://www.varnerlab.org.

Estimation of model parameters using multiobjective optimization.

We estimated a population of likely EMT model parameter sets (each consistent with the training data) using 41 data sets generated in DLD1 colon carcinoma, MDCKII, and A375 melanoma cells taken from Medici et al. [7]. In addition to rate and saturation constants appearing in the signal transduction equations, we estimate the transcript/protein specific correction factors [alpha] j, [Beta]j, [delta]j and [delta]j and control coefficients from these data. We used the Pareto Optimal Ensemble Technique (JuPOETs) multiobjective optimization framework in combination with leave-one-out cross-validation to estimate an ensemble of model parameters [20, 96]. Model parameter values were adjusted to minimize the residual between simulations and experimental measurements. Cross-validation was used to calculate both training and prediction error during the parameter estimation procedure [97]. The 41 intracellular protein and mRNA data-sets used for identification were organized into 11 objective functions. These 11 objective functions were then partitioned, where each partition contained ten training objectives and one validation objective. The training and validation data were Western blots. We achieved a biologically realistic concentration scale by establishing characteristic rates of transcription, translation, mRNA and protein degradation, as well as characteristic concentrations of ribosomes and RNAPs using the Bionumbers database [31], (supplemental materials). The overall concentration scale was nM, with proteins ranging from 10-1000nM and mRNA ranging from 0.01 to 1nM, reflecting the true abundances and ratios between each species. An initial nominal parameter set was established by inspection. JuPOETs was then allowed to search in a neighborhood of ±30% of this nominal set. The correction factors estimated by JuPOETs were in the range of 0.1-10, with the mean value of 1.68; thus, the characteristic rate parameters estimated from literature were consistent with measurements. The parameter ensemble estimated by JuPOETs is available with the model source code. JuPOETs is open source and freely available for download under an MIT software license from http://www.varnerlab.org.

Cell culture and experimental interrogation

DLD1 colon carcinoma and MCF10A cells were acquired from the American Tissue Culture Collection (Manassas, VA). Cells were grown in culture with RPMI-1640 medium with 10% fetal bovine serum and 1% penicillin/streptomycin or MGEM-2 supplemented with insulin, bovine pituitary extract, cholera toxin, hEGF, hydrocortisone, 5% horse serum, and 1% penicillin/streptomycin, respectively. Cells were passaged 1:3 or 1:4 every 3-6 d and used between passages 4 and 8. During treatment regime, DLD1 and MCF10A cells aggregated overnight in hanging drop culture (20ML; 20,000 cells). The spherical aggregates were placed on the surface of neutralized type I collagen hydrogels (1.5mg/mL) and allowed to adhere. Cultures were then serum starved (1% serum) for 24 hours. Recombinant human TGF-[Beta]2 (R & D Systems, Minneapolis, MN) was added to the culture medium at a concentration of 10 ng/ml and recombinant VEGFA165 at a concentration of (5ng/ml, 50ng/ml) for all relative experiments. NFAT inhibitor (VIVIT peptide) (EMDBiosciences, Darmstadt, Germany), was added to the culture medium at a concentration of 10MM for all relative experiments. mRNA was then harvested at both the 3hr and 48hr time-point.

RT-PCR.

RNA extractions were performed using a Qiagen total RNA purification kit (Qiagen, Valencia, CA) and RNA was reverse transcribed to cDNA using the SuperScript III RT-PCR kit with oligo(dT) primer (Invitrogen). Sufficient quality RNA was determined by an absorbance ratio A260/A280 of 1.8-2.1, while the quantity of RNA was determined by measuring the absorbance at 260nm (A260). Real-time PCR experiments were conducted using the SYBR Green PCR system (Biorad, Hercules, CA) on a Biorad CFX96 cycler, with 40 cycles per sample. Cycling temperatures were as follows: denaturing, 95C; annealing, 60C; and extension, 70C. Primers were designed to detect GAPDH, E-cadherin, vimentin, Slug, Sp1, and NFATc1 in cDNA clones: Sp1 (F-TTG AAA AAG GAG TTG GTG GC, R-TGC TGG TTC TGT AAG TTG GG, Accession NG030361.1), NFATc1 (F-GCA TCA CAG GGA AGA CCG TGT C, R-GAA GTT CAA TGT CGG AGT TTC TGA G, Accession NG029226.1). GAPDH, E-cadherin, vimentin, and Slug primers were taken from previously published literature [7].

Antibody staining.

Samples were fixed in 4% PFA overnight at 4C. Samples were then washed for 15 minutes on a rocker 3 times with PBS, permeabilized with 0.2% Triton-X 100 (VWR International, Radnor, PA) for 10 minutes, and washed another 3 times with PBS. Samples were incubated overnight at 4C in a 1% BSA (Rockland Immunochemicals, Inc., Gilbertsville, PA) blocking solution followed by another 4C overnight incubation with either rabbit anti-human E-cadherin 1:100 (Abcam, ab53033), mouse anti-human phospho-Sp1 1:100 (Abcam, ab37707), mouse anti-human vimentin 1:100 (Invitrogen, V9), and rabbit anti-human NFATc1 (Santa Cruz, sc-7294) 1:100. After 3 washes for 15 minutes with PBS, samples were exposed to Alexa Fluor 488 or 568 conjugated (Invitrogen), species specific secondary antibodies at 1:100 in 1% BSA for 2 hours at room temperature. Three more washes with PBS for 15 minutes were followed by incubation with either DRAQ5 far red nuclear stain (Enzo Life Sciences, Plymouth Meeting, PA) at 1:1000.

FACS.

Flow cytometry for E-cadherin 1:100 (Abcam) and vimentin 1:100 expressing cells was performed. Briefly, cells were trypsinized, fixed with 4% PFA for 10 min and then preserved in 50% methanol/PBS. Cells were kept in the -20C until antibody staining was preformed. Samples were divided into multiple aliquots in order to stain the proteins separately and compensate for secondary antibody non-specific binding. Cells were incubated for 24 hrs at 4 C in primary antibody diluted in either PBS (extracellular) or 0.2% saponin-PBS (intracellular). Cells were then washed 3 times with PBS and incubated with appropriate secondary antibodies and imaged using a Coulter Epics XL-MCL Flow Cytometer (Coulter). All samples were compensated using appropriate background subtraction and all samples were normalized using 7500 cells per flow condition.

Three-dimensional culture and tubulogenesis assays.

For invasion/migration assays, cells were resuspended in culture media, and allowed to aggregate overnight in hanging drop culture (20ML; 20,000 cells). The spherical aggregates were placed on the surface of neutralized type I collagen hydrogels (1.5mg/mL) and allowed to adhere for 2 hrs before adding treatments. Cultures were maintained for 72 hrs, after which they were fixed in 4% PFA and slowly rehydrated using PBS. For compaction assays, cells were pelleted via centrifugation and resuspended within a neutralized collagen hydrogel (1.5mg/mL) solution at a density of 400,000 cells/mL. 250ML of gel was inoculated into culture wells, which solidified after 60min. Treatments were then added within 800ML of the culture medium without serum. Gels were liberated from the surfaces of the culture wells the next day and cultured free floating for an additional 3-7 days, exchanging serum free media with appropriate factors every 48 hrs.

Tubulogenesis was defined as a typical nonmalignant acini structure. This includes a polarized epithelial cell, hollow lumen, and the basal sides of the cell are surrounded by ECM proteins (Fig 6A, Controls or VEGF treated). Previous work has shown that change in the morphological characteristics of nontumorigenic MCF10A epithelial acini occur over time and exploiting them to growth in 3D culture can be quantified. For example, using image segmentation, Chang et al. [98] examined the elongation of the MCF10A acini at 6, 12, and 96 hours after a particular treatment. Polizzotti et al. [99] also suggested a computational method to quantify acini structure based on morphological characteristics in nonmalignant, noninvasive, and invasive conditions. Adapted from these approaches, we first fluorescently labeled our cultures and captured the acini structures by 3D confocal microscopy. Next individual acini structures in the images were segmented by imageJ and labeled. We then extracted the number of ductal branches. Ductal branching was defined as any elongated cell cluster extending away from the total acini structure, which was manually segmented and counted using ImageJ. A total of 5 images for each condition were used, and approximately 12 acini were analyzed in each image. Total branching was normalized to the amount of acini present, and provides an overall general assessment to the extent of acini remodeling.

Statistics.

Results are expressed as mean ± standard error, n≥6. Data was analyzed with the GraphPad Prism version 4.00 for Windows (GraphPad Software, San Diego, CA) and SAS (Statistical Analysis Software, Cary, NC). A one-way ANOVA with Tukey's post hoc was used to compare differences between means and data was transformed when necessary to obtain equal sample variances. Differences between means were considered significant at p[precedes]0.05.

Supporting Information

S1 Text. Supplemental Materials and Methods.

(PDF)

S1 Fig. VEGF-A qPCR data used to hand fit VEGF enhancement of E-cadherin expression.

mRNA was harvested after 3hr and 48hr timepoint.

(TIF)

S2 Fig. We fit an additional objective demonstrating the activation of GSK3 through PI3K.

The model captured this activation through TGF-[Beta]3 signaling. LY294002 is a PI3K inhibitor.

(TIF)

S3 Fig. Training and prediction values as a function of condition for the 11 TGF-[Beta] objective functions versus a random parameter control.

(TIF)

S4 Fig. VEGF-A attenuates TGF-[Beta]1/2 to induce phenotype heterogeneity in DLD1.

(A) In DLD1, we found that 5ng/ml of VEGFA increased NFATc1 and E-cadherin gene expression via qPCR and 50ng/ml potentiated this effect at 48 hrs. (B-C) These findings were confirmed at the protein level via immunofluorescence, as ecadherin levels and nuclear localization of NFATc1 increased. (D) Treatment with (10ng/ml) TGF[Beta]2 resulted in mesenchymal transformation as measured via qPCR against target genes Slug, ecadherin, vimentin, Sp1, and NFATc1. (E-F) Immunofluorescence and nuclear localization revealed a strong presence of phospho-Sp1. (G) Combination of VEGFA (50ng/ml) and TGF[Beta]2 (10ng/ml) treatment resulted in increased Slug, NFATc1, and vimentin expression, while also increasing ecadherin levels compared to control. (H) Immunofluorescence confirmed these results, as both ecadherin and vimentin levels were elevated. (I) A significant increase in nuclear localization of both NFATc1 and phospho-Sp1 were also found. Magnification, 40x. Scale bars: 50Mm. C = Control, T = TGF[Beta]2, V = VEGFA, VI = NFAT inhibitor (VIVIT). Asterisks signify statistical differences from each other according to a one-way ANOVA with Tukey's post hoc (p[precedes]0.05).

(TIF)

S5 Fig. E-cadherin expression is dependent upon NFAT activity in DLD1.

(A) Treatment with VEGFA (50ng/ml) and NFAT inhibitory peptide VIVIT (10MM) resulted in significantly reduced ecadherin expression (qRT-PCR at 48hrs). Addition of TGF[Beta]2 (10ng/ml) and VIVIT resulted in increased Slug and vimentin expression, while inhibiting ecadherin levels. Combined TGF[Beta]2, VEGFA, and VIVIT treatment resulted in target genes Slug and vimentin expression increased, while inhibiting ecadherin levels. No change in Sp1 or NFATc1 expression was found. (B) These findings were confirmed via immunofluorescence as the VIVIT inhibitors was shown to inhibit ecadherin levels in all three cases. We also found no change in gene or nuclear localization of NFATc1 in all three cases, while phospho-Sp1 was found to increase in both TGF[Beta] conditions. (C) Quantitative flow cytometry also confirmed this trend. (D,E) TGF[Beta]2, VEGFA and VIVIT treatment in DLD1 and MCF10A resulted in no change of Sp1 expression or NFATc1 expression. (F) Likewise, no change in nuclear localization of NFAT in all three cases, however phospho-Sp1 was found to increase in both TGF[Beta] conditions. Magnification, 40x. Scale bars: 50Mm. C = Control, T = TGF[Beta]2, V = VEGFA, VI = NFAT inhibitor (VIVIT). Asterisks signify statistical differences from each other according to a one-way ANOVA with Tukey's post hoc (p[precedes]0.05).

(TIF)

Author Contributions

1. Conceptualization: RG JDV JB.

2. Data curation: RG DMB AC JDV JB.

3. Formal analysis: RG DMB JDV JB.

4. Funding acquisition: JDV JB.

5. Investigation: RG DMB.

6. Methodology: RG AC DMB JDV JB.

7. Project administration: JDV JB.

8. Resources: JDV JB.

9. Software: DMB AC JDV.

10. Supervision: JDV JB.

11. Validation: RG DMB JDV JB.

12. Visualization: RG DMB JDV JB.

13. Writing - original draft: RG DMB JDV JB.

14. Writing - review & editing: RG DMB AC JDV JB.

References

Thiery JP. Epithelial-mesenchymal transitions in development and pathologies. Curr Opin Cell Biol. 2003;15(6):740-6. doi: 10.1016/j.ceb.2003.10.006. pmid:14644200

Stahl PJ, Felsen D. Transforming growth factor-beta, basement membrane, and epithelial-mesenchymal transdifferentiation: implications for fibrosis in kidney disease. Am J Pathol. 2001;159(4):1187-92. doi: 10.1016/S0002-9440(10)62503-5. pmid:11583944

Eisenberg LM, Markwald RR. Molecular regulation of atrioventricular valvuloseptal morphogenesis. Circ Res. 1995;77(1):1-6. doi: 10.1161/01.RES.77.1.1. pmid:7788867

Zavadil J, Böttinger EP. TGF-beta and epithelial-to-mesenchymal transitions. Oncogene. 2005;24(37):5764-74. doi: 10.1038/sj.onc.1208927. pmid:16123809

Xu J, Lamouille S, Derynck R. TGF-beta-induced epithelial to mesenchymal transition. Cell Res. 2009;19(2):156-72. doi: 10.1038/cr.2009.5. pmid:19153598

Medici D, Hay ED, Goodenough DA. Cooperation between snail and LEF-1 transcription factors is essential for TGF-beta1-induced epithelial-mesenchymal transition. Mol Biol Cell. 2006;17(4):1871-9. doi: 10.1091/mbc.E05-08-0767. pmid:16467384

Medici D, Hay ED, Olsen BR. Snail and Slug promote epithelial-mesenchymal transition through beta-catenin-T-cell factor-4-dependent expression of transforming growth factor-beta3. Mol Biol Cell. 2008;19(11):4875-87. doi: 10.1091/mbc.E08-05-0506. pmid:18799618

Kim K, Lu Z, Hay ED. Direct evidence for a role of beta-catenin/LEF-1 signaling pathway in induction of EMT. Cell Biol Int. 2002;26(5):463-76. doi: 10.1006/cbir.2002.0901. pmid:12095232

Jiang YG, Luo Y, He Dl, Li X, Zhang Ll, Peng T, et al. Role of Wnt/beta-catenin signaling pathway in epithelial-mesenchymal transition of human prostate cancer induced by hypoxia-inducible factor-1alpha. Int J Urol. 2007;14(11):1034-9. doi: 10.1111/j.1442-2042.2007.01866.x. pmid:17956532

Huber MA, Azoitei N, Baumann B, Grünert S, Sommer A, Pehamberger H, et al. NF-kappaB is essential for epithelial-mesenchymal transition and metastasis in a model of breast cancer progression. J Clin Invest. 2004;114(4):569-81. doi: 10.1172/JCI21358. pmid:15314694

Hardy KM, Booth BW, Hendrix MJC, Salomon DS, Strizzi L. ErbB/EGF signaling and EMT in mammary development and breast cancer. J Mammary Gland Biol Neoplasia. 2010;15(2):191-9. doi: 10.1007/s10911-010-9172-2. pmid:20369376

Ahmed S, Nawshad A. Complexity in interpretation of embryonic epithelial-mesenchymal transition in response to transforming growth factor-beta signaling. Cells Tissues Organs. 2007;185(1-3):131-45. doi: 10.1159/000101314. pmid:17587819

Chung SW, Miles FL, Sikes RA, Cooper CR, Farach-Carson MC, Ogunnaike BA. Quantitative modeling and analysis of the transforming growth factor beta signaling pathway. Biophys J. 2009;96(5):1733-50. doi: 10.1016/j.bpj.2008.11.050. pmid:19254534

Vilar JMG, Jansen R, Sander C. Signal processing in the TGF-beta superfamily ligand-receptor network. PLoS Comput Biol. 2006;2(1):e3. doi: 10.1371/journal.pcbi.0020003. pmid:16446785

Steinway SN, Zañudo JGT, Ding W, Rountree CB, Feith DJ, Loughran TP Jr, et al. Network modeling of TGFB signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation. Cancer Res. 2014;74(21):5963-77. doi: 10.1158/0008-5472.CAN-14-0225. pmid:25189528

Cirit M, Haugh JM. Data-driven modelling of receptor tyrosine kinase signalling networks quantifies receptor-specific potencies of PI3K- and Ras-dependent ERK activation. Biochem J. 2012;441(1):77-85. doi: 10.1042/BJ20110833. pmid:21943356

Morris MK, Saez-Rodriguez J, Clarke DC, Sorger PK, Lauffenburger DA. Training signaling pathway maps to biochemical data with constrained fuzzy logic: quantitative analysis of liver cell responses to inflammatory stimuli. PLoS Comput Biol. 2011;7(3):e1001099. doi: 10.1371/journal.pcbi.1001099. pmid:21408212

Bailey JE. Complex biology with no parameters. Nat Biotechnol. 2001;19(6):503-4. doi: 10.1038/89204. pmid:11385433

Machta BB, Chachra R, Transtrum MK, Sethna JP. Parameter space compression underlies emergent theories and predictive models. Science. 2013;342(6158):604-7. doi: 10.1126/science.1238723. pmid:24179222

Song SO, Chakrabarti A, Varner JD. Ensembles of signal transduction models using Pareto Optimal Ensemble Techniques (POETs). Biotechnol J. 2010;5(7):768-80. doi: 10.1002/biot.201000059. pmid:20665647

Lequieu J, Chakrabarti A, Nayak S, Varner JD. Computational modeling and analysis of insulin induced eukaryotic translation initiation. PLoS Comput Biol. 2011;7(11):e1002263. doi: 10.1371/journal.pcbi.1002263. pmid:22102801

Luan D, Szlam F, Tanaka KA, Barie PS, Varner JD. Ensembles of uncertain mathematical models can identify network response to therapeutic interventions. Mol Biosyst. 2010;6(11):2272-86. doi: 10.1039/b920693k. pmid:20844798

Kalita MK, Sargsyan K, Tian B, Paulucci-Holthauzen A, Najm HN, Debusschere BJ, et al. Sources of cell-to-cell variability in canonical nuclear factor-kB (NF-kB) signaling pathway inferred from single cell dynamic images. J Biol Chem. 2011;286(43):37741-57. doi: 10.1074/jbc.M111.280925. pmid:21868381

Hasenauer J, Waldherr S, Doszczak M, Radde N, Scheurich P, Allgöwer F. Identification of models of heterogeneous cell populations from population snapshot data. BMC Bioinformatics. 2011;12:125. doi: 10.1186/1471-2105-12-125. pmid:21527025

Creixell P, Schoof EM, Erler JT, Linding R. Navigating cancer network attractors for tumor-specific therapy. Nat Biotechnol. 2012;30(9):842-8. doi: 10.1038/nbt.2345. pmid:22965061

Schoeberl B, Eichler-Jonsson C, Gilles ED, Muller G. Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002;20(4):370-375. doi: 10.1038/nbt0402-370. pmid:11923843

Chen WW, Schoeberl B, Jasper PJ, Niepel M, Nielsen UB, Lauffenburger DA, et al. Input-output behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data. Mol Syst Biol. 2009;5:239. doi: 10.1038/msb.2008.74. pmid:19156131

Tasseff R, Nayak S, Song SO, Yen A, Varner JD. Modeling and analysis of retinoic acid induced differentiation of uncommitted precursor cells. Integr Biol (Camb). 2011;3(5):578-591. doi: 10.1039/c0ib00141d

Gadkar KG, Varner J, Doyle FJ 3rd. Model identification of signal transduction networks from data using a state regulator problem. Syst Biol (Stevenage). 2005;2(1):17-30. doi: 10.1049/sb:20045029

Lian Yg, Zhou Qg, Zhang Yj, Zheng Fl. VEGF ameliorates tubulointerstitial fibrosis in unilateral ureteral obstruction mice via inhibition of epithelial-mesenchymal transition. Acta Pharmacol Sin. 2011;32(12):1513-21. doi: 10.1038/aps.2011.111. pmid:21986574

Milo R, Jorgensen P, Moran U, Weber G, Springer M. BioNumbers-the database of key numbers in molecular and cell biology. Nucleic Acids Res. 2010;38(Database issue):D750-3. doi: 10.1093/nar/gkp889. pmid:19854939

Park SY, Lee HE, Li H, Shipitsin M, Gelman R, Polyak K. Heterogeneity for stem cell-related markers according to tumor subtype and histologic stage in breast cancer. Clin Cancer Res. 2010;16(3):876-87. doi: 10.1158/1078-0432.CCR-09-1532. pmid:20103682

Hemavathy K, Ashraf SI, Ip YT. Snail/slug family of repressors: slowly going into the fast lane of development and cancer. Gene. 2000;257(1):1-12. doi: 10.1016/S0378-1119(00)00371-1. pmid:11054563

Hemavathy K, Guru SC, Harris J, Chen JD, Ip YT. Human Slug is a repressor that localizes to sites of active transcription. Mol Cell Biol. 2000;20(14):5087-95. doi: 10.1128/MCB.20.14.5087-5095.2000. pmid:10866665

Dhasarathy A, Phadke D, Mav D, Shah RR, Wade PA. The transcription factors Snail and Slug activate the transforming growth factor-beta signaling pathway in breast cancer. PLoS One. 2011;6(10):e26514. doi: 10.1371/journal.pone.0026514. pmid:22028892

Guaita S, Puig I, Franci C, Garrido M, Dominguez D, Batlle E, et al. Snail induction of epithelial to mesenchymal transition in tumor cells is accompanied by MUC1 repression and ZEB1 expression. J Biol Chem. 2002;277(42):39209-16. doi: 10.1074/jbc.M206400200. pmid:12161443

Jackstadt R, Röh S, Neumann J, Jung P, Hoffmann R, Horst D, et al. AP4 is a mediator of epithelial-mesenchymal transition and metastasis in colorectal cancer. J Exp Med. 2013;210(7):1331-50. doi: 10.1084/jem.20120812. pmid:23752226

Eastman Q, Grosschedl R. Regulation of LEF-1/TCF transcription factors by Wnt and other signals. Curr Opin Cell Biol. 1999;11(2):233-40. doi: 10.1016/S0955-0674(99)80031-3. pmid:10209158

Aramburu J, Yaffe MB, López-Rodríguez C, Cantley LC, Hogan PG, Rao A. Affinity-driven peptide selection of an NFAT inhibitor more selective than cyclosporin A. Science. 1999;285(5436):2129-33. doi: 10.1126/science.285.5436.2129. pmid:10497131

Dhimolea E, Maffini MV, Soto AM, Sonnenschein C. The role of collagen reorganization on mammary epithelial morphogenesis in a 3D culture model. Biomaterials. 2010;31(13):3622-3630. doi: 10.1016/j.biomaterials.2010.01.077. pmid:20149444

Butcher JT, Penrod AM, García AJ, Nerem RM. Unique morphology and focal adhesion development of valvular endothelial cells in static and fluid flow environments. Arterioscler Thromb Vasc Biol. 2004;24(8):1429-34. doi: 10.1161/01.ATV.0000130462.50769.5a. pmid:15117733

Hong JP, Li XM, Li MX, Zheng FL. VEGF suppresses epithelial-mesenchymal transition by inhibiting the expression of Smad3 and miR-192 a Smad3-dependent microRNA. Int J Mol Med. 2013;31(6):1436-42. doi: 10.3892/ijmm.2013.1337. pmid:23588932

Nagy JA, Dvorak AM, Dvorak HF. VEGF-A and the induction of pathological angiogenesis. Annu Rev Pathol. 2007;2:251-75. doi: 10.1146/annurev.pathol.2.010506.134925. pmid:18039100

Ferrara N. VEGF and the quest for tumour angiogenesis factors. Nat Rev Cancer. 2002;2(10):795-803. doi: 10.1038/nrc909. pmid:12360282

Willis BC, Borok Z. TGF-beta-induced EMT: mechanisms and implications for fibrotic lung disease. Am J Physiol Lung Cell Mol Physiol. 2007;293(3):L525-34. doi: 10.1152/ajplung.00163.2007. pmid:17631612

Polyak K, Weinberg RA. Transitions between epithelial and mesenchymal states: acquisition of malignant and stem cell traits. Nat Rev Cancer. 2009;9(4):265-273. doi: 10.1038/nrc2620. pmid:19262571

Strauss R, Li ZY, Liu Y, Beyer I, Persson J, Sova P, et al. Analysis of epithelial and mesenchymal markers in ovarian cancer reveals phenotypic heterogeneity and plasticity. PLoS One. 2011;6(1):e16186. doi: 10.1371/journal.pone.0016186. pmid:21264259

Neve RM, Chin K, Fridlyand J, Yeh J, Baehner FL, Fevr T, et al. A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes. Cancer Cell. 2006;10(6):515-27. doi: 10.1016/j.ccr.2006.10.008. pmid:17157791

Welch-Reardon KM, Wu N, Hughes CCW. A Role for Partial Endothelial-Mesenchymal Transitions in Angiogenesis? Arterioscler Thromb Vasc Biol. 2014; doi: 10.1161/ATVBAHA.114.303220. pmid:25425619

Zajchowski DA, Bartholdi MF, Gong Y, Webster L, Liu HL, Munishkin A, et al. Identification of gene expression profiles that predict the aggressive behavior of breast cancer cells. Cancer Res. 2001;61(13):5168-78. pmid:11431356

Mancini M, Toker A. NFAT proteins: emerging roles in cancer progression. Nat Rev Cancer. 2009;9(11):810-820. doi: 10.1038/nrc2735. pmid:19851316

Singh G, Singh SK, König A, Reutlinger K, Nye MD, Adhikary T, et al. Sequential activation of NFAT and c-Myc transcription factors mediates the TGF-beta switch from a suppressor to a promoter of cancer cell proliferation. J Biol Chem. 2010;285(35):27241-50. doi: 10.1074/jbc.M110.100438. pmid:20516082

Debnath J, Muthuswamy SK, Brugge JS. Morphogenesis and oncogenesis of MCF-10A mammary epithelial acini grown in three-dimensional basement membrane cultures. Methods. 2003;30(3):256-268. doi: 10.1016/S1046-2023(03)00032-X. pmid:12798140

Pearson GW, Hunter T. Real-time imaging reveals that noninvasive mammary epithelial acini can contain motile cells. J Cell Biol. 2007;179(7):1555-67. doi: 10.1083/jcb.200706099. pmid:18166657

O'Brien LE, Tang K, Kats ES, Schutz-Geschwender A, Lipschutz JH, Mostov KE. ERK and MMPs sequentially regulate distinct stages of epithelial tubule development. Dev Cell. 2004;7(1):21-32. doi: 10.1016/j.devcel.2004.06.001. pmid:15239951

Debnath J, Brugge JS. Modelling glandular epithelial cancers in three-dimensional cultures. Nat Rev Cancer. 2005;5(9):675-88. doi: 10.1038/nrc1695. pmid:16148884

Moles CG, Mendes P, Banga JR. Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res. 2003;13(11):2467-74. doi: 10.1101/gr.1262503. pmid:14559783

Rodriguez-Fernandez M, Rehberg M, Kremling A, Banga JR. Simultaneous model discrimination and parameter estimation in dynamic models of cellular systems. BMC Syst Biol. 2013;7:76. doi: 10.1186/1752-0509-7-76. pmid:23938131

Villaverde AF, Banga JR. Reverse engineering and identification in systems biology: strategies, perspectives and challenges. J R Soc Interface. 2014;11(91):20130505. doi: 10.1098/rsif.2013.0505. pmid:24307566

Choi M, Shi J, Jung SH, Chen X, Cho KH. Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to DNA damage. Sci Signal. 2012;5(251):ra83. doi: 10.1126/scisignal.2003363. pmid:23169817

Terfve C, Cokelaer T, Henriques D, MacNamara A, Goncalves E, Morris MK, et al. CellNOptR: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms. BMC Syst Biol. 2012;6:133. doi: 10.1186/1752-0509-6-133. pmid:23079107

Hyduke DR, Palsson BØ. Towards genome-scale signalling network reconstructions. Nat Rev Genet. 2010;11(4):297-307. doi: 10.1038/nrg2750. pmid:20177425

Wayman J, Varner J. Biological systems modeling of metabolic and signaling networks. Curr Opin Chem Eng. 2013;2:365-372. doi: 10.1016/j.coche.2013.09.001

Sainani KL. Meet the Skeptics: Why some doubt biomedical models--and what it takes to win them over. Biomedical Computation Review. 2012; p. 12-18

Lee MW, Vassiliadis VS, Park JM. Individual-based and stochastic modeling of cell population dynamics considering substrate dependency. Biotechnol Bioeng. 2009;103(5):891-9. doi: 10.1002/bit.22327. pmid:19365871

Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci U S A. 2002;99(20):12795-800. doi: 10.1073/pnas.162041399. pmid:12237400

Handl J, Kell DB, Knowles J. Multiobjective optimization in bioinformatics and computational biology. IEEE/ACM Trans Comput Biol Bioinform. 2007;4(2):279-92. doi: 10.1109/TCBB.2007.070203. pmid:17473320

Vega S, Morales AV, Ocaña OH, Valdés F, Fabregat I, Nieto MA. Snail blocks the cell cycle and confers resistance to cell death. Genes Dev. 2004;18(10):1131-43. doi: 10.1101/gad.294104. pmid:15155580

Nawshad A, Hay ED. TGFbeta3 signaling activates transcription of the LEF1 gene to induce epithelial mesenchymal transformation during mouse palate development. J Cell Biol. 2003;163(6):1291-301. doi: 10.1083/jcb.200306024. pmid:14691138

Cano A, Pérez-Moreno MA, Rodrigo I, Locascio A, Blanco MJ, del Barrio MG, et al. The transcription factor snail controls epithelial-mesenchymal transitions by repressing E-cadherin expression. Nat Cell Biol. 2000;2(2):76-83. doi: 10.1038/35000025. pmid:10655586

Arce L, Pate KT, Waterman ML. Groucho binds two conserved regions of LEF-1 for HDAC-dependent repression. BMC Cancer. 2009;9:159. doi: 10.1186/1471-2407-9-159. pmid:19460168

Grumolato L, Liu G, Haremaki T, Mungamuri SK, Mong P, Akiri G, et al. [Beta]-Catenin-Independent Activation of TCF1/LEF1 in Human Hematopoietic Tumor Cells through Interaction with ATF2 Transcription Factors. PLoS Genet. 2013;9(8):e1003603. doi: 10.1371/journal.pgen.1003603. pmid:23966864

Wu Y, Deng J, Rychahou PG, Qiu S, Evers BM, Zhou BP. Stabilization of snail by NF-kappaB is required for inflammation-induced cell migration and invasion. Cancer Cell. 2009;15(5):416-28. doi: 10.1016/j.ccr.2009.03.016. pmid:19411070

Sullivan NJ, Sasser AK, Axel AE, Vesuna F, Raman V, Ramirez N, et al. Interleukin-6 induces an epithelial-mesenchymal transition phenotype in human breast cancer cells. Oncogene. 2009;28(33):2940-7. doi: 10.1038/onc.2009.180. pmid:19581928

Niessen K, Fu Y, Chang L, Hoodless PA, McFadden D, Karsan A. Slug is a direct Notch target required for initiation of cardiac cushion cellularization. J Cell Biol. 2008;182(2):315-25. doi: 10.1083/jcb.200710067. pmid:18663143

Liu Y, El-Naggar S, Darling DS, Higashi Y, Dean DC. Zeb1 links epithelial-mesenchymal transition and cellular senescence. Development. 2008;135(3):579-88. doi: 10.1242/dev.007047. pmid:18192284

Hong T, Watanabe K, Ta CH, Villarreal-Ponce A, Nie Q, Dai X. An Ovol2-Zeb1 Mutual Inhibitory Circuit Governs Bidirectional and Multi-step Transition between Epithelial and Mesenchymal States. PLoS Comput Biol. 2015;11(11):e1004569. doi: 10.1371/journal.pcbi.1004569. pmid:26554584

Bullock MD, Sayan AE, Packham GK, Mirnezami AH. MicroRNAs: critical regulators of epithelial to mesenchymal (EMT) and mesenchymal to epithelial transition (MET) in cancer progression. Biol Cell. 2012;104(1):3-12. doi: 10.1111/boc.201100115. pmid:22188537

Lu M, Jolly MK, Levine H, Onuchic JN, Ben-Jacob E. MicroRNA-based regulation of epithelial-hybrid-mesenchymal fate determination. Proc Natl Acad Sci U S A. 2013;110(45):18144-9. doi: 10.1073/pnas.1318192110. pmid:24154725

Zhang J, Tian XJ, Zhang H, Teng Y, Li R, Bai F, et al. TGF-beta-induced epithelial-to-mesenchymal transition proceeds through stepwise activation of multiple feedback loops. Sci Signal. 2014;7(345):ra91. doi: 10.1126/scisignal.2005304. pmid:25270257

Tian XJ, Zhang H, Xing J. Coupled reversible and irreversible bistable switches underlying TGFb-induced epithelial to mesenchymal transition. Biophys J. 2013;105(4):1079-89. doi: 10.1016/j.bpj.2013.07.011. pmid:23972859

Linding R, Jensen LJ, Ostheimer GJ, van Vugt MATM, Jørgensen C, Miron IM, et al. Systematic discovery of in vivo phosphorylation networks. Cell. 2007;129(7):1415-26. doi: 10.1016/j.cell.2007.05.052. pmid:17570479

Jensen LJ, Kuhn M, Stark M, Chaffron S, Creevey C, Muller J, et al. STRING 8-a global view on proteins and their functional interactions in 630 organisms. Nucleic Acids Res. 2009;37(Database issue):D412-6. doi: 10.1093/nar/gkn760

Zhou BP, Deng J, Xia W, Xu J, Li YM, Gunduz M, et al. Dual regulation of Snail by GSK-3beta-mediated phosphorylation in control of epithelial-mesenchymal transition. Nat Cell Biol. 2004;6(10):931-40. doi: 10.1038/ncb1173. pmid:15448698

Larue L, Bellacosa A. Epithelial-mesenchymal transition in development and cancer: role of phosphatidylinositol 3' kinase/AKT pathways. Oncogene. 2005;24(50):7443-54. doi: 10.1038/sj.onc.1209091. pmid:16288291

Derynck R, Zhang YE. Smad-dependent and Smad-independent pathways in TGF-beta family signalling. Nature. 2003;425(6958):577-84. doi: 10.1038/nature02006. pmid:14534577

Massagué J, Seoane J, Wotton D. Smad transcription factors. Genes Dev. 2005;19(23):2783-810. doi: 10.1101/gad.1350705. pmid:16322555

Nawshad A, Medici D, Liu CC, Hay ED. TGFbeta3 inhibits E-cadherin gene expression in palate medial-edge epithelial cells through a Smad2-Smad4-LEF1 transcription complex. J Cell Sci. 2007;120(Pt 9):1646-53. doi: 10.1242/jcs.003129. pmid:17452626

Medici D, Potenta S, Kalluri R. Transforming growth factor-[Beta]2 promotes Snail-mediated endothelial-mesenchymal transition through convergence of Smad-dependent and Smad-independent signalling. Biochem J. 2011;437(3):515-520. doi: 10.1042/BJ20101500. pmid:21585337

Suehiro Ji, Kanki Y, Makihara C, Schadler K, Miura M, Manabe Y, et al. Genome-wide approaches reveal functional vascular endothelial growth factor (VEGF)-inducible nuclear factor of activated T cells (NFAT) c1 binding to angiogenesis-related genes in the endothelium. J Biol Chem. 2014;289(42):29044-59. doi: 10.1074/jbc.M114.555235. pmid:25157100

Hu CY, Varner JD, Lucks JB. Generating Effective Models and Parameters for RNA Genetic Circuits. ACS Synth Biol. 2015;4(8):914-26. doi: 10.1021/acssynbio.5b00077. pmid:26046393

Wayman JA, Sagar A, Varner JD. Dynamic Modeling of Cell-Free Biochemical Networks Using Effective Kinetic Models. Processes. 2015;3(1):138. doi: 10.3390/pr3010138

Sagar A, Varner JD. Dynamic Modeling of the Human Coagulation Cascade Using Reduced Order Effective Kinetic Models. Processes. 2015;3(1):178. doi: 10.3390/pr3010178

Bezanson J, Edelman A, Karpinski S, Shah VB. Julia: A Fresh Approach to Numerical Computing. CoRR. 2014;abs/1411.1607

Hindmarsh A, Brown P, Grant K, Lee S, Serban R, Shumaker D, et al. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software. 2005;31:363-396. doi: 10.1145/1089014.1089020

Bassen D, Vilkhovoy M, Minot M, Butcher JT, Varner JD. JuPOETs: A Constrained Multiobjective Optimization Approach to Estimate Biochemical Model Ensembles in the Julia Programming Language. bioRxiv. 2016; doi: 10.1101/056044

Kohavi R. A study of cross-validation and bootstrap for accuracy estimation and model selection. In: International joint Conference on artificial intelligence. vol. 14. Citeseer; 1995. p. 1137-1145

Chang H, Park C, Parvin B. Quantitative Representation of Three-dimensional Cell Culture Models. In: Proceedings of the 2007 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Washington, DC, USA, April 12-16, 2007; 2007. p. 89-92. Available from: http://dx.doi.org/10.1109/ISBI.2007.356795

Polizzotti L, Basak O, Bjornsson C, Shubert K, Yener B, Plopper G. Novel Image Analysis Approach Quantifies Morphological Characteristics of 3D Breast Culture Acini with Varying Metastatic Potentials. J Biomed Biotech. 2012;2012:1-16. doi: 10.1155/2012/102036