Understanding large scale sequencing datasets through changes to protein folding

The expansion of high-quality, low-cost sequencing has created an enormous opportunity to understand how genetic variants alter cellular behaviour in disease. The high diversity of mutations observed has however drawn a spotlight onto the need for predictive modelling of mutational effects on phenotype from variants of uncertain significance. This is particularly important in the clinic due to the potential value in guiding clinical diagnosis and patient treatment. Recent computational modelling has highlighted the importance of mutation induced protein misfolding as a common mechanism for loss of protein or domain function, aided by developments in methods that make large computational screens tractable. Here we review recent applications of this approach to different genes, and how they have enabled and supported subsequent studies. We further discuss developments in the approach and the role for the approach in light of increasingly high throughput experimental approaches.

Interpreting the effect of mutations to protein binding sites from large-scale genomic screens

Predicting the functionality of missense mutations is extremely difficult. Large-scale genomic screens are commonly performed to identify mutational correlates or drivers of disease and treatment resistance, but interpretation of how these mutations impact protein function is limited. One such consequence of mutations to a protein is to impact its ability to bind and interact with partners or small molecules such as ATP, thereby modulating its function. Multiple methods exist for predicting the impact of a single mutation on protein-protein binding energy, but it is difficult in the context of a genomic screen to understand if these mutations with large impacts on binding are more common than statistically expected. We present a methodology for taking mutational data from large-scale genomic screens and generating functional and statistical insights into their role in the binding of proteins both with each other and their small molecule ligands. This allows a quantitative and statistical analysis to determine whether mutations impacting protein binding or ligand interactions are occurring more or less frequently than expected by chance. We achieve this by calculating the potential impact of any possible mutation and comparing an expected distribution to the observed mutations. This method is applied to examples demonstrating its ability to interpret mutations involved in protein-protein binding, protein-DNA interactions, and the evolution of therapeutic resistance.

Large-scale computational modelling of the M1 and M2 synovial macrophages in rheumatoid arthritis

Macrophages play an essential role in rheumatoid arthritis. Depending on their phenotype (M1 or M2), they can play a role in the initiation or resolution of inflammation. The M1/M2 ratio in rheumatoid arthritis is higher than in healthy controls. Despite this, no treatment targeting specifically macrophages is currently used in clinics. Thus, devising strategies to selectively deplete proinflammatory macrophages and promote anti-inflammatory macrophages could be a promising therapeutic approach. State-of-the-art molecular interaction maps of M1 and M2 macrophages in rheumatoid arthritis are available and represent a dense source of knowledge; however, these maps remain limited by their static nature. Discrete dynamic modelling can be employed to study the emergent behaviours of these systems. Nevertheless, handling such large-scale models is challenging. Due to their massive size, it is computationally demanding to identify biologically relevant states in a cell- and disease-specific context. In this work, we developed an efficient computational framework that converts molecular interaction maps into Boolean models using the CaSQ tool. Next, we used a newly developed version of the BMA tool deployed to a high-performance computing cluster to identify the models' steady states. The identified attractors are then validated using gene expression data sets and prior knowledge. We successfully applied our framework to generate and calibrate the M1 and M2 macrophage Boolean models for rheumatoid arthritis. Using KO simulations, we identified NFkB, JAK1/JAK2, and ERK1/Notch1 as potential targets that could selectively suppress proinflammatory macrophages and GSK3B as a promising target that could promote anti-inflammatory macrophages in rheumatoid arthritis.

SMAD4 and KCNQ3 alterations are associated with lymph node metastases in oesophageal adenocarcinoma

Metastasis in oesophageal adenocarcinoma (OAC) is an important predictor of survival. Radiological staging is used to stage metastases in patients, and guide treatment selection, but is limited by the accuracy of the approach. Improvements in staging will lead to improved clinical decision making and patient outcomes. Sequencing studies on primary tumours and pre-cancerous tissue have revealed the mutational landscape of OAC, and increasingly cheap and widespread sequencing approaches offer the potential to improve staging assessment. In this work we present an analysis of lymph node metastases found by radiological and pathological sampling, identifying new roles of the genes SMAD4 and KCNQ3 in metastasis. Through transcriptomic analysis we find that both genes are associated with canonical Wnt pathway activity, but KCNQ3 is uniquely associated with changes in planar cell polaritiy associated with non-canonical Wnt signalling. We go on to validate our observations in KCNQ3 in cell line and xenograph systems, showing that overexpression of KCNQ3 reduces wound closure and the number of metastases observed. Our results suggest both genes as novel biomarkers of metastatic risk and offer new potential routes to drug targeting.

KCNQ potassium channels modulate Wnt activity in gastro-oesophageal adenocarcinomas

Voltage-sensitive potassium channels play an important role in controlling membrane potential and ionic homeostasis in the gut and have been implicated in gastrointestinal (GI) cancers. Through large-scale analysis of 897 patients with gastro-oesophageal adenocarcinomas (GOAs) coupled with in vitro models, we find KCNQ family genes are mutated in ∼30% of patients, and play therapeutically targetable roles in GOA cancer growth. KCNQ1 and KCNQ3 mediate the WNT pathway and MYC to increase proliferation through resultant effects on cadherin junctions. This also highlights novel roles of KCNQ3 in non-excitable tissues. We also discover that activity of KCNQ3 sensitises cancer cells to existing potassium channel inhibitors and that inhibition of KCNQ activity reduces proliferation of GOA cancer cells. These findings reveal a novel and exploitable role of potassium channels in the advancement of human cancer, and highlight that supplemental treatments for GOAs may exist through KCNQ inhibitors.

Somatic mutations in facial skin from countries of contrasting skin cancer risk

The incidence of keratinocyte cancer (basal cell and squamous cell carcinomas of the skin) is 17-fold lower in Singapore than the UK1-3, despite Singapore receiving 2-3 times more ultraviolet (UV) radiation4,5. Aging skin contains somatic mutant clones from which such cancers develop6,7. We hypothesized that differences in keratinocyte cancer incidence may be reflected in the normal skin mutational landscape. Here we show that, compared to Singapore, aging facial skin from populations in the UK has a fourfold greater mutational burden, a predominant UV mutational signature, increased copy number aberrations and increased mutant TP53 selection. These features are shared by keratinocyte cancers from high-incidence and low-incidence populations8-13. In Singaporean skin, most mutations result from cell-intrinsic processes; mutant NOTCH1 and NOTCH2 are more strongly selected than in the UK. Aging skin in a high-incidence country has multiple features convergent with cancer that are not found in a low-risk country. These differences may reflect germline variation in UV-protective genes.

Mutations observed in somatic evolution reveal underlying gene mechanisms

Highly sensitive DNA sequencing techniques have allowed the discovery of large numbers of somatic mutations in normal tissues. Some mutations confer a competitive advantage over wild-type cells, generating expanding clones that spread through the tissue. Competition between mutant clones leads to selection. This process can be considered a large scale, in vivo screen for mutations increasing cell fitness. It follows that somatic missense mutations may offer new insights into the relationship between protein structure, function and cell fitness. We present a flexible statistical method for exploring the selection of structural features in data sets of somatic mutants. We show how this approach can evidence selection of specific structural features in key drivers in aged tissues. Finally, we show how drivers may be classified as fitness-enhancing and fitness-suppressing through different patterns of mutation enrichment. This method offers a route to understanding the mechanism of protein function through in vivo mutant selection.

Notch1 mutations drive clonal expansion in normal esophageal epithelium but impair tumor growth

NOTCH1 mutant clones occupy the majority of normal human esophagus by middle age but are comparatively rare in esophageal cancers, suggesting NOTCH1 mutations drive clonal expansion but impede carcinogenesis. Here we test this hypothesis. Sequencing NOTCH1 mutant clones in aging human esophagus reveals frequent biallelic mutations that block NOTCH1 signaling. In mouse esophagus, heterozygous Notch1 mutation confers a competitive advantage over wild-type cells, an effect enhanced by loss of the second allele. Widespread Notch1 loss alters transcription but has minimal effects on the epithelial structure and cell dynamics. In a carcinogenesis model, Notch1 mutations were less prevalent in tumors than normal epithelium. Deletion of Notch1 reduced tumor growth, an effect recapitulated by anti-NOTCH1 antibody treatment. Notch1 null tumors showed reduced proliferation. We conclude that Notch1 mutations in normal epithelium are beneficial as wild-type Notch1 favors tumor expansion. NOTCH1 blockade may have therapeutic potential in preventing esophageal squamous cancer.

The NALCN channel regulates metastasis and nonmalignant cell dissemination

We identify the sodium leak channel non-selective protein (NALCN) as a key regulator of cancer metastasis and nonmalignant cell dissemination. Among 10,022 human cancers, NALCN loss-of-function mutations were enriched in gastric and colorectal cancers. Deletion of Nalcn from gastric, intestinal or pancreatic adenocarcinomas in mice did not alter tumor incidence, but markedly increased the number of circulating tumor cells (CTCs) and metastases. Treatment of these mice with gadolinium-a NALCN channel blocker-similarly increased CTCs and metastases. Deletion of Nalcn from mice that lacked oncogenic mutations and never developed cancer caused shedding of epithelial cells into the blood at levels equivalent to those seen in tumor-bearing animals. These cells trafficked to distant organs to form normal structures including lung epithelium, and kidney glomeruli and tubules. Thus, NALCN regulates cell shedding from solid tissues independent of cancer, divorcing this process from tumorigenesis and unmasking a potential new target for antimetastatic therapies.

Heterogeneity of the cancer cell line metabolic landscape

The unravelling of the complexity of cellular metabolism is in its infancy. Cancer-associated genetic alterations may result in changes to cellular metabolism that aid in understanding phenotypic changes, reveal detectable metabolic signatures, or elucidate vulnerabilities to particular drugs. To understand cancer-associated metabolic transformation, we performed untargeted metabolite analysis of 173 different cancer cell lines from 11 different tissues under constant conditions for 1,099 different species using mass spectrometry (MS). We correlate known cancer-associated mutations and gene expression programs with metabolic signatures, generating novel associations of known metabolic pathways with known cancer drivers. We show that metabolic activity correlates with drug sensitivity and use metabolic activity to predict drug response and synergy. Finally, we study the metabolic heterogeneity of cancer mutations across tissues, and find that genes exhibit a range of context specific, and more general metabolic control.

p53 mutation in normal esophagus promotes multiple stages of carcinogenesis but is constrained by clonal competition

Aging normal human oesophagus accumulates TP53 mutant clones. These are the origin of most oesophageal squamous carcinomas, in which biallelic TP53 disruption is almost universal. However, how p53 mutant clones expand and contribute to cancer development is unclear. Here we show that inducing the p53R245W mutant in single oesophageal progenitor cells in transgenic mice confers a proliferative advantage and clonal expansion but does not disrupt normal epithelial structure. Loss of the remaining p53 allele in mutant cells results in genomically unstable p53R245W/null epithelium with giant polyaneuploid cells and copy number altered clones. In carcinogenesis, p53 mutation does not initiate tumour formation, but tumours developing from areas with p53 mutation and LOH are larger and show extensive chromosomal instability compared to lesions arising in wild type epithelium. We conclude that p53 has distinct functions at different stages of carcinogenesis and that LOH within p53 mutant clones in normal epithelium is a critical step in malignant transformation.

HOXA9 has the hallmarks of a biological switch with implications in blood cancers

Blood malignancies arise from the dysregulation of haematopoiesis. The type of blood cell and the specific order of oncogenic events initiating abnormal growth ultimately determine the cancer subtype and subsequent clinical outcome. HOXA9 plays an important role in acute myeloid leukaemia (AML) prognosis by promoting blood cell expansion and altering differentiation; however, the function of HOXA9 in other blood malignancies is still unclear. Here, we highlight the biological switch and prognosis marker properties of HOXA9 in AML and chronic myeloproliferative neoplasms (MPN). First, we establish the ability of HOXA9 to stratify AML patients with distinct cellular and clinical outcomes. Then, through the use of a computational network model of MPN, we show that the self-activation of HOXA9 and its relationship to JAK2 and TET2 can explain the branching progression of JAK2/TET2 mutant MPN patients towards divergent clinical characteristics. Finally, we predict a connection between the RUNX1 and MYB genes and a suppressive role for the NOTCH pathway in MPN diseases.

Tumor-Derived Lactic Acid Modulates Activation and Metabolic Status of Draining Lymph Node Stroma

Communication between tumors and the stroma of tumor-draining lymph nodes (TDLN) exists before metastasis arises, altering the structure and function of the TDLN niche. Transcriptional profiling of fibroblastic reticular cells (FRC), the dominant stromal population of lymph nodes, has revealed that FRCs in TDLNs are reprogrammed. However, the tumor-derived factors driving the changes in FRCs remain to be identified. Taking an unbiased approach, we have shown herein that lactic acid (LA), a metabolite released by cancer cells, was not only secreted by B16.F10 and 4T1 tumors in high amounts, but also that it was enriched in TDLNs. LA supported an upregulation of Podoplanin (Pdpn) and Thy1 and downregulation of IL7 in FRCs of TDLNs, making them akin to activated fibroblasts found at the primary tumor site. Furthermore, we found that tumor-derived LA altered mitochondrial function of FRCs in TDLNs. Thus, our results demonstrate a mechanism by which a tumor-derived metabolite connected with a low pH environment modulates the function of fibroblasts in TDLNs. How lymph node function is perturbed to support cancer metastases remains unclear. The authors show that tumor-derived LA drains to lymph nodes where it modulates the function of lymph node stromal cells, prior to metastatic colonization.

Mutant clones in normal epithelium outcompete and eliminate emerging tumours

Human epithelial tissues accumulate cancer-driver mutations with age1-9, yet tumour formation remains rare. The positive selection of these mutations suggests that they alter the behaviour and fitness of proliferating cells10-12. Thus, normal adult tissues become a patchwork of mutant clones competing for space and survival, with the fittest clones expanding by eliminating their less competitive neighbours11-14. However, little is known about how such dynamic competition in normal epithelia influences early tumorigenesis. Here we show that the majority of newly formed oesophageal tumours are eliminated through competition with mutant clones in the adjacent normal epithelium. We followed the fate of nascent, microscopic, pre-malignant tumours in a mouse model of oesophageal carcinogenesis and found that most were rapidly lost with no indication of tumour cell death, decreased proliferation or an anti-tumour immune response. However, deep sequencing of ten-day-old and one-year-old tumours showed evidence of selection on the surviving neoplasms. Induction of highly competitive clones in transgenic mice increased early tumour removal, whereas pharmacological inhibition of clonal competition reduced tumour loss. These results support a model in which survival of early neoplasms depends on their competitive fitness relative to that of mutant clones in the surrounding normal tissue. Mutant clones in normal epithelium have an unexpected anti-tumorigenic role in purging early tumours through cell competition, thereby preserving tissue integrity.

Simulations reveal that different responses to cell crowding determine the expansion of p53 and Notch mutant clones in squamous epithelia

During ageing, normal epithelial tissues progressively accumulate clones carrying mutations that increase mutant cell fitness above that of wild-type cells. Such mutants spread widely through the tissues, yet despite this cellular homeostasis and functional integrity of the epithelia are maintained. Two of the genes most commonly mutated in human skin and oesophagus are p53 and Notch1, both of which are also recurrently mutated in cancers of these tissues. From observations taken in human and mouse epithelia, we find that clones carrying p53 and Notch pathway mutations have different clone dynamics which can be explained by their different responses to local cell crowding. p53 mutant clone growth in mouse epidermis approximates a logistic curve, but feedbacks responding to local crowding are required to maintain tissue homeostasis. We go on to show that the observed ability of Notch pathway mutant cells to displace the wild-type population in the mouse oesophageal epithelium reflects a local density feedback that affects both mutant and wild-type cells equally. We then show how these distinct feedbacks are consistent with the distribution of mutations observed in human datasets and are suggestive of a putative mechanism to constrain these cancer-associated mutants.

SARS-CoV-2 Variants Are Selecting for Spike Protein Mutations That Increase Protein Stability

The emergence of variants of SARS-CoV-2 with mutations in their spike protein are a major cause for concern for the efficacy of vaccines and control of the pandemic. We show that mutations in the spike protein of SARS-CoV-2 are selecting for amino acid changes that result in a more thermodynamically stable protein than expected from background. We suggest that the computationally efficient analysis of mutational stability may aid in early screening of variants.

Cancer-causing BRCA2 missense mutations disrupt an intracellular protein assembly mechanism to disable genome maintenance

Cancer-causing missense mutations in the 3418 amino acid BRCA2 breast and ovarian cancer suppressor protein frequently affect a short (∼340 residue) segment in its carboxyl-terminal domain (DBD). Here, we identify a shared molecular mechanism underlying their pathogenicity. Pathogenic BRCA2 missense mutations cluster in the DBD’s helical domain (HD) and OB1-fold motifs, which engage the partner protein DSS1. Pathogenic - but not benign – DBD mutations weaken or abolish DSS1-BRCA2 assembly, provoking mutant BRCA2 oligomers that are excluded from the cell nucleus, and disable DNA repair by homologous DNA recombination (HDR). DSS1 inhibits the intracellular oligomerization of wildtype, but not mutant, forms of BRCA2. Remarkably, DSS1 expression corrects defective HDR in cells bearing pathogenic BRCA2 missense mutants with weakened, but not absent, DSS1 binding. Our findings identify a DSS1-mediated intracellular protein assembly mechanism that is disrupted by cancer-causing BRCA2 missense mutations, and suggest an approach for its therapeutic correction.

Methods for analysing lineage tracing datasets

A single population of progenitor cells maintains many epithelial tissues. Transgenic mouse cell tracking has frequently been used to study the growth dynamics of competing clones in these tissues. A mathematical model (the 'single-progenitor model') has been argued to reproduce the observed progenitor dynamics accurately. This requires three parameters to describe the growth dynamics observed in transgenic mouse cell tracking-a division rate, a stratification rate and the probability of dividing symmetrically. Deriving these parameters is a time intensive and complex process. We compare the alternative strategies for analysing this source of experimental data, identifying an approximate Bayesian computation-based approach as the best in terms of efficiency and appropriate error estimation. We support our findings by explicitly modelling biological variation and consider the impact of different sampling regimes. All tested solutions are made available to allow new datasets to be analysed following our workflows. Based on our findings, we make recommendations for future experimental design.

Selection of Oncogenic Mutant Clones in Normal Human Skin Varies with Body Site

Skin cancer risk varies substantially across the body, yet how this relates to the mutations found in normal skin is unknown. Here we mapped mutant clones in skin from high- and low-risk sites. The density of mutations varied by location. The prevalence of NOTCH1 and FAT1 mutations in forearm, trunk, and leg skin was similar to that in keratinocyte cancers. Most mutations were caused by ultraviolet light, but mutational signature analysis suggested differences in DNA-repair processes between sites. Eleven mutant genes were under positive selection, with TP53 preferentially selected in the head and FAT1 in the leg. Fine-scale mapping revealed 10% of clones had copy-number alterations. Analysis of hair follicles showed mutations in the upper follicle resembled adjacent skin, but the lower follicle was sparsely mutated. Normal skin is a dense patchwork of mutant clones arising from competitive selection that varies by location. SIGNIFICANCE: Mapping mutant clones across the body reveals normal skin is a dense patchwork of mutant cells. The variation in cancer risk between sites substantially exceeds that in mutant clone density. More generally, mutant genes cannot be assigned as cancer drivers until their prevalence in normal tissue is known.See related commentary by De Dominici and DeGregori, p. 227.This article is highlighted in the In This Issue feature, p. 211.

Spatial competition shapes the dynamic mutational landscape of normal esophageal epithelium

During aging, progenitor cells acquire mutations, which may generate clones that colonize the surrounding tissue. By middle age, normal human tissues, including the esophageal epithelium (EE), become a patchwork of mutant clones. Despite their relevance for understanding aging and cancer, the processes that underpin mutational selection in normal tissues remain poorly understood. Here, we investigated this issue in the esophageal epithelium of mutagen-treated mice. Deep sequencing identified numerous mutant clones with multiple genes under positive selection, including Notch1, Notch2 and Trp53, which are also selected in human esophageal epithelium. Transgenic lineage tracing revealed strong clonal competition that evolved over time. Clone dynamics were consistent with a simple model in which the proliferative advantage conferred by positively selected mutations depends on the nature of the neighboring cells. When clones with similar competitive fitness collide, mutant cell fate reverts towards homeostasis, a constraint that explains how selection operates in normal-appearing epithelium.

A single-progenitor model as the unifying paradigm of epidermal and esophageal epithelial maintenance in mice

In adult skin epidermis and the epithelium lining the esophagus cells are constantly shed from the tissue surface and replaced by cell division. Tracking genetically labelled cells in transgenic mice has given insight into cell behavior, but conflicting models appear consistent with the results. Here, we use an additional transgenic assay to follow cell division in mouse esophagus and the epidermis at multiple body sites. We find that proliferating cells divide at a similar rate, and place bounds on the distribution cell cycle times. By including these results in a common analytic approach, we show that data from eight lineage tracing experiments is consistent with tissue maintenance by a single population of proliferating cells. The outcome of a given cell division is unpredictable but, on average, the likelihood of producing proliferating and differentiating cells is equal, ensuring cellular homeostasis. These findings are key to understanding squamous epithelial homeostasis and carcinogenesis.

Constructing and Analyzing Computational Models of Cell Signaling with BioModelAnalyzer

BioModelAnalyzer (BMA) is an open-source graphical tool for the development of executable models of protein and gene networks within cells. Based upon the Qualitative Networks formalism, the user can rapidly construct large networks, either manually or by connecting motifs selected from a built-in library. After the appropriate functions for each variable are defined, the user has access to three analysis engines to test the model. In addition to standard simulation tools, BMA includes an interface to the stability-testing algorithm and to a graphical Linear Temporal Logic (LTL) editor and analysis tool. Alongside this, we have developed a novel ChatBot to aid users constructing LTL queries and to explain the interface and run through tutorials. Here we present worked examples of model construction and testing via the interface. As an initial example, we discuss fate decisions in Dictyostelium discoidum and cAMP signaling. We go on to describe the workflow leading to the construction of a published model of the germline of C. elegans. Finally, we demonstrate how to construct simple models from the built-in network motif library. © 2020 by John Wiley & Sons, Inc. Basic Protocol 1: Modeling the signaling network of Dictyostelium discoidum Basic Protocol 2: Modeling the germline progression of Caenorhabditis elegans Basic Protocol 3: Constructing a model of the cell cycle using motifs.

Cellular survival over genomic perfection

DNA repair pathways permit some damage, leading to mutagenesis but not always cancer

Relating evolutionary selection and mutant clonal dynamics in normal epithelia

Cancer develops from mutated cells in normal tissues. Whether somatic mutations alter normal cell dynamics is key to understanding cancer risk and guiding interventions to reduce it. An analysis of the first incomplete moment of size distributions of clones carrying cancer-associated mutations in normal human eyelid skin gives a good fit with neutral drift, arguing mutations do not affect cell fate. However, this suggestion conflicts with genetic evidence in the same dataset that argues for strong positive selection of a subset of mutations. This implies cells carrying these mutations have a competitive advantage over normal cells, leading to large clonal expansions within the tissue. In the normal epithelium, clone growth is constrained by the limited size of the proliferating compartment and competition with surrounding cells. We show that if these factors are taken into account, the first incomplete moment of the clone size distribution is unable to exclude non-neutral behaviour. Furthermore, experimental factors can make a non-neutral clone size distribution appear neutral. We validate these principles with a new experimental dataset showing that when experiments are appropriately designed, the first incomplete moment can be a useful indicator of non-neutral competition. Finally, we discuss the complex relationship between mutant clone sizes and genetic selection.

Epidermal Tissue Adapts to Restrain Progenitors Carrying Clonal p53 Mutations

Aging human tissues, such as sun-exposed epidermis, accumulate a high burden of progenitor cells that carry oncogenic mutations. However, most progenitors carrying such mutations colonize and persist in normal tissue without forming tumors. Here, we investigated tissue-level constraints on clonal progenitor behavior by inducing a single-allele p53 mutation (Trp53R245W; p53∗/wt), prevalent in normal human epidermis and squamous cell carcinoma, in transgenic mouse epidermis. p53∗/wt progenitors initially outcompeted wild-type cells due to enhanced proliferation, but subsequently reverted toward normal dynamics and homeostasis. Physiological doses of UV light accelerated short-term expansion of p53∗/wt clones, but their frequency decreased with protracted irradiation, possibly due to displacement by UV-induced mutant clones with higher competitive fitness. These results suggest multiple mechanisms restrain the proliferation of p53∗/wt progenitors, thereby maintaining epidermal integrity.

Exploring the role of stromal osmoregulation in cancer and disease using executable modelling

Osmotic regulation is a vital homoeostatic process in all cells and tissues. Cells initially respond to osmotic stresses by activating transmembrane transport proteins to move osmotically active ions. Disruption of ion and water transport is frequently observed in cellular transformations such as cancer. We report that genes involved in membrane transport are significantly deregulated in many cancers, and that their expression can distinguish cancer cells from normal cells with a high degree of accuracy. We present an executable model of osmotic regulation and membrane transport in mammalian cells, providing a mechanistic explanation for phenotype change in varied disease states, and accurately predicting behaviour from single cell expression data. We also predict key proteins involved in cellular transformation, SLC4A3 (AE3), and SLC9A1 (NHE1). Furthermore, we predict and verify a synergistic drug combination in vitro, of sodium and chloride channel inhibitors, which target the osmoregulatory network to reduce cancer-associated phenotypes in fibroblasts.

A toolbox for discrete modelling of cell signalling dynamics

We present a library of network motifs for the development of complex and realistic biological network models using the BioModelAnalyzer, and demonstrate their wider value by using them to construct a model of the cell cycle.

Tumor-induced stromal reprogramming drives lymph node transformation

Lymph node (LN) stromal cells, particularly fibroblastic reticular cells (FRCs), provide critical structural support and regulate immunity, tolerance and the transport properties of LNs. For many tumors, metastasis to the LNs is predictive of poor prognosis. However, the stromal contribution to the evolving microenvironment of tumor-draining LNs (TDLNs) remains poorly understood. Here we found that FRCs specifically of TDLNs proliferated in response to tumor-derived cues and that the network they formed was remodeled. Comparative transcriptional analysis of FRCs from non-draining LNs and TDLNs demonstrated reprogramming of key pathways, including matrix remodeling, chemokine and/or cytokine signaling, and immunological functions such as the recruitment, migration and activation of leukocytes. In particular, downregulation of the expression of FRC-derived chemokine CCL21 and cytokine IL-7 were accompanied by altered composition and aberrant localization of immune-cell populations. Our data indicate that following exposure to tumor-derived factors, the stroma of TDLNs adapts on multiple levels to exhibit features typically associated with immunosuppression.

Emergent stem cell homeostasis in the C. elegans germline is revealed by hybrid modeling

The establishment of homeostasis among cell growth, differentiation, and apoptosis is of key importance for organogenesis. Stem cells respond to temporally and spatially regulated signals by switching from mitotic proliferation to asymmetric cell division and differentiation. Executable computer models of signaling pathways can accurately reproduce a wide range of biological phenomena by reducing detailed chemical kinetics to a discrete, finite form. Moreover, coordinated cell movements and physical cell-cell interactions are required for the formation of three-dimensional structures that are the building blocks of organs. To capture all these aspects, we have developed a hybrid executable/physical model describing stem cell proliferation, differentiation, and homeostasis in the Caenorhabditis elegans germline. Using this hybrid model, we are able to track cell lineages and dynamic cell movements during germ cell differentiation. We further show how apoptosis regulates germ cell homeostasis in the gonad, and propose a role for intercellular pressure in developmental control. Finally, we use the model to demonstrate how an executable model can be developed from the hybrid system, identifying a mechanism that ensures invariance in fate patterns in the presence of instability.

Sidekick for Membrane Simulations: Automated Ensemble Molecular Dynamics Simulations of Transmembrane Helices

The interactions of transmembrane (TM) α-helices with the phospholipid membrane and with one another are central to understanding the structure and stability of integral membrane proteins. These interactions may be analyzed via coarse grained molecular dynamics (CGMD) simulations. To obtain statistically meaningful analysis of TM helix interactions, large (N ca. 100) ensembles of CGMD simulations are needed. To facilitate the running and analysis of such ensembles of simulations, we have developed Sidekick, an automated pipeline software for performing high throughput CGMD simulations of α-helical peptides in lipid bilayer membranes. Through an end-to-end approach, which takes as input a helix sequence and outputs analytical metrics derived from CGMD simulations, we are able to predict the orientation and likelihood of insertion into a lipid bilayer of a given helix of a family of helix sequences. We illustrate this software via analyses of insertion into a membrane of short hydrophobic TM helices containing a single cationic arginine residue positioned at different positions along the length of the helix. From analyses of these ensembles of simulations, we estimate apparent energy barriers to insertion which are comparable to experimentally determined values. In a second application, we use CGMD simulations to examine the self-assembly of dimers of TM helices from the ErbB1 receptor tyrosine kinase and analyze the numbers of simulation repeats necessary to obtain convergence of simple descriptors of the mode of packing of the two helices within a dimer. Our approach offers a proof-of-principle platform for the further employment of automation in large ensemble CGMD simulations of membrane proteins.

Probing the Solution Structure of IκB Kinase (IKK) Subunit γ and Its Interaction with Kaposi Sarcoma-associated Herpes Virus Flice-interacting Protein and IKK Subunit β by EPR Spectroscopy

Viral flice-interacting protein (vFLIP), encoded by the oncogenic Kaposi sarcoma-associated herpes virus (KSHV), constitutively activates the canonical nuclear factor κ-light-chain-enhancer of activated B cells (NF-κB) pathway. This is achieved through subversion of the IκB kinase (IKK) complex (or signalosome), which involves a physical interaction between vFLIP and the modulatory subunit IKKγ. Although this interaction has been examined both in vivo and in vitro, the mechanism by which vFLIP activates the kinase remains to be determined. Because IKKγ functions as a scaffold, recruiting both vFLIP and the IKKα/β subunits, it has been proposed that binding of vFLIP could trigger a structural rearrangement in IKKγ conducive to activation. To investigate this hypothesis we engineered a series of mutants along the length of the IKKγ molecule that could be individually modified with nitroxide spin labels. Subsequent distance measurements using electron paramagnetic resonance spectroscopy combined with molecular modeling and molecular dynamics simulations revealed that IKKγ is a parallel coiled-coil whose response to binding of vFLIP or IKKβ is localized twisting/stiffening and not large-scale rearrangements. The coiled-coil comprises N- and C-terminal regions with distinct registers accommodated by a twist: this structural motif is exploited by vFLIP, allowing it to bind and subsequently activate the NF-κB pathway. In vivo assays confirm that NF-κB activation by vFLIP only requires the N-terminal region up to the transition between the registers, which is located directly C-terminal of the vFLIP binding site.

Drug target optimization in chronic myeloid leukemia using innovative computational platform

Chronic Myeloid Leukemia (CML) represents a paradigm for the wider cancer field. Despite the fact that tyrosine kinase inhibitors have established targeted molecular therapy in CML, patients often face the risk of developing drug resistance, caused by mutations and/or activation of alternative cellular pathways. To optimize drug development, one needs to systematically test all possible combinations of drug targets within the genetic network that regulates the disease. The BioModelAnalyzer (BMA) is a user-friendly computational tool that allows us to do exactly that. We used BMA to build a CML network-model composed of 54 nodes linked by 104 interactions that encapsulates experimental data collected from 160 publications. While previous studies were limited by their focus on a single pathway or cellular process, our executable model allowed us to probe dynamic interactions between multiple pathways and cellular outcomes, suggest new combinatorial therapeutic targets, and highlight previously unexplored sensitivities to Interleukin-3.

Logic programming to predict cell fate patterns and retrodict genotypes in organogenesis

Caenorhabditis elegans vulval development is a paradigm system for understanding cell differentiation in the process of organogenesis. Through temporal and spatial controls, the fate pattern of six cells is determined by the competition of the LET-23 and the Notch signalling pathways. Modelling cell fate determination in vulval development using state-based models, coupled with formal analysis techniques, has been established as a powerful approach in predicting the outcome of combinations of mutations. However, computing the outcomes of complex and highly concurrent models can become prohibitive. Here, we show how logic programs derived from state machines describing the differentiation of C. elegans vulval precursor cells can increase the speed of prediction by four orders of magnitude relative to previous approaches. Moreover, this increase in speed allows us to infer, or 'retrodict', compatible genomes from cell fate patterns. We exploit this technique to predict highly variable cell fate patterns resulting from dig-1 reduced-function mutations and let-23 mosaics. In addition to the new insights offered, we propose our technique as a platform for aiding the design and analysis of experimental data.

Computing Clinically Relevant Binding Free Energies of HIV-1 Protease Inhibitors

The use of molecular simulation to estimate the strength of macromolecular binding free energies is becoming increasingly widespread, with goals ranging from lead optimization and enrichment in drug discovery to personalizing or stratifying treatment regimes. In order to realize the potential of such approaches to predict new results, not merely to explain previous experimental findings, it is necessary that the methods used are reliable and accurate, and that their limitations are thoroughly understood. However, the computational cost of atomistic simulation techniques such as molecular dynamics (MD) has meant that until recently little work has focused on validating and verifying the available free energy methodologies, with the consequence that many of the results published in the literature are not reproducible. Here, we present a detailed analysis of two of the most popular approximate methods for calculating binding free energies from molecular simulations, molecular mechanics Poisson-Boltzmann surface area (MMPBSA) and molecular mechanics generalized Born surface area (MMGBSA), applied to the nine FDA-approved HIV-1 protease inhibitors. Our results show that the values obtained from replica simulations of the same protease-drug complex, differing only in initially assigned atom velocities, can vary by as much as 10 kcal mol-1, which is greater than the difference between the best and worst binding inhibitors under investigation. Despite this, analysis of ensembles of simulations producing 50 trajectories of 4 ns duration leads to well converged free energy estimates. For seven inhibitors, we find that with correctly converged normal mode estimates of the configurational entropy, we can correctly distinguish inhibitors in agreement with experimental data for both the MMPBSA and MMGBSA methods and thus have the ability to rank the efficacy of binding of this selection of drugs to the protease (no account is made for free energy penalties associated with protein distortion leading to the over estimation of the binding strength of the two largest inhibitors ritonavir and atazanavir). We obtain improved rankings and estimates of the relative binding strengths of the drugs by using a novel combination of MMPBSA/MMGBSA with normal mode entropy estimates and the free energy of association calculated directly from simulation trajectories. Our work provides a thorough assessment of what is required to produce converged and hence reliable free energies for protein-ligand binding.

Primary and secondary dimer interfaces of the fibroblast growth factor receptor 3 transmembrane domain: characterization via multiscale molecular dynamics simulations

Receptor tyrosine kinases are single-pass membrane proteins that form dimers within the membrane. The interactions of their transmembrane domains (TMDs) play a key role in dimerization and signaling. Fibroblast growth factor receptor 3 (FGFR3) is of interest as a G380R mutation in its TMD is the underlying cause of ~99% of the cases of achondroplasia, the most common form of human dwarfism. The structural consequences of this mutation remain uncertain: the mutation shifts the position of the TMD relative to the lipid bilayer but does not alter the association free energy. We have combined coarse-grained and all-atom molecular dynamics simulations to study the dimerization of wild-type, heterodimer, and mutant FGFR3 TMDs. The simulations reveal that the helices pack together in the dimer to form a flexible interface. The primary packing mode is mediated by a Gx3G motif. There is also a secondary dimer interface that is more highly populated in heterodimer and mutant configurations that may feature in the molecular mechanism of pathology. Both coarse-grained and atomistic simulations reveal a significant shift of the G380R mutant dimer TMD relative to the bilayer to allow interactions of the arginine side chain with lipid headgroup phosphates.

Coarse-grain simulations reveal movement of the synaptobrevin C-terminus in response to piconewton forces

Fusion of neurosecretory vesicles with the plasma membrane is mediated by SNARE proteins, which transfer a force to the membranes. However, the mechanism by which this force transfer induces fusion pore formation is still unknown. The neuronal vesicular SNARE protein synaptobrevin 2 (syb2) is anchored in the vesicle membrane by a single C-terminal transmembrane (TM) helix. In coarse-grain molecular-dynamics simulations, self-assembly of the membrane occurred with the syb2 TM domain inserted, as expected from experimental data. The free-energy profile for the position of the syb2 membrane anchor in the membrane was determined using umbrella sampling. To predict the free-energy landscapes for a reaction pathway pulling syb2 toward the extravesicular side of the membrane, which is the direction of the force transfer from the SNARE complex, harmonic potentials were applied to the peptide in its unbiased position, pulling it toward new biased equilibrium positions. Application of piconewton forces to the extravesicular end of the TM helix in the simulation detached the synaptobrevin C-terminus from the vesicle's inner-leaflet lipid headgroups and pulled it deeper into the membrane. This C-terminal movement was facilitated and hindered by specific mutations in parallel with experimentally observed facilitation and inhibition of fusion. Direct application of such forces to the intravesicular end of the TM domain resulted in tilting motion of the TM domain through the membrane with an activation energy of ∼70 kJ/mol. The results suggest a mechanism whereby fusion pore formation is induced by movement of the charged syb2 C-terminus within the membrane in response to pulling and tilting forces generated by C-terminal zippering of the SNARE complex.

How lipid headgroups sense the membrane environment: an application of ¹⁴N NMR

The orientation of lipid headgroups may serve as a powerful sensor of electrostatic interactions in membranes. As shown previously by (2)H NMR measurements, the headgroup of phosphatidylcholine (PC) behaves like an electrometer and varies its orientation according to the membrane surface charge. Here, we explored the use of solid-state (14)N NMR as a relatively simple and label-free method to study the orientation of the PC headgroup in model membrane systems of varying composition. We found that (14)N NMR is sufficiently sensitive to detect small changes in headgroup orientation upon introduction of positively and negatively charged lipids and we developed an approach to directly convert the (14)N quadrupolar splittings into an average orientation of the PC polar headgroup. Our results show that inclusion of cholesterol or mixing of lipids with different length acyl chains does not significantly affect the orientation of the PC headgroup. In contrast, measurements with cationic (KALP), neutral (Ac-KALP), and pH-sensitive (HALP) transmembrane peptides show very systematic changes in headgroup orientation, depending on the amount of charge in the peptide side chains and on their precise localization at the interface, as modulated by varying the extent of hydrophobic peptide/lipid mismatch. Finally, our measurements suggest an unexpectedly strong preferential enrichment of the anionic lipid phosphatidylglycerol around the cationic KALP peptide in ternary mixtures with PC. We believe that these results are important for understanding protein/lipid interactions and that they may help parametrization of membrane properties in computational studies.

Accommodation of a central arginine in a transmembrane peptide by changing the placement of anchor residues

Both Trp and Arg in transmembrane protein domains make important interactions with lipids at the membrane/water interface, but at different depths. Derivatives of the designed peptide GWALP23, acetyl-GGALW(5)LALALALALALALW(19)LAGA-amide, with single Trp anchors, have proven useful for characterizing such interactions. Indeed, previous work revealed quite different effects emanating from Arg substitutions at positions 12 and 14 within GWALP23, with the R12 peptide exhibiting multiple positions and orientations with respect to DOPC bilayer membranes (Vostrikov et al. J. Am. Chem. Soc. 2010, 132, 5803-5811). To gain further understanding of the multistate behavior, we moved the Trp "anchor" residues to more outer positions 3 and 21 in GWALP23 itself, and in the R12 and R14 derivatives. The locations and orientations of the peptides with respect to lipid bilayer membranes of differing thickness were investigated by means of solid-state (2)H NMR spectroscopy, using labeled alanines, and coarse-grained molecular dynamics simulations. Interestingly, relatively intense and narrow (2)H resonances from selected backbone C(α) deuterons were observed over quite narrow ranges of frequency and sample orientation. The backbone resonances reflect dynamic complexities and at the same time provide important contributions for the analysis of peptide transmembrane orientation. With the Trp(3,21) anchors relatively far from the peptide and bilayer center, the results indicate significantly large apparent tilt angles, for example, close to 30° for the new R12 and R14 peptides with respect to the bilayer normal of DLPC membranes. The R12 side chain indeed is "rescued" to a stable position, where it is accommodated within the transmembrane helix, when the Trp anchors are moved outward and to another face of the helix. At the same time, the R14 side chain of transmembrane GW(3,21)ALP23 also retains a stable favored position.

Mechanism of bacterial signal transduction revealed by molecular dynamics of Tsr dimers and trimers of dimers in lipid vesicles

Bacterial chemoreceptors provide an important model for understanding signalling processes. In the serine receptor Tsr from E. coli, a binding event in the periplasmic domain of the receptor dimer causes a shift in a single transmembrane helix of roughly 0.15 nm towards the cytoplasm. This small change is propagated through the ∼22 nm length of the receptor, causing downstream inhibition of the kinase CheA. This requires interactions within a trimer of receptor dimers. Additionally, the signal is amplified across a 53,000 nm² array of chemoreceptor proteins, including ∼5,200 receptor trimers-of-dimers, at the cell pole. Despite a wealth of experimental data on the system, including high resolution structures of individual domains and extensive mutagenesis data, it remains uncertain how information is communicated across the receptor from the binding event to the downstream effectors. We present a molecular model of the entire Tsr dimer, and examine its behaviour using coarse-grained molecular dynamics and elastic network modelling. We observe a large bending in dimer models between the linker domain HAMP and coiled-coil domains, which is supported by experimental data. Models of the trimer of dimers, built from the dimer models, are more constrained and likely represent the signalling state. Simulations of the models in a 70 nm diameter vesicle with a biologically realistic lipid mixture reveal specific lipid interactions and oligomerisation of the trimer of dimers. The results indicate a mechanism whereby small motions of a single helix can be amplified through HAMP domain packing, to initiate large changes in the whole receptor structure.

To understand cell signalling events requires a physical model of the structure and behaviour of the signalling proteins involved. The methyl-accepting chemoreceptor proteins direct bacterial movement towards food sources and away from toxins. Based on experimental data we have built structural models of the serine chemoreceptor (Tsr) as a dimer, which is incapable of activating the downstream kinase CheA, and as a trimer of dimers, which can activate CheA. We have performed molecular dynamics simulation to reveal the behaviour of these two forms in a planar lipid bilayer and in a 70 nm diameter lipid vesicle with a mixture of lipids mimicking the E. coli inner membrane. We show that in isolation the dimers undergo a bending movement around the central HAMP domain, whereas the trimer-of-dimers model does not. Comparison with published experimental data suggests that these bending motions are real, and that they occur in the trimer of dimers only in response to ligand binding. Drawing together these observations with studies showing that the signalling event involves small piston motions in the transmembrane helices suggests that the bending motion is frustrated in the unliganded trimer of dimers, and that ligand binding induces bending by repacking the HAMP interface.

Quantized water access to the HIV-1 protease active site as a proposed mechanism for cooperative mutations in drug affinity

The development of resistance to different drugs remains a major problem for a wide range of infections. In particular, combinations of specific mutations, which individually demonstrate no effect, exhibit significant cooperativity. Here we show that changes to the energy of ligand binding in different resistant HIV-1 proteases are correlated with the creation of water binding sites in the active site. This correlation is conserved across two drugs (ritonavir and lopinavir). We propose that individual mutations induce changes in flap packing that are insufficient to allow water binding but in combination allow access, leading to the observed cooperative resistance.