ABSTRACT. Perturbation of cancer cells often leads to heterogeneous outcomes, in that most cells exhibit a dominant phenotype, but the rest appear resistant or hypersensitive to the perturbation.  If the penetrance of such a phenotype is heritably incomplete, then it becomes extremely difficult to decipher the upstream molecular events that heterogenize the population and cause response variability.  By combining quantitative measurements with dynamical models, systems approaches should be useful if provided with a core network of important biomolecules.  The daunting hurdle lies in identifying phenotype-relevant regulatory heterogeneities that define the network for penetrance at the single-cell level.  Here, I will introduce a new approach, called stochastic frequency matching (SFM), for elaborating the molecular networks upstream of incompletely penetrant phenotypes.  SFM identifies and parameterizes single-cell heterogeneities—which emerge after a uniform perturbation but before the appearance of a variable phenotype—to hone in on regulatory states corresponding to future penetrance.  For an onco-phenotype incompletely triggered by ErbB receptor tyrosine kinase signaling in 3D cultured breast epithelia, we implemented SFM using microarrays to uncover a network of critical nucleocytoplasmic regulators.

ABSTRACT. An individual tumor is composed of a heterogeneous collection of single cells with distinct molecular and phenotypic features, a phenomenon termed intratumoral heterogeneity. Intratumoral heterogeneity poses serious challenges to cancer treatment, motivating the need for combination therapies. We optimize drug combination by accounting for intratumoral heterogeneity through the analysis of single cell signaling perturbations when an individual tumor sample is screened by a drug panel. Mass Cytometry Time-of-Flight (CyTOF) is a high throughput single cell technology that enables the simultaneous measurements of multiple (>40) intracellular and surface markers at the level of single cells for hundreds of thousands of cells in a sample, analyzed pre- and post-treatment. We developed a computational framework, entitled DRUG-NEM, to analyze  CyTOF-based single cell drug perturbation data for the purpose of individualizing drug combinations. In its current implementation, DRUG-NEM optimizes the drug combinations by choosing the minimum number of drugs that produces the maximal desired intercellular effects based on nested-effects modeling. We demonstrate the performance of  DRUG-NEM for leveraging single cell perturbation data to identify  optimal drug combinations on tumor cell lines and primary leukemia samples.

ABSTRACT. The gut microbiome is now widely recognized as a dynamic ecosystem that plays an essential role in health and disease. While current sequencing technologies make it possible to estimate relative abundances of host-associated microbiota over time, the processes governing their temporal dynamics remain poorly understood due to significant bacterial diversity and interaction complexity in the gut. Consequently, as in other ecological systems, it is important to investigate global statistical relationships that describe microbiota dynamics. Here we perform such an ecological analysis using several high-resolution time series data sets from humans and mice, finding that microbiota dynamics can be described by robust scaling relationships spanning several orders of magnitude. Specifically, we observe power laws governing changes in population abundance, species local residence times, and mean versus variance scaling of individual taxa abundances over time. Interestingly, the observed patterns are highly similar to those describing multiple other ecological communities and economic systems, including temporal fluctuations of animal populations and performance of publicly traded companies. We find that these scaling relationships are altered in mice receiving different diets, and identify individual taxa whose dynamics deviate from overall trends in each group. Collectively, our results provide a systematic statistical framework for understanding complex ecological processes in the gut microbiome.

ABSTRACT. Acute kidney injury (AKI) is associated with a high risk of death. The death rate from AKI is much higher than that from prostate cancer, breast cancer, heart failure, diabetes. In response to renal injury, a complex wound-healing program is triggered to minimize the damage. If the damage is small, the function of the kidney is completely recovered. However, severe AKI or repeated episodes of AKI leads to kidney fibrosis and Chronic kidney disease (CKD). Interestingly, the long-term outcome of AKI patients after discharge from the hospital varies from person to person. While some of them are completely recovered from AKI, the others progress to fibrosis. However, the underlying mechanism is not fully understood.

Here, we first built a cell-cell communication mathematic model for renal homeostasis and fibrosis in response to injury. We found that depending on the level of injury, the outcomes can be death, perfective-adaptive AKI, imperfective-adaptive AKI, and maladaptive CKD. We verified this prediction in mice with different duration of ischemia-reperfusion injury (IRI) treatment. Furthermore, we demonstrated computationally and experimentally that imperfective adaptive AKI has a double-edge effect for the subsequent injury and significant increase the risk of fibrosis. On one hand, that imperfective adaptive AKI functions as a priming factor to reduce the death risk. On the other hand, it increases the fibrosis risk. Using combined mathematical modeling and mouse model studies, we not only recapitulates multiple-objective optimization of the renal system but also elucidates lots of seemly controversial experimental results. Most importantly, we proposed and tested a new treatment design by targeting on the dynamics of the Wnt signaling which can both reduce the death and fibrosis risk, especially under the circumstance of severe renal damage.

ABSTRACT. Cell cycle and metabolism are coupled networks. For example, cell growth and division require synthesis of macromolecules which is dependent on metabolic cues. Conversely, metabolites involved in nucleotide and protein synthesis are fluctuating periodically as a function of cell cycle progression. Although computational models of these networks are being developed for some time, to date no effort has been made to integrate these two systems in any organism. We aim to investigate cell cycle robustness by generating the first multi-scale model that integrates cell cycle with metabolism, and investigating their bidirectional regulation. Connections among these two biochemical networks have been recently elucidated in budding yeast. However, high-throughput and manually curated studies point at many more physical interaction, which relevance for precise cell cycle timing remains unknown. A framework is presented that integrates a Boolean cell cycle model with a constraint-based model of metabolism, incorporating mechanistic and high-throughput interactions. Directionality and effect are incorporated for the mechanistic interactions. Conversely, as this information is unknown for the high-throughput interactions, an informed optimization algorithm has been developed to generate models that can incorporate it iteratively. To verify the results of the informed optimization algorithm against metabolomic data, changes in flux through a number of metabolic pathways are compared to metabolic pathway enrichment time-series. The multi-scale model predicts expected changes in a number of pathways, ranging from amino-acid to pentose phosphate to lipid metabolism. Many model variants that differ in number and directionality of interactions robustly predict the effect of definite cell cycle-metabolism pairs. Furthermore, the integrative model shows a temporal export of acetate, pyruvate and alanine, reminiscent of yeast metabolic oscillations. Altogether, our multi-scale framework is able to integrate computer models of biological networks with high-throughput data, to capture the functional connectivity among their elements that ultimately results in systems robustness.

ABSTRACT. Multicellular biological organisms possess the extraordinary ability to grow and maintain intricate body shapes and forms. However, the processes regulating the growth of exact shapes in a developing organism are not well understood due to the non-linearity of gene regulatory networks and the complex, systemic interactions between tissue and regulatory signals. It is now clear that robust mathematical approaches combined with automated reverse-engineering methods are necessary to describe, infer, and understand these dynamic processes controlling tissue growth and shape formation. To this end, we have developed a mathematical framework based on continuous dynamics and a machine learning methodology based on evolutionary computation to accurately reverse-engineer and predict the growth and regulation of multicellular biological forms and shapes. Our formulation is based on diffusion-advection-reaction partial differential equations, permitting the fast and efficient simulation of multicellular growth. In a continuous fashion, this mathematical framework models the proliferation of cells, their chemotactic migrations, and the adhesive forces between them, in addition to the genetic regulatory mechanisms controlling these processes. Importantly, our computational methodology can automatically infer de novo these complex dynamic models directly from quantitative experimental data, including the precise genes, their interactions, and model parameters necessary and sufficient to develop and maintain specific biological forms and shapes. As a proof of concept, we have applied our novel framework to infer the genetic network controlling the allometric growth in regenerating planarian worms. Using a dataset of planarian experiments and their resultant morphologies, we have discovered a suitable, dynamic genetic model that accurately recapitulates the regeneration of shape and form after surgical amputations. This work paves the way for the understanding and reverse engineering directly from experimental data the dynamic regulation of biological shapes, an essential step towards the much sought-after quantitative and predictive models in developmental and regenerative biology.

ABSTRACT. After developing from pluripotent embryo stem cells, multipotent intestinal stem cells continuously differentiate into mature functional cells (e.g. Enterocytes, Goblet cells, Enteroendocrine cells and Paneth cells). In presence of oncogenetic mutations, intestinal cells might evolve into cancer cells that are killed during the treatment with chemotherapy drugs. The self-renewal, differentiation and death of intestinal cells are complex dynamical processes that are orchestrated by the molecular control networks in these cells. Previous work has discovered valuable pieces of the molecular control network, but an integrated, systemic understanding of the overall dynamical process remains unknown. Without such an integrated understanding of how the cells normally coordinate their self-renewal and differentiation, it is hard to understand how these processes are disputed in diseases such as Inflammatory Bowel Disease (IBD) and cancer. Consequently, it is hard to rationally design optimal clinical protocols to treat these diseases. In order to cope with this challenge, we have converted several key pathways into computational models. These models have been constrained with available data reported in the literature, and the novel predictions generated by these models are being tested in intestinal enteroids. Testing of these predictions will indicate how well the system is currently understood and push forward the boundary of our understanding.

ABSTRACT. In cancer, nearly all aspects of cancer pathophysiology, including cancer initiation, development, progression and metastasis are driven by the dysregulation of one or more signalling networks. These signalling networks consist of sets of genes that are dynamically organised. In response to a perturbation the activity of the network changes over time. In this way, signalling pathways exert finely tuned control over cell fate decisions that ultimately determine the behaviour of cancer cells. However, we barely understand how these dynamic activation patterns are shaped by differences between cells and patients: How does cell-to-cell variability affect drug responses? How does the genomic background of patients affect the network's input/output behaviour? Here, I will address these questions focusing on the on the dynamic network around the JNK and p53 stress- and DNA-damage responses in neuroblastoma and breast-cancer. Firstly, I will present a generally applicable method for integrating tumour data into patient-specific dynamic models of cancer signalling. Secondly, I will show how patient-specific modelling of the JNK-p53 response network can stratify neuroblastoma patients. Thirdly, I will present our progress on modelling the cell-to-cell variability during the development of drug resistance. Both our theoretical and experimental results indicate that so-called cell-ensemble modelling can be used model the selective pressure of chemotherapy on cancer-cell populations.

[1] Fey, D. et al. (2015) "Signaling pathway models as biomarkers: Patient-specific simulations of JNK activity predict the survival of neuroblastoma patients" Science Signalling 8(408): ra130. [2] Kim, J., B. Schoeberl (2015). "Beyond static biomarkers - The dynamic response potential of signaling networks as an alternate biomarker?" Science Signaling 8(408): fs21. [3] Fey D, Kuehn A, Kholodenko BN (2016) "On the personalised modelling of cancer signalling" IFAC-PapersOnLine 49(26):312-317

ABSTRACT. Prostate cancer is a leading male cancer in the world, especially in industrialized countries. Nevertheless, proteome-wide analysis of human prostate cancer cells remains scarce, hindering our understanding of prostate cancer development and progression at the systems level. Here, we report the most comprehensive proteogenomic profiling study of human prostate cancer cells undertaken to date. A total of 11,759 protein groups, corresponding to 10,561 human genes, were identified with a false discovery rate of ≤1%. Of these, 51 “missing proteins” were identified with high-stringency mass spectrometry evidence. Absolute quantification of all protein groups suggested that the abundance of proteins identified as relevant to prostate cancer spans five orders of magnitude, and that the identified “missing proteins” are of very low abundance (median < 2,000 copies per cell). Owing to the relatively high (median = 45.0%) protein sequence coverage, protein isoforms encoded by 755 genes were distinguished. In addition, our deep proteogenomic analysis identified 291 proteogenomic peptides, including 227 novel peptides derived from 165 protein variants. Through targeted mass spectrometry quantification of proteogenomic peptides of interest in seven prostate cell lines, the N-terminal extension form of ZDHHC20 was found to be much more abundantly expressed in bone metastasis-derived PC3 cells than in other cell lines. The proteogenomics data are expected to be a valuable resource for proteogenomics and prostate cancer research.

ABSTRACT. Epithelial-to-Mesenchymal Transition (EMT) and its reverse Mesenchymal-to-Epithelial Transition (MET) often play crucial roles in cancer metastasis and drug resistance. Recent reports highlight that EMT and MET are not ‘all-or-none’ processes, instead cells can attain a hybrid epithelial/mesenchymal (E/M) phenotype. But, a hybrid E/M phenotype has been tacitly assumed to be ‘metastable’ that can be attained only transiently en route to EMT/MET, and remains poorly characterized. Rapid progress in mapping the regulatory networks for EMT/MET has enabled developing computational systems biology models to characterize a hybrid E/M phenotype. Here, using mechanism-based mathematical modeling, we identify a set of ‘phenotypic stability factors’ (PSFs) – OVOL2 and GRHL2 – that can help maintain cells in a hybrid E/M state. Next, we identify H1975 cells as stably maintaining a hybrid E/M state over multiple passages, and validate the role of these PSFs experimentally. We show that the knockdown of these PSFs that act as a ‘brake’ on full EMT drives cells to a fully mesenchymal phenotype. Finally, we devise a statistical model built upon gene expression profiles that can quantitatively predict where a given sample lies on a scale of 0 (fully epithelial) to 2 (fully mesenchymal). Intriguingly, GRHL2 and OVOL2 were identified among the top predictors that could resolve a hybrid E/M phenotype, through an unsupervised screening, thereby reinforcing their suggested roles as PSFs. This model can recapitulate the experimentally observed behavior for multiple scenarios such as EMT induction, and unravels the association of a hybrid E/M phenotype with poor clinical outcomes across multiple tumor types. Collectively, our integrated theoretical-experimental approach enables a quantitative understanding of the role of a hybrid E/M state in tumor progression, and reinforces the emerging notion that cells in a hybrid E/M state may be more aggressive than cells in a full EMT state.

ABSTRACT. Understanding the basis of complex cellular diseases such as cancer requires a multifaceted picture of cellular function and dysfunction. Here, we describe recent work integrating proteomic, expression, interaction and drug perturbation data to provide a broader view of the architecture of the breast cancer kinome and its behavior in response to targeted perturbations. Specifically, multiplexed small molecule inhibitors covalently bound to Sepharose beads (MIBs) were used to capture functional kinases in luminal, HER2-enriched and triple negative, basal-like and claudin-low breast cancer cell lines and tumors. Kinase MIB-binding profiles at baseline without perturbation were significantly uncorrelated to transcript abundance for many kinases and proteomically distinguished the four breast cancer subtypes. Understudied kinases were highly represented in the MIB-binding taxonomies and shown to be integrated in kinase signaling subnetworks with characterized kinases. Furthermore, we were able to identify regions of the kinome that are poorly characterized and/or poorly targeted by existing inhibitor therapies. Computationally it was possible to define subtypes using profiles of less than 50 of the more than 300 kinases bound to MIBs that included understudied as well as metabolic and lipid kinases. MIB-binding profiles readily defined subtype-selective differential adaptive kinome reprogramming in response to targeted kinase inhibition. Comprehensive MIBs-based capture of kinases provides a unique proteomics-based method for defining functional kinome dynamics and subnetworks in cells and tumors that integrates poorly characterized kinases of the understudied kinome that is not possible using genomic strategies. Together, this systems view of the kinome presents potential opportunities for disease classification, identification of potential drug targets and the broader design of rational combination therapies.

ABSTRACT. Quantitative description and classification of aberrant pathway activities in tumors can inform design of effective cancer therapies as most targeted agents inhibit tumor cell proliferation by blocking oncogenic signaling. However, development of such therapies has been a challenge since multiple oncogenic pathways can be co-activated in a given tumor and the pathway activation patterns vary substantially even within similar tumor types. Therefore, tumor-specific combination therapies are required to block multiple aberrant pathways. Our strategy involves an algorithmic approach to classify actionable oncogenic pathway signatures in large tumor/cell line cohorts and experimental testing of combination therapies specific to each oncogenic signature. For this purpose, we developed an integrated bioinformatics pipeline and a high-throughput experimental validation platform. We analyzed the expression and phosphorylation level changes of > 200 proteins in > 7000 tumor samples available from the TCGA project and > 600 cell lines. The phosphoproteomic data was collected at MDACC using the reverse phase protein array technology. We employed an iterative machine-learning algorithm that couples feature selection with clustering to identify a combination of discriminant and actionable protein biomarkers shared within each tumor subcohort. We identified the actionable targets within each subcohort specific oncogenic signaling signature in collaboration with domain experts and through database searches. Our results uncovered potentially actionable combinations of protein targets shared among subcohorts of tumors and cell lines over a large number of lineages. We are testing our predictions with drug combinations in cell lines that share the target oncogenic signaling signature. The most promising combination therapy candidates will be tested in patient derived xenograft models. We expect our strategy will expedite the global efforts for precision therapy development as the experimentally validated drug combinations will be nominated for basket clinical trials at MDACC and elsewhere.

ABSTRACT. Interferons (IFNs) signalling is a key mechanism to coordinate antiviral, anti-proliferative and immunomodulatory effects (1). A substantial amount of molecular details is known regarding IFNs signalling pathways, even though understanding, how information about complex mixture of IFNs is processed and translated into distinct cellular responses remains elusive (5,6). A good illustration is a sensitising effect of IFN type-I. Although the presence of this phenomenon is well known, its impact on signalling fidelity and biochemical mechanism that lead to these changes has not been recognised so far. Our experimental studies on IFNs signalling on mouse embryonic fibroblasts have shown that prior exposure to IFN type-I modify cellular response to IFN type-II stimulation. Precisely, information-theoretic analysis indicate higher sensitivity of pre-stimulated cells to the presence of IFN type-II in the intercellular environment. However, due to the complexity of signalling networks identification of origins of this mechanism cannot be addressed with solely experimental methods (7). Here we propose integration of high-throughput experimental single-cell measurements with a stochastic modelling in order to provide better understanding of mediation between IFNs type-I and -II signalling. Our solution is based on analysing intrinsic and extrinsic sources of heterogenous cellular response using unscented transformation and Sequential Monte Carlo methods (8,9). We have shown that origins of increasing sensitivity of pre-stimulated cells are changes in the initial cellular concentration of signal transducer and activator of transcription (STAT) proteins. Deciphering the mechanism that lead to more sensitive cellular response informs further research on novel therapeutic and diagnostic strategies to utilise the clinical potential of IFNs (5). 1. R.J.Critchley-Thorne et al. (2009) Impaired interferon signaling is a common immune defect in human cancer, Proceedings of the National Academy of Sciences.National Acad Sciences, 106:9010–5. 2. D.S.Aaronson (2002) A Road Map for Those Who Don't Know JAK-STAT, Science, 296:1653–5. 3. I.M.Kerr et al. (1994) Jak-STAT pathways and transcriptional activation in response to IFNS and other extracellular, Science. American Association for the Advancement of Science, 264:1415–21. 5. B.S.Parker et al. (2016) Antitumour actions of interferons: implications for cancer therapy, Nat Rev Cancer, 16(3):131-144.. 6. L.Zitvogel et al. (2015) Type I-interferons in anticancer immunity, Nature Rev Immunol, 7:405-414.. 7. B.N.Kholodenko (2006) Cell-signalling dynamics in time and space, Nat Rev Mol Cell Biol, 7:165–76. 8. T.Toni, B.Tidor (2013) Combined Model of Intrinsic and Extrinsic Variability for Computational Network Design with Application to Synthetic Biology, PLoS Comput Biol, 9:e1002960. 9. S.Filippi et al. (2016) Robustness of mek-erk dynamics and origins of cell-to-cell variability in mapk signaling, Cell reports, 15:2524-2535.

ABSTRACT. Humans are constantly exposed to complex mixtures of environmental chemicals. Several ongoing efforts seek to increase our knowledge of chemical effects. For example, the EPA/NIH-funded ToxCast and Tox21 initiatives monitor the effect of chemicals on selected proteins using high-throughput screening assays. Toxicogenomic databases store gene expression profiles after chemical exposure. However, each dataset probes a different dimension of cell's response. Moreover, these experiments ignore complex networks through which proteins interact and computational methods to integrate these data are underdeveloped. These major barriers limit their usefulness of these data.
We propose computing toxicant signaling networks to address these challenges. For a chemical, such a network is composed of regulatory, signaling and physical interactions connecting proteins perturbed as a result of exposure to that chemical. We describe algorithms that connect the responding proteins in ToxCast/Tox21 assays in the context of the underlying network of regulatory and physical interactions. For well-studied chemicals, our networks are enriched in biological processes that are known to be perturbed by the chemicals. Toxicant signaling networks promise to reveal important intermediate proteins that have not have been tested and physiological processes that have not been previously implicated in connection with the chemical.

ABSTRACT. The proper transition of mammalian cells between quiescence and proliferation is critical to tissue homeostasis and differentiation. Its deregulation is commonly found in many human diseases. Cryptochrome (Cry) is a transcriptional factor that is responsible for generating the negative feedback loop to maintain the circadian clock; its deregulated expression has been found to affect many physiological processes. Here, we focused on the role of circadian clock protein Cry in regulating the depth of quiescent cells. Experimentally, we showed that the overexpression or stabilization of Cry resulted in upregulation of c-MYC in quiescent rat embryo fibroblasts; however, cells did not go to a ‘shallower’ quiescent state as expected but a 'deeper' state. That is, a higher instead of lower serum concentration was required to drive cells out of quiescence in cells with upregulated Cry and Myc activities. Through systematic modeling of an array of possible regulatory network topologies between Cry and the Rb-E2F bistable switch that controls the quiescence-to-proliferation transition, we found that the ‘deeper’ quiescent state in response to increased Cry activity may be explained by upregulation of CDK inhibitors which counteract the c-MYC effect on the activation threshold of the Rb-E2F bistable switch. We further confirmed this model prediction in follow-up experiments. Our findings suggest a mechanistic role for circadian clock protein Cry in modulating the depth of cellular quiescence, which may have implications in varying potentials of tissue repair and regeneration in different times of the day.

ABSTRACT. In the last few years, several models of the spatio-temporal replication program in eukaryotic cells were proposed in the literature. In these models, the frequency of new replication origin firing per length of unreplicated DNA along the S-phase, I(t), is a fundamental quantity characterizing the DNA replication program dynamics. The I(t) curves have been shown to present a universal bell shape in eukaryotes (Goldar, 2009), increasing until it reaches a maximum after mid-S-phase, and then decreasing to reach zero at the end of S-phase. Analytical modeling in Xenopus (Gautier, 2009) explained the initial increase by the progressive import of a factor required for origin activation, and the final decrease by a sub-diffusive motion of this factor. A simulation-based model (Goldar, 2008) also concluded for the need of a constant import rate of initiation factors and explained the final decrease by a coupling between the frequency of initiation and the fork density. A recent simulation-based model in human cell (Gindin, 2014) was able to reproduce the I(t) bell shape assuming an increasing reactivity of the initiation factors. After a careful analysis of these different models that enlighten their advantages and limitations, we propose a model where (i) origin firing results from a second-order reaction between unreplicated origins and unbound initiation factors and (ii) initiation factors remain bound to the replication forks until replication termination. This simple model that we corroborated by 3D simulations, fully accounts for I(t) universal bell shape. The initial increase results from the recycling of the initiation factors and the final decrease is obtained when the time between two initiations become shorter than the time to replicate the mean distance between two origins. This modeling provide some prediction of the maximum value of I(t) as a function of the replication fork speed and the density of potential origins.

Goldar, Arach, Marie-Claude Marsolier-Kergoat, and Olivier Hyrien. “Universal Temporal Profile of Replication Origin Activation in Eukaryotes.” PLOS ONE 4, no. 6 (June 12, 2009): e5899. doi:10.1371/journal.pone.0005899.

Goldar, Arach, Hélène Labit, Kathrin Marheineke, and Olivier Hyrien. “A Dynamic Stochastic Model for DNA Replication Initiation in Early Embryos.” PLOS ONE 3, no. 8 (August 6, 2008): e2919. doi:10.1371/journal.pone.0002919.

Gauthier, Michel G., and John Bechhoefer. “Control of DNA Replication by Anomalous Reaction-Diffusion Kinetics.” Physical Review Letters 102, no. 15 (April 16, 2009): 158104. doi:10.1103/PhysRevLett.102.158104.

Gindin, Yevgeniy, Manuel S Valenzuela, Mirit I Aladjem, Paul S Meltzer, and Sven Bilke. “A Chromatin Structure‐based Model Accurately Predicts DNA Replication Timing in Human Cells.” Molecular Systems Biology 10, no. 3 (March 28, 2014). doi:10.1002/msb.134859.

ABSTRACT. G protein-coupled receptors (GPCRs) are critical cell signaling molecules that also comprise the largest class of therapeutic drug targets. GPCRs are well known to signal upon ligand binding via production of the second messenger molecule, cyclic AMP (cAMP), at the plasma membrane, but recent evidence has indicated that various GPCRs may also signal after ligand-induced internalization. We investigated the functional consequences of compartmentalized GPCR signaling, using the beta2-adrenoceptor (β2-AR) as a model system. Global profiling of β2-AR activation identified a core set of transcriptional target genes, and revealed that endocytosis is required for the full repertoire of downstream cAMP-dependent transcriptional responses. We then developed and applied an orthogonal optogenetic approach to definitively establish that the location of cAMP production is indeed the critical variable controlling the transcriptional response to β2-AR signaling. Further, we found that β2-AR endocytosis may provide a signaling ‘checkpoint’ that enables cells to respond uniformly to chemically distinct ligands acting on the same receptor while limiting spurious responses from non-cognate ligands. We are currently investigating the molecular mechanisms underlying location bias of signaling through functional genome-wide CRISPR interference-based gene silencing. Altogether, our findings establish a novel principle for functional diversification of signaling, based on the location of second messenger production, which underlies cellular discrimination of chemically distinct ligands.

ABSTRACT. Cells live in ever-changing environment and need to reliably receive, decode, and transmit information of extracellular signals such as their strength and duration, so cells can respond appropriately. It is still an open question how cells distribute and integrate information through an intertwined network, and generate either adaptive or sustained responses.

TGF-beta is a multifunctional cytokine that regulates many important cellular processes. For a long time people are puzzled by the paradoxical roles of TGF-beta, for example, it functions as both tumor suppressor and tumor growth/metastasis promotor. Through integrated modeling and quantitative measurements we uncovered a novel nested relay mechanism that allows cells to read the temporal information of TGF-beta, and generate Snail1 expressions with different temporal profiles. We expect the mechanism to be general for signal transduction. We also identified a modified positive feedback loop network motif that can achieve seeming opposite tasks of being robustly against noises and having fast response time. The motif again appears frequently in networks regulating processes such as cell differentiation and immune responses.

Given the importance of TGF-beta signaling and the amount of existing studies on it, the new discovery is rather surprising. We argue that this work is a good example why an integrated approach including modeling and quantitative experimental studies is needed to study complex cellular process. Some key features of the system, such as the temporal and spatial evolution of different phosphorylation forms GSK3, have escaped from observation in the past. The work also emphasizes the importance of “temporal dynamics” on understanding biological systems and on biomedical interventions.

ABSTRACT. Directed or “polarized” growth and the detection of chemical gradients are two fundamental cellular processes. Here we combine mathematical modeling and analysis with various experimental approaches to investigate the molecular mechanisms that underline both processes during the mating response of yeast. Our analysis reveals a novel method for gradient sensing and insight into the biochemical mechanisms that either ensure the establishment of a unique site for polarized growth or allow multiple sites to coexist.

ABSTRACT. The downstream signaling through insulin receptor and insulin-like growth factor I (IGF1) receptor are different in normal and disease states. Under normal conditions, IGF1 is a proliferation agent whereas insulin is one of the regulators of glucose homeostasis. Multiple investigations have shown evidence of similar functions and associations of the two hormones in cancer progression, cell proliferation, and evasion of apoptosis. Here, new insights on the mechanisms of differential MAPK and Akt activation are revealed by an iterative systems biology approach. The mechanistic network modeling here provided a framework to elucidate experimental targets downstream of two receptors, which were treated as indistinguishable in previous models. The model included cascades of both mitogen-activated protein kinase (MAPK) and Akt signaling, as well as the crosstalk and feedback loops in between. The parameter perturbation scanning yielded new experimental hypotheses on how differential responses of MAPK and Akt originate. Complementary to our previous efforts, the results suggested that SOS activation through IRS is critical to inducing greater MAPK activation in IGF1 stimulated cells. Secondly, the negative feedback from Akt to IRS is predicted to play a key role in the enhancement of Akt activation, seen in IGF1 stimulated cells with E-Cadherin knock-down. Third, the feedback from ribosomal protein S6 kinase (p70S6K) on IRS is predicted to differentially affect Akt activation under IGF1 and insulin stimulated cells. The experimental validation of the last prediction showed that there indeed is a difference in the regulation of Akt activity in response to different stimuli. Our results, computational and experimental together, showed the importance of interactions of the adaptor protein IRS in activation/inhibition of these specific cascades in breast cancer cells.

ABSTRACT. Progression throughout the eukaryotic cell cycle is governed by waves of cyclin-dependent kinase (cyclin/Cdk) activity, which rise and fall at a specific timing. By integrating computation and experimentation, we have recently unraveled a mechanism involving the Fkh2 transcription factor in the synchronization of waves of mitotic (Clb) cyclin expression in budding yeast. Fkh2 promotes transcription of both CLB3 and CLB2 mitotic genes, and it is activated by Clb/Cdk1-mediated phosphorylation. Conversely, the stoichiometric inhibitor Sic1 and the Anaphase Promoting Complex (APC) abolish the Cdk catalytic activity by binding to and inhibiting the Clb/Cdk complexes, and degrading their regulatory (Clb) subunits, respectively. A trade-off exists between (i) the accuracy by which system’s complexity is modeled and (ii) the computational cost to analyze a model. Here we explore one aspect of this trade-off, in order to identify a minimal cell cycle model able to generate sustained oscillations in the form of limit cycles. We updated our previously published Clb/Cdk1 network by including new interactions, yielding three structurally different kinetic models. Furthermore, one of the models is simulated under a quasi-steady state assumption, which implies a fast equilibrium of ternary complex formation between Clb/Cdk1 complexes and Sic1. Additionally, the number of parameters is progressively reduced, allowing the models to be manageable for comprehensive parameter scans. For each step of model reduction, we highlight the model predictions for transient oscillations, illustrate the ability to generate sustained oscillations in the form of limit cycles, and perform a sensitivity analysis on the limit cycle. Through perturbation analyses of the limit cycles, we identify crucial parameters in the models that control properties of oscillations, which are conserved throughout progressive model reductions. Altogether, we uncover a definite network motif able to generate autonomous oscillations of a minimal cell cycle network with respect to the formation of Clb cyclin waves.

ABSTRACT. As a well-known tumor suppressor, p53 plays a key role in cell fate decision in cellular response to various stresses. It has been reported that p53 dynamics are stimulus-dependent, and different modes of p53 dynamics lead to distinct cellular outcomes. In this work, we developed a network model of p53 to explore how p53 dynamics modulates cell fate in the DNA damage response. We found that p53 exhibits three modes of dynamics in response to DNA double-strand breaks (DSBs). Upon mild damage, p53 shows persistent pulses and induce cell cycle arrest. For severe damage, p53 exhibits two-phase dynamics: a series of pulses appears to induce cell cycle arrest in the early phase, and p53 levels switch to high constant levels to trigger apoptosis in the late phase. Moreover, the number of pulses in the first phase drops with increasing DNA damage. When DNA damage is extremely severe, p53 directly rises to high levels and induces apoptosis quickly. We proposed that the alternation in the predominance of positive and negative feedback loops results in the transition of p53 dynamics from pulses to high constant levels. Our results suggested that different modes of p53 dynamics directs cells toward distinct cellular outcomes.

ABSTRACT. The Spindle Assembly Checkpoint (SAC) is a complex surveillance mechanism that ensures the fidelity of chromosome segregation during mitosis. To probe the mechanism by which SAC signal is modulated during mitosis, the Joglekar lab has engineered an ectopic, kinetochore-independent SAC activator, the “eSAC”, that stimulates the signaling cascade of the SAC, and arrests cells in mitosis. The time spent in mitosis by individual cells in the presence of eSAC activation shows nonlinear dependencies on the concentration of eSAC (dose-response curves).

By combining a previous SAC model (He et. al., PNAS, 108(24), 2011) with details of eSAC, we have developed a stochastic mathematical model which explains the complex dose-response characteristics of eSAC. By doing so, we have also arrived at a plausible mechanism of how the SAC signal is maintained as the number of unattached kinetochores decreases over time during mitosis. Briefly, in the presence of multiple unattached kinetochores, SAC proteins are sparsely distributed among KNL1 molecules on unattached kinetochores. This leads to a steady SAC signal generated by multiple, weakly signaling kinetochores. In contrast, when only a few kinetochores are unattached, multiple SAC proteins are recruited per KNL1. By a synergetic effect, a small number of unattached kinetochores produce a strong anaphase-inhibitory signal. Together, these mechanisms ensure a steady SAC signal for the full duration of mitosis.

ABSTRACT. Quiescence is a “sleep-like” non-proliferative cellular state. Reactivating quiescent cells (e.g., fibroblasts, lymphocytes, and stem cells) to proliferate is fundamental to tissue repair and regeneration. Often described as the “G0 phase”, quiescence is in fact not a homogeneous state. As cells remain quiescent for longer durations, they move progressively “deeper” into quiescence, exiting from which requires prolonged and stronger growth stimulation. Nevertheless, deep quiescent cells can still re-enter the cell cycle under physiological conditions, distinguishing them from senescent cells. Underlying mechanisms of quiescence heterogeneity remain an enigma, and represent a currently underappreciated layer of complexity in growth control. Previously, we have shown that the retinoblastoma (Rb)-E2F pathway functions as a bistable switch that controls the all-or-none transition from quiescence to proliferation. Here by coupling modeling and single-cell measurements, we show that quiescence depth is controlled by the serum threshold to activate the Rb-E2F switch, which is in turn modulated by different Rb-E2F pathway components with different efficacies. We also found that the Rb-E2F activation threshold can be modulated by other cellular pathways including Notch and circadian pathways as well as metabolic responses. Such pathways crosstalk with the Rb-E2F pathway and together they form a quiescence regulatory network that determines the final quiescence states. Further elucidating the control of quiescence heterogeneity should provide the basis for future therapeutic strategies against hypo- and hyper-proliferative diseases by counteracting abnormal quiescence depths in diseased cells.

ABSTRACT. Non-small cell lung carcinoma (NSCLC) is one of the leading causes of death for both men and women. Early stage NSCLC is often treated with surgical resection alone, yet a significant fraction of patients’ progress with local-regional and distant metastatic disease. An active research area is to identify which patients are likely to progress based on gene expression signature of their initial tumor. To date, many different signatures have been proposed. These signatures are likely related through underlying molecular processes of the tumor microenvironment. A better understanding of the tumor microenvironment, achieved by dissecting tumor-stromal interactions, promises a more robust interpretation and evaluation of genomic signatures of prediction and prognosis. Recent studies have demonstrated that stromal cells within the tumor collaborate to drive cancer progression. To elucidate the role of stromal cells, we want to pinpoint the highly active biological networks within each cell-type. Current approaches of the tumor microenvironment include simply finding differentially expressed genes. The gene lists returned often have high false positive rates and are hard to interpret. Instead, we developed a new algorithm to extract protein-protein interaction subnetworks using gene expression data. This algorithm was applied to both RNA-Seq data from individually sorted stromal and epithelial cells within the NSCLC tumor microenvironment and from the Cancer Genome Atlas (TCGA). The resulting subnetworks elucidated cell-type specific functions within the tumor microenvironment with the goal of providing new insights into tumor-stromal mechanisms associated with tumor progression. This will aid with the interpretation of prognostic signatures and identification of tumor-stromal interaction drug targets.

ABSTRACT. When chemotherapeutics are applied to tumor cells with the same or similar genotypes, some tumor cells are killed while others survive. This so called fractional killing can contribute to drug resistance in some cancer types. Drugs such as cisplatin not only activate p53 to induce cell death, but they also promote the expression of apoptosis inhibitors such as cIAP that protects cells through an incoherent feedforward loop. Due to this regulation, some cells activate their p53 quickly and undergo apoptosis, and cells that activate p53 slowly, survive. The incoherent feedforward loop along with the essential role of p53 activation timing makes fractional killing a complex dynamical challenge, which is difficult to understand using intuition alone. To better understand this process, we have constructed a representative model by integrating the control of apoptosis with the relevant signaling pathways. The model was then trained to recapture the observed properties of fractional killing and analyzed with nonlinear dynamical tools. Our analysis suggested a simple formula for fractional killing, where cell fate is a function of the bifurcation geometry and cell the trajectory. This formula predicts that cell fate can be altered in three possible ways: alteration of bifurcation geometry, alteration of cell trajectories or both. These predicted categories are able to explain existing strategies known to combat fractional killing and allow us to design novel strategies. In conclusion, we propose that our simple formula can be used as a powerful tool for understanding and designing novel strategies that can combat fractional killing.

ABSTRACT. Many studies in the systems biology of cancer treat cell signaling networks as “molecular circuits” and aim to identify key elements responsible for disease. Alternatively, one may consider normal and cancer cells as different dynamic states of the molecular network and to treat disease onset as a global “phase transition”. Following this premise, we focus on dynamic network entropy – a statistical measure related to stability and robustness of network states. We investigate whether it can be used to address the problem of drug response in cancer cells. We assume that cellular homeostasis corresponds to a dynamically stable state of the molecular network. In cells where this state is not very robust the networks are likely to respond to drug-induced perturbations by undergoing a significant irreversible transition, resulting in a “drug-sensitive” phenotype. To the contrary, robust networks will relax back to the original state, exhibiting drug resistance. The likelihood of such transition can be approximated by the so-called “network entropy rates” which characterize network’s ability to dissipate perturbations. To test this hypothesis, we calculated network entropy rates for a set of approximately 100 cancer cell lines using gene expression profiles and curated data on protein-protein interactions. We investigated relations between entropy rates and sensitivity to several cancer drugs and found that for multiple network modules they significantly correlate with IC50 values of drug response. Among those we found pathways related to targets of the drugs used in the study, indicating that dynamic stability of these elements is one of the determinants of drug sensitivity. We also demonstrated that entropy rates can effectively discriminate between sensitive and resistant cell lines and therefore can be used as variables in building predictive models of drug response.

ABSTRACT. Intravenous delivery of pharmacological ascorbate has shown to be a promising adjuvant in the treatment of patients with pancreatic ductal adenocarcinomasi. Administered as a series of infusions, pharmacological ascorbate generates high fluxes of extracellular hydrogen peroxide (H2O2), which is toxic to some cancer cells while not affecting normal cellsi. In vitro studies also indicate that cancer cells have a wide range of susceptibility to pharmacological ascorbateii and subsequently to extracellular H2O2. The resulting H2O2 flux into cancer cells is believed to accumulate differently within the intracellular space when compared to normal cells. We hypothesize that internal H2O2 concentration has a steady-state value that is significant for cell susceptibility and independent of cell type. Although this has been alluded to, this value has yet to be quantified. Quantification of resulting intracellular H2O2 concentrations during pharmacological ascorbate therapy is necessary for understanding the relevant transport and reactions associated with this cancer treatment. Here, we develop a lumped parameter model for intracellular H2O2 quantification for any given cell type using experimental procedures to obtain necessary cell-specific parameters. We show the variations in cell parameters (i.e. membrane permeability via peroxiporin expression, catalase rates, cell size etc.) for various cells types are significant enough to alter the internal H2O2 concentration during ascorbate therapy, thereby impacting cell susceptibility. Further we found that the intracellular H2O2 is highly sensitive to the membrane permeability of the plasma and peroxisome membranes. In addition, flow cytometry displays elevated peroxiporin expression on ascorbate-susceptible cells. These results imply that plasma membrane permeability to H2O2 via peroxiporins is an important factor in the surviving fraction outcomes for susceptible cell lines. This work suggests that relatively high H2O2-permeability of the plasma membrane of cells, either via peroxiporins or other mechanisms, is a critical factor in the success of therapeutic ascorbate in cancer treatment.

ABSTRACT. Use of methotrexate (MTX), a widely used anti-cancer drug is being limited by the emergence of resistance. MTX inhibits several enzymes in the folate pathway to different extents. Despite a large number of studies, a quantitative understanding of target pathway dynamics in resistant cancers is majorly lacking. In this work, we constructed context specific kinetic models by integrating gene expression data for folate pathway enzymes from untreated, sensitive and resistant variants of 7 cancer cell lines. The kinetic model contained 11 enzymatic reactions, 1 non-enyzmatic reaction and 6 folate metabolites. Differences in pathway activity were identified from steady state fluxes and metabolite concentrations obtained for the different models. Fluxes and metabolite concentrations varied among untreated cancers suggesting inherent differences in pathway activity in different tissue types. More interestingly, the response of sensitive and resistant variants of each cancer type was also seen to vary in simulations of MTX-inhibition. However, accumulation of dihydrofolate at steady state in all sensitive models and decrease towards normal levels for their resistant counterparts was commonly observed for most cancers. Further, metabolic control analysis was used to identify crucial flux controlling enzymes in the folate pathway which can be targeted in MTX-resistant cancers.

ABSTRACT. Observational studies have revealed a complex relationship between fluid overload, ultrafiltration rate (UFR) profile, systolic blood pressure (SBP), and mortality in end-stage renal disease (ESRD) patients. We present a black box model describing the relationship between the UFR profile during a hemodialysis (HD) session and SBP. We developed the model based on 5-year treatment data for approximately 10,000 hemodialyzed patients, provided by Visonex, LLC (Green Bay, WI). The data included patient demographics, outcomes (surviving, deceased, received transplant), and in-session SBP and rate of fluid removal (ultrafiltration rate, UFR) measurements. We modeled the effect of UFR changes on the patient’s SBP by estimating a single input single output (SISO) black box non-linear model (Systems Identification Toolbox, MATLAB). Changes in the extracellular fluid volume, obtained by integrating UFR over time, represented the input signal, whereas the SBP fluctuations about the mean represented the output signal. This model can inform how modifying the UFR profile may favorably influence intradialytic SBP fluctuations, now tested as follows. After determining model parameters using measured data, various UFR profiles were simulated and fed into the model until a desirable output was obtained, under the constraint that the same volume of fluid was removed. The process yielded a new UFR profile that allowed for lower amplitude fluctuations of SBP, thereby reducing circulatory stress. This method has potential therapeutic applications in that it can be used to minimize acute and chronic cardiovascular complications of fluid overload by determining an optimal UFR profile based on the patient’s treatment history, including factors predicting switching from one type of intra-session SBP response to another.

ABSTRACT. A number of devastating human disorders, for example Alzheimer's disease (AD), Hungtington's diseases, type 2 diabetes and transmissible spongiform encephalopathies (TSEs), are associated with the abnormal folding and assembly of proteins. The net result of this misfolding is the formation of large insoluble protein deposits and small toxic and transmissible protein particles in a state called amyloid. What are the molecular mechanisms that govern the amyloid fibrils’ potential to seed the formation of new aggregates, to propagate the amyloid state as prion particles, and to damage cells in amyloid-associated diseases? Here, a generic systems model describing the lifecycle of amyloid is presented in the context of experimental data derived from time dependent fluorescence kinetic assays and atomic force microscopy analysis of amyloid fibrils. We have developed AFM imaging approaches that are capable of resolving the fibril particle concentrations, their length distributions, as well as their toxic and infective potential to cells. With these approaches, we have shown that the disease-associated properties of amyloid can be linked to small nano-sized amyloid particles created through the breakage of amyloid fibrils. The approaches we have developed offer new opportunities to determine, quantify, and predict the course and the consequences in amyloid assembly of cytotoxic, infectious as well as functional amyloid systems.

References: Eugene, Xue, Robert, Doumic-Jauffret, Insights into the variability of nucleated amyloid polymerization by a minimalistic model of stochastic protein assembly, J Chem Phys 14 (2016) 175101 Xue and Radford, An imaging and systems modeling approach to fibril breakage enables prediction of amyloid behavior. Biophys J 105 (2013) 2811-2819. Xue, Homans, Radford, Systematic analysis of nucleation-dependent polymerization reveals new insights into the mechanism of amyloid self-assembly. Proc Natl Acad Sci U S A 105 (2008) 8926-8931

ABSTRACT. Motivation: Processes dynamics in heterogeneous cell populations are encoded in high dimensional single cell snap-shot data. Reconstruction of the embedded low dimensional process manifold is successfully achieved by algorithms like wanderlust or diffusion maps. However, the derived process scales e.g. wanderlust axis or pseudotime are in general not equal to the true scale on which the observed process evolves. Results: We describe a universal and at the same time simple transformation scheme to recover the true process scale from any pseudotime algorithm. The transformation is based on the knowledge of the expected distribution of cells on the process scale which may be temporal, spatial or something else. We applied this method successfully to reconstruct the dynamics of the cell cycle machinery. Furthermore, we demonstrate the power of the method by reconstructing the spatial composition of multicellular tumor spheroids from snapshot data of dissociated cells. By reconstructing spatial spheroid composition from single cell data, we provide insights to the evolution of cell heterogeneity within tumor microenvironments and the associated responsiveness to therapeutics. Outlook: The method may be applied to lineage dynamics in cell differentiation or more sophisticated spatial geometries e.g. intestinal epithelial tissue. By transforming different psoudotime scales to the true process scale, the method can be used to integrate different experiments into a unified scale.

ABSTRACT. In transcriptional regulatory networks (TRNs), a coherent type 1 feed forward loop (C1-FFL) is a structure in which one transcription factor induces the expression of another, and the two induce the expression of a target gene. An AND-gated C1-FFL can suppress responses to spurious and short-lived signals, and adaptation for this purpose is hypothesized to contribute to the abundance of C1-FFLs seen in real TRNs. In order to test this and similar adaptive hypotheses, we develop a null model that can computationally capture how the topology of TRNs is shaped by non-adaptive factors during evolution. Gene duplication and deletion shape TRN topology very differently from the gain and loss of individual regulatory connections; our model captures the ratio between these mutation types in a realistic way. In addition, background stochastic “burstiness” in the expression of other genes might attenuate signals into forms quite different from the idealized transient signal previously considered in the C1-FFL context. We simulated the evolution of TRNs under i) selection for filtering out a brief square-wave spurious signal, ii) non-adaptive evolution to respond to an error-free signal in the presence of background stochasticity in gene expression, and iii) neutral evolution. We find that AND-gated C1-FFLs evolve readily under selection for a filter against the idealized spurious signal, but rarely under the other conditions. We also find that intrinsic noise in gene expression can promote the evolution of OR-gated C1-FFLs, suggesting that OR-gated C1-FFLs adapted for intrinsic noise might also explain the observed overrepresentation of C1-FFLs in real TRNs.

ABSTRACT. A common way to analyze gene expression data is by using pathway networks from databases such as KEGG, along with network visualization and analysis tools like Cytoscape. Applying gene expression data onto a KEGG network highlights changes in how the pathway is regulated, and can help understanding the biology. A major challenge with KEGG networks is that multiple genes are often found to be responsible for the same enzymatic reaction, and Cytoscape only shows the first listed gene in a given position, which can lead to incorrect interpretations. By manipulating the network files, it is possible to show all possible genes for a given reaction as a common box, and display fold change and p-values for each of these genes. However, displaying all possible genes gives a more complicated and messy analysis, which can be hard to interpret. We suggest a solution to simplify both analysis and interpretation, by using the counts from RNA-seq to collapse all genes in a KEGG box into one measurement. We show that the number of counts for genes can be used as a proxy for relative expression level between genes. Assuming that activity for an enzymatic reaction mainly depends upon the gene with the highest number of reads, and by weighting the reads on gene length and gene ratio, a new expression value is calculated for the KEGG box as a whole. The weighted counts for each box are then used for differential expression analysis, and these new up-/down-regulations can be applied to the networks. Using prostate cancer as a model, we integrate RNA-seq data from two prostate cancer patient cohorts with data and prostate cancer metabolism data from literature. We show how using the reduced boxes instead of original genes gives more reliable pathways that are easier to interpret biologically.

ABSTRACT. Introduction: The dynamics of biochemical networks can be modeled by systems of ordinary differential equations. However, these networks are typically large and contain many parameters. A model reduction procedure is presented in [1] where metabolites are reduced one at the time minimizing the difference between the original and reduced model. We extended this model reduction procedure using a new way to compare and select models.

Methods: Instead of considering reductions of only one metabolite at the time, we consider reductions of several metabolites simultaneously for different parameter sets sampled from a distribution. We clustered the dynamics of metabolites using a symmetric dissimilarity measure together with single linkage clustering. For the model consisting of four metabolites shown in the left panel of the figure, we show the full reduction procedure.

Results: We found that the best model reduction depended on the parameter and initial values of the full model. For some values, the iterative method in [1] did not find the best reduced model. For example, in the middle panel of the figure removing two metabolites is better than only removing one. We found that for sampled parameter values, the models with both intermediate metabolites removed form a cluster together with the original models as seen in the right panel of the figure.

Discussion: Our analysis shows that considering all possible model reductions simultaneously can lead to better reduced models since the model reduction algorithm in [1] can be sensitive to parameter and initial values. Cluster analysis enables us to include parameter uncertainty in the model reduction, where the clusters consist of models with similar dynamics.

Acknowledgements: This research was supported by the Research Council of Norway through grant 248840.

References: [1] Shodhan Rao, Arjan van der Schaft, Karen van Eunen, Barbara M. Bakker, and Bayu Jayawardhana. A model reduction method for biochemical reaction networks. BMC Systems Biology, 8(1):52, 2014. URL http://dx.doi.org/10.1186/1752-0509-8-52.

ABSTRACT. Infection, replication and mutation are the key mechanisms governing the population dynamics of viruses and driving their evolution. In particular, RNA viruses have the high mutation rates which enable them to form highly diverse populations within a single host, evade immune responses and develop resistances to drugs. The advent of next generation sequencing platforms opens up the possibility of understanding virus evolution at the molecular level. However, virus evolution is a multi-scale process, involving multiple entities, e.g. heterogeneous host cells and receptors, and intracellular and extracellular mechanisms. Understanding this complexity necessitates the development of a multi-scale model that is computationally challenging and thus requires high-performance computing. We present an exploratory simulation model to study the evolution of heterogeneous virus populations in heterogeneous cell environments. This is a unique model that operates at three scales capturing the core of the evolutionary process. To the best of our knowledge, this is the first HPC-based simulation of its kind. Our model algorithmically represents key known or hypothesized biological mechanisms in a computationally feasible way. This can help biologists to test what-if scenarios and gain insights into the evolutionary dynamics of RNA viruses in a cell culture with less effort and resources in a shorter amount of time. The model does not aim to substitute for lab experiments, but to guide them and to speed up knowledge discovery. We simultaneously simulate the extracellular and the intracellular activities of RNA viruses, explicitly represented as nucleotide (NT) and amino acid (AA) sequences, as well as the diffusion of virus particles among cells. The simulation produces the quasispecies population evolved from initial viral and cell populations as a result of genotype-specific replication. Given NT/AA frequencies available in public databases, our model tentatively scores each genotype’s replication capabilities. This is a novel data-driven method, designed to explore hypotheses.

ABSTRACT. Catalyzing biochemical reactions with high specificity has traditionally been considered a hallmark of metabolic enzymes, which are thought to have evolved from catalytically inefficient and unspecific ancestors. Accumulating evidence, however, suggests that many extant enzymes have retained their ability to act on multiple substrates or catalyze different reactions at physiologically relevant rates. Here, we report our results of testing all known metabolic enzymes in Escherichia coli for catalytic activity towards hundreds of potential substrates by incubating them in a mixture of complex metabolome extracts. Our data reveal that enzyme promiscuity extends far beyond what is currently appreciated, and suggest that promiscuous enzymes played a key role for organisms to acquire new metabolic capabilities when confronted with adaptation to changing environments. We expressed and purified 1,043 known metabolic enzymes from E. coli using the ASKA collection. Each enzyme was separately incubated in a complex metabolome cocktail consisting of combined cellular extracts of Escherichia coli grown in different nutrient combinations supplemented with general enzyme cofactors. Over the time course of 15 minutes, aliquots of each reaction mixture were sampled and quenched. A total of 13,000 samples were analyzed using nontargeted flowinjection time-of-flight mass spectrometry, revealing temporal dynamics of 10,000 detected ions annotated as up to 1,500 metabolites based on accurate mass. Discovered promiscuous reactions were integrated into a genome-scale metabolic model of E. coli to investigate their contribution to the adaptability of metabolism to changing nutrient availability in the environment. For 60% of the analyzed enzymes we detected ions depleting or accumulating over time, corresponding to substrates and products of catalyzed biochemical reactions. Among these active enzymes, 70% affected at least one ion that did not correspond to any of their previously known reactants, suggesting promiscuous activity. By comparing the molecular structures of annotated novel reactants with known reactants, we noted that for some enzymes the structures were highly similar, consistent with the concept of “substrate promiscuity” according to which enzymes can catalyze the same biochemical reaction with closely related substrates. In other cases, novel and previously known reactants were structurally dissimilar, suggesting that these enzymes were able to catalyze mechanistically different reactions (also referred to as “reaction promiscuity”). The newly discovered promiscuous activities of several enzymes were (and are being) validated using enzyme assays and deletion mutants. Analyzing potential roles of this wide-spread promiscuity on metabolism within a genome-scale metabolic model suggest that promiscuous enzymes expand the metabolic capabilities in adapting to novel environments (i.e. accessing novel nutrient sources). We will further discuss the implications of wide-spread enzyme promiscuity for biotechnology and medical research. We present the first global enzyme activity screen, revealing unprecedented insights into the unexpectedly diverse catalytic capabilities of metabolic enzymes.

ABSTRACT. The design of biological oscillators is an intriguing topic in Synthetic Biology. Although there have been many artificially designed biological oscillators, none of them has a period of around 24 hours like a circadian clock oscillator, especially solely based on proteins. The core circadian oscillator of cyanobacteria consists of three proteins, KaiA, KaiB, and KaiC. This circadian oscillator could be functionally reconstituted in vitro with these three proteins, and therefore has been a very important model in circadian rhythm research. KaiA can bind to KaiC and then stimulate its phosphorylation, whereas KaiB antagonizes KaiA’s function leading to the de-phosphorylation of KaiC. Using this protein-based circadian clock oscillator as a model, we aim to find its controlling point and re-design this oscillator. To this end, we combined the tools in bioinformatics, evolutionary biology, protein design, molecular biology, and mathematical modeling to study the interaction between KaiA and KaiC. We found that there exist complicated but critical structural movements during the binding of KaiA and KaiC, and these movements are determinant to the oscillation of the KaiABC system. We further revealed that the KaiA has an asymmetric structural flexibility, which regulates its auto-inhibition and the interaction with KaiC. Based on our findings, it could be possible to redesign this oscillator with more interesting functionalities, which would provide useful insights to the design of protein-based oscillators.

ABSTRACT. The importance of p53 dynamics is evident in its control over differential cellular responses to various stress stimuli. It was previously identified that a bimodal regulation of p53 dynamics modulated by DNA damage strength is crucial for cell fate control. Further single cell analysis revealed significant cell-type variation in this bimodal switch that correlated with the dose-dependent DNA damage response. While mammalian cell lines activated similar periodic pulsing of p53 followed by cell-cycle arrest at low DNA damage, they switched to distinct dynamics at high damage, i.e. monotonic increase or a single pulse of p53. Cell lines with monotonically increasing p53 underwent faster and more extensive damage-induced apoptosis than cell lines exhibiting a single p53 pulse. By combining single cell imaging with computational modeling, we uncovered a regulatory module of coupled feed-forward and feedback, involving ATM, p53, Mdm2 and Wip1, which differentially activates bimodal p53 dynamics, with expression of ATM key to determine the collective control strength. We also identified an effective strategy for re-sensitizing the apoptotic response in the cell lines resistant to DNA damage signaling by combinatorial inhibition of Mdm2 and Wip1. Our work results not only elucidate a cell-type dependent dynamic control of p53, but also point to p53 pulsing as a suppressive mechanism that renders resistance to DNA damage signaling.

ABSTRACT. Pro-inflammatory cytokines such as Tumor Necrosis Factor alpha (TNF) have been implicated in the pathogenesis of skeletal muscle atrophy, a phenomenon commonly associated with many chronic systemic inflammatory diseases such as cancer and AIDS. Although the regulation of mRNA levels in response to TNF signaling has been studied extensively in the context of immune activation, such as in macrophages, the mechanisms that control the gene repression program, which is of relevance in muscle cells, are still poorly understood. Here, we examine quantitatively the genome-wide gene expression effects of TNF exposure, both in the short and long term. We characterize the regulatory strategies of dynamic gene induction and repression programs in skeletal muscle cells by measuring both mRNA synthesis and decay rates, and connecting these via mathematical modeling. Our data points to a dominant role of synthesis control in the regulation of both gene induction and repression in response to TNF, but that mRNA half-life control determines the majority of temporal profiles of gene expression. Furthermore, our analysis unveils an unexpected gene expression strategies. We found a cluster of extracellular matrix genes whose sustained TNF-triggered downregulation is actually preceded by transient production overshoot. This suggests that beneficial and detrimental effects of TNF are linked but distinguished by a signaling dynamics. Our fine-grained data highlights the importance of signaling dynamics in mediating TNF effect on muscle cells, and critical interplay between synthesis and degradation control in shaping dynamic gene expression programs.

ABSTRACT. Phages and their host bacteria form ecosystems, which typically involve intriguing spatiotemporal dynamics. Previous studies on spatiotemporal dynamics of phage-bacteria systems mostly focused on growth of circular phage plaques on the bacterial lawn. Recent experiments demonstrated intriguing asymmetrical patterns when phages were inoculated at the edge of an expanding bacterial colony on the agar medium that allows bacterial motility. The bacteria-clear zone grows into sectorial shapes which could persist straightly, flare out, or close up, depending on the infection efficiency. We developed a model to address the pattern dynamics of the phage-bacteria system. We found that the observed sectorial patterns critically depend on negative impacts of the bacterial density on (1) bacterial motility, and (2) phage reproduction. With appropriate parameters, these two relationships cause “freezing” of bacteria-phage plaque boundaries as the bacterial colony expands radially outward, and lead to the sectorial pattern. Our model reveals how spatial niche partitioning emerges in a host-pathogen system, in which the pathogens rely on host, not only to reproduce, but also to spread spatially.

ABSTRACT. Medium chain fatty alcohols (C8-C12) are commonly used in surfactants, detergents, biofuels, and cosmetics. Sustainable production of medium chain fatty alcohols using renewable feedstock such as xylose has been recently achieved by expressing a xylose utilization pathway and a peroxisome-targeted fatty acyl-CoA reductase in yeast Saccharomyces cerevisiae. It was found that much higher yield of fatty alcohol was achieved by using xylose as feedstock rather than glucose. However, the reprogramming of metabolic fluxes of Saccharomyces cerevisiae under xylose-based fatty alcohol production remains largely unknown. To this end, we applied a systems biology approach, namely 13C metabolic flux analysis, to rigorously quantify the differences of carbon fluxes in central metabolic pathways of Saccharomyces cerevisiae between glucose-based fatty alcohol production and xylose-based fatty alcohol production. We found that metabolic flux of the pentose phosphate pathway dramatically increased in xylose-based fatty alcohol production while metabolic flux towards the byproduct ethanol significantly decreased, which accounted for the improved yield of fatty alcohol under xylose fermentation. In addition, the cell maintenance energy was recognized as a rate-limiting step of xylose-based fatty alcohol production. Based on these discoveries, we next optimized the bioprocesses and successfully improved the titer of fatty alcohol to 0.6 g/L, which represents the highest production of medium chain fatty alcohols from xylose.

ABSTRACT. Transcription factors usually function as co-regulators to synergistically induce or inhibit expression of their targets genes. In recent years, genome-scale gene regulatory networks have been generated for multiple organisms including both human and Arabidopsis. Existing module-finding algorithms fail to capture transcription co-regulators in these large-scale networks because these algorithms usually search for groups of densely connected genes (nodes) rather than co-regulating genes. In this study, we developed a new computational tool, CoReg, to identify transcription co-regulators in large-scale gene regulatory networks. CoReg groups genes based on similarities of shared targets between regulators in a network. We applied hierarchical clustering followed by dynamic tree cut to identify co-regulatory modules. We tested our approach in Arabidopsis thaliana (A. thaliana), Escherichia coli (E. coli) and Homo sapiens (H. sapiens) gene regulatory networks. Using these network data sets, we explored the performance of different similarity indices and compared them to existing module-finding algorithms (Walk Trap, Edge Betweenness and Label Propagation). We conducted network-rewiring simulations and found that CoReg+jaccard similarity index performed better than other methods in identifying true co-regulators from rewired networks. Furthermore, we integrated a cell type-specific gene expression data set for Arabidopsis root and a large-scale transcription network of Arabidopsis generated by DAP-seq. We applied CoReg to this integrated data set to identify cell type-specific co-regulators. Our study provides a new tool for dissecting the architecture of gene regulatory networks.

ABSTRACT. A key feature of ageing, in both multi- and unicellular organisms, is the accumulation of damage such as malfunctioning proteins. This evolutionary well-conserved feature can be studied in the baker’s yeast Saccharomyces cerevisiae where one of the main aspects of damage accumulation is the distribution of damage between the mother and daughter cells after cell division. This distribution is partially characterised by the fact that the mother cell retains damage, and thus prevents it from leaking over to the daughter, in order to generate a young and healthy offspring. In this work we aim to elucidate the advantage of dynamic damage retention on a population level by using large scale simulations of a dynamical model representing the whole pedigree. The model consists of ordinary differential equations (ODEs) describing the formation of intact and damaged proteins for each cell in the population coupled with two discrete events corresponding to cell division and death respectively. After division, a new cell represented with a specific set of ODEs is introduced and when cell death occurs the corresponding set of ODEs are removed from the model. The expanding system of ODEs is reinitialised by using the damage retention parameter and the various so-called ageing strategies constitute different means by which cells can alter their retention capacity. By using this novel large scale setup in combination with efficient numerical simulations, we investigate the most optimal ageing strategy in terms of replicative potential, and viable population size. Our results suggest that dynamically changing the retention throughout the lifetime of an individual cell is the preferred division and damage segregation strategy giving cells a direct selective advantage.

ABSTRACT. Interactions of cancer cells with their microenvironment play a critical role in their survival, differentiation, and progression. The cell microenvironment composed of the extracellular matrix (ECM), tumor interstitium and surrounding cells influences cell phenotype through a concoction of physical, mechanical and biochemical factors. Association of growth factors (GF) with the extracellular matrix (ECM) in vivo enhances the duration and potency of GF signaling. Contributions of the 3D solid-phase microenvironment are not reflected in conventional in vitro cell culture techniques. The present study sought to develop a 3D solid-phase microenvironment using heparin-based hydrogel with epidermal growth factor (EGF) for slow, controlled release of the GF and to compare the effects of soluble and solid-phase EGF on a model breast cancer cell line. To reflect the in vivo tumor characteristics, we developed a 3D microenvironment that promotes self-assembly of cells intro spheroids, demonstrated using breast cancer cell line MDA-MB-231. RNA-seq analyses of the molecular signatures of breast cancer can be linked to the phenotype quantified in our hydrogel platform with systems biology and analysis of the changes in parameters that lead to phenotype changes in the cancer cells. We aim to optimize the 3D tumor microenvironment to culture biopsies and isolated CTCs in the clinic for further studies on underlying molecular mechanisms driving site-specific metastasis, and other diagnostic and prognostic markers and therapeutic screening. Computational modeling approaches using our experimental platform may be used to provide further insight into tumor dynamics influenced by the microenvironment.

ABSTRACT. During sepsis, the current leading cause of death in hospitals, neutrophils migrate and accumulate in healthy organs instead of migrating toward the bacterial infection. This dysregulation of neutrophil migratory phenotype as seen in loss of directionality and oscillatory migration has been reported by us in human burn patients with sepsis. The goal of the proposed work is to investigate the heterogeneity in neutrophil migration phenotypes during sepsis. We have developed a microfluidic device to measure neutrophil chemotaxis in an opposing chemoattractant gradient to quantify decision-making. We use two chemoattractants: a pro-resolution and pro-inflammatory chemoattractant to model how a cell makes a decision toward a bacterial infection or an inflammatory signal. Our hypothesis is that low-level pro-inflammatory ‘programming or training’ signals, such as lipopolysaccharide, have a central role in determining the final neutrophil phenotype and in the development of sepsis. Despite tremendous advances in the understanding of signaling molecules and pathways within neutrophils, our understanding of the directional decision-making process is limited, and consequently, our abilities to modulate the activity of neutrophils restricted. Using an ODE-based dynamical framework, we model the interaction of the mutually inhibitory GRK2 and GRK5 proteins and its role in decision-making. Our model results show a bimodal switch between high and low levels of GRK2, indicating that GRK2 may play a role in neutrophil decision-making process. We show that unstimulated cells preferentially migrate toward a pro-resolution signal over a pro-inflammatory signal. However, when primed with a low-dose of a pro-inflammatory mediator, higher ratios of cells migrated toward the pro-inflammatory chemoattractant with a higher velocity. In the future, we will extract cells based on their migration decisions and measure receptor levels to determine the underlying molecular mechanism that drives neutrophil decision-making to validate our computational model. Using this experimental data, we will derive molecular parameters for our model.

ABSTRACT. The surge in antimicrobial resistance requires urgent development of innovative approaches to address the numerous bacterial pathogen threats outlined by the CDC and WHO. Notably, a growing body of evidence suggests that the presumed fitness disadvantages of resistant pathogens conferred by expression of resistance genes is not fully accurate. Compounding this issue, the metabolic responses of pathogens to antibiotic stress surprisingly remain poorly understood despite our great appreciation of specific drug-target interactions. Numerous omics-driven studies focused on treatment of diverse bacterial species with antibiotics have noted a clear shift in metabolism without deeper computational examination. Arising from these data is the increasingly attractive hypothesis that modification of metabolism is a key component of cell death of susceptible bacteria after antibiotic treatment. Further exploration of the relationship between metabolism, antibiotic stress, and resistance is clearly needed.

To address these gaps in our fundamental understanding, we have compared the metabolic behaviors of wild type and resistant strains Escherichia coli. Specifically, we have performed RNAseq to measure mRNA abundance wild-type, ampicillin resistant, and kanamycin resistant E. coli. Transcriptomic data were integrated with the genome-scale metabolic model for E. coli to generate metabolic flux predictions for the stressed and unstressed states. To compare these unique phenotypes, we are developing a computational pipeline that combines the strategies of numerous algorithms, yielding robust and testable predictions for differential states. Our preliminary findings highlight that pathogenic bacteria may reductively constrain their metabolism upon antimicrobial challenge. To validate and improve these predictions, we have recently initialized ¹³C fluxomics experiments, which will refine our view of the central carbon metabolism behaviors of these pathogens.

ABSTRACT. The asymmetric cell division cycle in Caulobacter crescentus is controlled by a cascade of cellular processes from DNA replication and segregation to cytokinesis and cell division to ensure the progression of cell growth and production. The timing of these cellular functions is regulated by the biochemical and genetic logic circuitry which initiate precise temporal activation of a host of proteins. Though a full map of how the cell directs and coordinates diverse mechanisms and cell cycle modules is still unclear, the genetic circuit alone predominately controls cell cycle regulation.

With recent discoveries and more information such as temporal dynamics for mRNAs and spatial distribution for proteins made available, we established an elaborate gene-protein regulatory network to achieve the oscillating cell-type gene expression coordination, centering on five master regulatory proteins: CcrM, DnaA, GcrA, CtrA and SciP, which directly controls approximately 60% of the cell genes. We explicitly characterized the control mechanisms for three key aspects of the genetic circuit: DNA replication, the function of DNA Methyltransferases CcrM, and the interactions between the five regulators themselves. Particularly, we have introduced mRNA variables which play a role in the successive translation of the five regulatory proteins over the course of the cell cycle. The genetic circuit network is cast into a deterministic computational model to explore how it involves cell cycle progression and cell differentiation as an integrated system. 

Simulated results successfully reproduce the genetic circuit with an oscillating pattern and fits well with experimental data of both the proteins’ and the corresponding mRNAs’ temporal behavior. We hope to incorporate the genetic model with other separate cellular modules to build the whole picture of Caulobacter Cell Cycle in the future. A similar methodology applied to other bacterial genera may be useful to fields such as biomedical sciences and natural resources.

ABSTRACT. Capacity of plants of converting light to different valuable products makes them highly relevant to “green factories”. To reduce burden of arable land usage for valuable product production there is need to increase significant plant productivity. In this study to we use Arabidopsis thaliana genome scale model in a synthetic biology based approach to increase production of z-abienol and isoprene depending on light intensity and diurnal cycle. Genome scale model is primarily derived from annotation aracyc 13.0 in the Aracyc database. We are using linear programming to demonstrate model capability to produce biomass components (amino acids, nucleotides, lipids, starch, cellulose). Model structure uses cytosol, mitochondria, plastid, vacuole and peroxisome compartments. Compartmented reactions are mostly used for interaction analysis of light and diurnal cycle involved reactions. Kinetic model of MEP as plant isoprenyl precursor pathway located in plastids will be used in combination with genome scale model to assess the feasibility of kinetic model steady states suggested by optimizations.

ABSTRACT. Quantitative modeling of complex biological systems requires a large number of parameters such as kinetic constants and experimental conditions. Most of the data is buried in a constantly and quickly growing number of scientific publications, thus not being available for computation. It is practically impossible for researchers to keep up with this data flood without support of computers and databases.

Our group maintains a database for the storage of such information namely SABIO-RK, a resource for biochemical reactions and their kinetic properties (http://sabiork.h-its.org/).

SABIO-RK is part of the data management node NBI-SysBio within de.NBI (German Network of Bioinformatics Infrastructure, http://www.denbi.de/) program which is a newly-established BMBF-funded initiative having the mission to provide comprehensive first-class bioinformatics services to users in life sciences. In this context our group provides a couple of services including the manual extraction of kinetics data from the literature upon user requests. The information extracted comprises the origin of the data (publication), kinetic parameters, kinetic formula and kinetic type, reaction, organism, tissue, strain, compartment, enzyme classification and experimental conditions. Extracted data are provided in a structured way either accessible via our web interface, or as spreadsheets, or as SBML files for direct import into modelling and simulation tools.

Currently, this service is still free of charge and not restricted to German scientists. Any requests can be made using the corresponding form (http://sabiork.h-its.org/contactFormSabio) or by email (sabiork@h-its.org).

ABSTRACT. In spite of increasing evidence showing that cellular drug transport is mainly, if not exclusively, carrier-mediated rather than via passive diffusion [1], most chemical compounds still lack an associated transporter that explains their entry and distribution in cells and tissues. Currently known drug transporters correspond to two main protein families [2]: the ATP-Binding Cassette (ABC) transporter family, whose members are often involved in xenobiotic efflux and drug resistance, and the large and heterogeneous family of SoLute Carriers (SLCs), which includes various cases of drug uptake. We recently argued that SLCs are a highly neglected gene group [3], with most of its members still poorly characterized, and thus likely to include many yet-to-be-discovered cases of drug transport.

We therefore mined a publicly available pharmacogenomics dataset involving 1001 molecularly annotated cancer cell lines and their response to 256 anti-cancer compounds [4] in order to prioritize new SLC/ABC-drug associations. To this end, regularized linear regression models (Elastic Net) were generated to predict drug response based on SLC and ABC data (expression levels, Single Nucleotide Variants, Copy Number Variations), and their predictive performance assessed. The best predictive models included some known transporter-drug pairs, such as SLC35F2-YM155 [5] or several cases involving the multidrug resistance protein ABCB1, together with other associations not yet described. We are currently carrying out experimental validation of these by using CRISPR-Cas9-based single genetic knockouts in haploid cells (HAP1) and in a panel of cancer cell lines.

[1] Kell et al., Nat.Rev.Drug.Disc., 2011 [2] Giacomini, Nat.Rev.Drug.Disc., 2010 [3] César-Razquin, Snijder, et al., Cell, 2015 [4] Iorio et al., Cell, 2016 [5] Winter et al, Nat.Chem.Biol., 2014

ABSTRACT. In response to a concentration gradient of chemo-attractant, {\sl E. coli} bacterium modulates the rotational bias of flagellar motors that control its run-and-tumble motion, to migrate towards regions of high chemo-attractant concentration. Presence of stochastic noise in the biochemical pathway of the cell has important consequence on the switching mechanism of motor bias, which in turn affects the runs and tumbles of the cell in a significant way. We model the intra-cellular reaction network in terms of coupled time-evolution of three stochastic variables, kinase activity, methylation level and CheY-P protein level, and study the effect of methylation noise on the chemotactic performance of the cell. A good performance consists of reaching the favorable region quickly and localizing there in the long time limit. Our simulations show that the best performance is obtained at an optimal noise strength. While it is expected that chemotaxis will be weaker for very large noise, it is counter-intuitive that the performance worsens even when noise level falls below a certain value. We explain this striking result by detailed analysis of CheY-P protein level statistics for different noise strengths. We show that when the CheY-P level falls below a certain (noise- dependent) threshold, the cell tends to move down the concentration gradient of the nutrient, which impairs its chemotactic response. This threshold value decreases as noise is increased, and this effect is responsible for noise-induced enhancement of chemotactic performance. In a harsh chemical environment, when the amount of nutrient depletes with time, the amount of nutrient intercepted by the cell trajectory, is an effective performance criterion. In this case also, depending on the nutrient lifetime, we find an optimum noise strength when the performance is at its best.

ABSTRACT. Estrogen receptor positive (ER+) breast cancer is the most common type of breast cancer today. The major treatment of ER+ cancer is endocrine therapy that targets the ER signaling or estrogen production, using selective estrogen-receptor response modulators to block the effect of estrogen, aromatase inhibitors to stop the production of estrogen, or estrogen-receptor down-regulators to reduce the number of ERs. Although these therapies can improve overall survival, breast cancers often develop resistance to endocrine therapies, and this resistance represents a significant impediment to successful treatment. The unfolded protein response (UPR) pathway, which is elevated in endocrine-resistant cells, plays a role in mediating endocrine resistance by influencing the balance between apoptosis and autophagy. To better understand the role of the UPR in the development of endocrine resistance, we built an ordinary differential equation model of the UPR for three cell lines: MCF7 (endocrine therapy sensitive, estrogen dependent), LCC1 (endocrine therapy responsive, estrogen independent) and LCC9 (endocrine therapy resistant, estrogen independent). The UPR of these three cell lines is induced by the strong reducing agent, dithiothreitol (DTT), which blocks disulfide-bond formation and leads to endoplasmic reticulum stress. DTT stress is modeled as an increased flow of unfolded proteins into the endoplasmic reticulum. The mathematical model is constructed from quantitative measurements of key UPR signaling molecules XBP1, GRP78/BiP, eIF2α and phosphorylated eIF2α. The models of the three cell lines have the same structure, and a minimal change of parameter values in the LCC1 model can adequately recapitulate the data from MCF7 and LCC9 cell lines. This result leads to the hypothesis that UPR signaling is the same for these three ER+ cell lines, and that the acquisition of resistance is due to changes in the balance of UPR signaling. The model also can contribute to identifying new drug targets for reversing endocrine resistance.

ABSTRACT. Hepatitis B is a major disease that chronically infects millions of people in the world, especially in developing countries. Currently, one of the effective vaccines to prevent Hepatitis B is the Hepatitis B Small Antigen (HBsAg), which is mainly produced by the recombinant yeast Saccharomyces cerevisiae. In order to bring down the price, which is still too high for people in developing countries to afford, it is important to understand key cellular processes that limit protein expression. In this study, we took advantage of yeast knockout collection and systematically screened 194 S. cerevisiae strains with single genes knocked out in four major steps of the protein secretory pathway, i.e., endoplasmic-reticulum (ER)-associated protein degradation, protein folding, unfolded protein response and translocation and exocytosis. The screening showed that the single deletion of YPT32, SBH1, and HSP42 led to the most significant increase of HBsAg expression over the wild type while the deletion of IRE1 led to a profound decrease of HBsAg expression. The synergistic effects of gene knockout and gene overexpression were next tested. We found that simultaneously deleting YPT32 and overexpressing IRE1 led to a 2.12-fold increase in HBsAg expression over the wild type strain. The results of this study revealed novel genetic targets of protein secretory pathways that could potentially improve the manufacturing of broad scope vaccines in a cost-effective way using recombinant S. cerevisiae.

ABSTRACT. Glioblastomas are highly diffuse tumors with very poor prognosis. Their critical anatomical location restricts therapeutic options, reflected in tumor recurrence typically within months after primary resection. How glioblastomas evolve under treatment and which role tumor heterogeneity plays in this process is of major interest in the search for new therapeutic approaches. To reconstruct glioblastoma evolution from tumor samples, we analyzed 22 pairs of matched primary and relapse tumors with whole genome sequencing and investigated their degree of genomic heterogeneity. This was achieved by modeling the sequencing read count data as a sampling result from a mixture of genetically distinct subclones which correspond to the tips of a phylogenetic tree. Our analysis, restricted to a finite set of candidate trees, suggests that glioblastomas are heterogeneous tumors with two to three dominating subclones per tumor sample. The subclonal distributions support a model of branched evolution in which a common stem of mutations is shared by the entire tumor before additional cycles of mutation and selection shape genetically distinct subpopulations. Our data further indicates that relapse tumors do not directly evolve from subclones identified in the primary sample, but from common ancestor populations which are not recovered in the primary sample. Our findings shed light on the clonal evolution of glioblastomas, both between tumor initiation and primary lesion, and between primary lesion and recurrence. A thorough analysis of the mutations assigned to the tumor stem is ongoing and will help us to distinguish early from late mutations and thus gain insight into the early development of glioblastoma.

ABSTRACT. Quantitative characterizations of horizontal gene transfer mechanisms are needed for understanding and predicting the dynamics of gene distribution in natural and engineered systems. In this study, we developed a mathematical characterization of plasmid conjugation between two bacterial populations (filter mating).

We mated two E.coli strains. The donors harboured the self conjugative, GFP-coding plasmid PKJK10. The recipients expressed RFP from the plasmid PSB1C3. Time series assays were made by flow cytometry to quantify the distribution of the three subpopulations involved in the filter mating process (GFP+/RFP- donors, GFP-/RFP+ recipients, and GFP+/RFP+ transconjugants). Corresponding measures of optical density determined the temporal variation in the abundance of each population.

We used the data to fit ordinary differential equation models of the process, based on previously published model frameworks. Model comparison tools were applied to arrive at an optimal model formulation, and the accuracy of the best-fit parameter estimates was assessed via uncertainty analysis. We tested the model’s predictive power by comparing model simulation to experimental results that demanded extrapolation from the training data. These comparisons provide evidence that the model can be successfully used as a predictive tool for characterizing horizontal gene transfer mechanisms in natural or synthetic systems. This work recently appeared in Frontiers in Microbiology.

ABSTRACT. Molecular biology is experiencing a revolution, in one part thanks to new technologies to measure and perturb biological systems in vitro, and also due to the growing importance of mathematical modeling which enables us to understand biological mechanisms in a more profound way. However, one crucial point in this transforming field is the need to provide completely new tools, which should be computationally efficient, versatile, and compatible.

To this end, we developed SigNetSim, a web platform written in Python and using the Django framework. SigNetSim uses Bootstrap as a graphical front-end, which makes it usable on most devices. It is designed to be installed on computation servers, with all the work being executed server-side.

Users can create and edit biological models in the standard SBML format. Additionally, the platform supports the SBML comp package, which allows the user to write hierarchical models, where models contains other SBML models as sub-models. This allows reuse of models, and simplify the writing of large models. SigNetSim performs simulation both for time series and steady states, and plots the results using JavaScript interactive libraries. It is partially compatible with SED-ML format, which stores simulation settings and allows users to easily reproduce simulations from literature. Moreover, SigNetSim can use combine archive format, which stores both the SBML model and the SED-ML simulation file into one file. The platform includes a simple database to store experimental data, which can be used to simulate models according to a set of initial conditions, or to fit models using a parallelized simulated annealing algorithm. Users can perform various types of dynamical analysis, including bifurcation analysis. SigNetSim also supports model annotation, following the MIRIAM guidelines.

Finally, SigNetSim is available on GitHub under GPLv3 license. As a case study, we present a model for mitogenic pathways in mouse Y1 tumor cells.

ABSTRACT. Colorectal cancer is becoming a major threat these days whereas only a part of patients respond to targeted anti-cancer therapies. Such heterogeneous responses of colorectal cancer patients are primarily due to the complex molecular interactions of cancer cells that are differentially wired according to the mutational profiles. Analysis of such diverse cancer networks is therefore required to develop precision medicine. In this study, we have reconstructed a colorectal cancer signaling network including four frequently mutated canonical pathways: mitogen-activated protein kinase, DNA damage, Wnt, and transforming growth factor beta signalling pathways. Based on this, we have developed colorectal cancer cell-type specific network models by employing discrete logic-based Boolean network modeling. We revealed the distinct dynamical characteristics of various cell-types by investigating the attractor landscape of each network model. Finally, we have identified an optimal drug combination for each network model by evaluating the drug responses on the basis of efficacy, toxicity, and potency. Our study provides a new insight into the optimized therapy for colorectal cancer.

ABSTRACT. T cells scan body cells for pathogen-derived peptides that are presented by major histocompatibility complex (pMHC) molecules. To discriminate between self and foreign peptides, cytotoxic T cells express the T cell receptor (TCR) and the co-receptor CD8, both binding to pMHC. Several models of this process have been proposed, however, lacking molecular detail. Here, we combine dose response measurements for soluble pMHC tetramers of graded affinity to the TCR in the absence or presence of CD8 with mathematical modeling. Data simulation selected one particular mechanistic model of the TCR-pMHC-CD8 interactions. This model shows that CD8 enhances binding of pMHC to T cells. However, at the same time CD8 reduces the capability of low affinity pMHC to stimulate TCRs by hindering pMHC from binding to the TCR. In sharp contrast, CD8 increases the capability of high affinity pMHC to stimulate TCRs by allowing few pMHC to simultaneously bind to TCR and CD8 thereby recruiting the CD8-associated kinase Lck to those TCRs. T cell activation data confirm this prediction, leading us to describe a detailed molecular mechanism in which CD8 enhances both the specificity and sensitivity of affinity-based antigen discrimination by T cells.

ABSTRACT. Constraint-based metabolic models have become fundamental and trusted tools in systems biology. Several layers of biological information are combined in a compact format in order to describe a metabolic model. A richly annotated model is required for its various areas of application and represents a veritable knowledge base about an organism's metabolism. However, coherently describing a complex interlinked system such as metabolism is a challenge in and of itself that is only aggravated by the current lack of cohesive, widely-accepted, testable, and modern standards.

Here, we introduce memote (Metabolic Model Tests {https://github.com/biosustain/memote}), a Python package designed to run a given model through a set of hard and soft tests and generate a report that reflects model integrity. Soft tests focus on aspects that do not influence the performance of the model, such as syntactic conventions whereas hard tests determine whether a model is fully functional.

While memote can be run locally as a stand-alone testing suite, it shows its full potential when combined with web-based version controlling (Github) and continuous integration tools (Travis CI). Every tracked edit of a model automatically triggers the memote test suite, and generates a corresponding report that facilitates factual debate of model changes.

Thus, memote not only allows researchers to more quickly iterate through the design-build-test cycle but also provides the scientific community with a measure of quality that is consistent across setups, as well as an opportunity to interact and collaborate by establishing workflows for publicly hosted and version controlled models.

ABSTRACT. Introduction: Bisphenol A (BPA) is a chemical used in the production of polycarbonate plastics and is notable for its endocrine-disrupting effects acting as a xenoestrogen. Its ubiquitous nature in the environment is highlighted by the fact that 92.6 % of adults excrete BPA. Recently estrogens have been implicated as potential agents in the development and progression of prostate cancer. miRNAs act as gatekeepers in transcription modules, increasing the robustness of transcription networks. Objectives: We hypothesize that BPA exposure negatively impacts transcriptional programs via alterations in miRNA expression, resulting in less robust biological circuits that are more prone to unstable outputs that can ultimately lead to prostate cancer. Materials and Methods: We extracted RNA from Human Prostate Epithelial Cells (CloneticsTM), derived from a 23-year-old male, and cultured in the presence or absence of two BPA doses (5 and 25 nM) and ethinylestradiol EE2 (0.1 nM) for 24 hours. RNAseq was performed on an Illumina HiSeq2500. DE mRNAs and miRNAs were determined using Limma. Correlation and systems analyses of DE mRNA and miRNAs were performed using multiMir and iPathwayGuide. Results and Discussion: With both BPA exposure conditions, hsa-miR-335-3p, hsa-miR-543, hsa-miR-424-5p, hsa-miR-548h-5p, hsa-miR-493-5p were the top 5 miRNAs when ranked by the number of affected mRNA targets. EE2 exhibited similar results, with the exception of hsa-miR-493-5p, which is ranked 6th, and hsa-miR-548d-5p, which is ranked third and is only up-regulated with EE2. All miRNAs noted impact cell cycle at a systems level. hsa-miR-335-3p is the top ranked miRNA across all three conditions and is linked to obesity and ER-α inhibition, well-known effects of BPA exposure Conclusion: BPA exposure, even at low doses, disrupts miRNAs involved in cell cycle regulation. This disruption suggests a loss of robustness in this system, which can potentially facilitate carcinogenesis.

ABSTRACT. Introduction: Exposures of Atlantic cod (Gadus morhua) to environmental contaminants such as PCB153, increase the levels of enzymes involved in fatty acid synthesis in cod hepatocytes. Mathematical modeling using non-linear differential equations may help to identify how this affects the synthesis rate of fatty acids and triglyceride stores in cod hepatocytes. Such models can also describe the dynamics of fatty acid synthesis and the flux at different boundary conditions. In this study, we have used mathematical modeling to investigate the impact of feed-forward and feed-back regulatory interactions in the pathway as well as possible time delay for the reaction catalyzed by the multi-enzyme complex fatty acid synthase (FAS).

Method: Using mass action and Hill type kinetics, we constructed a model involving non-linear controlled differential equations. A synthesis pathway from citrate to triglycerides is assumed, including also palmitoyl-CoA and the mitochondrial import for subsequent beta-oxidation. In this pathway, palmitoyl-CoA is a feedback inhibitor of acetyl-CoA carboxylase and citrate is an allosteric activator. A non-competitive inhibition with regard to acetyl-CoA is assumed. Insulin and glucagon triggers are assumed as control variables to consider the effects of external disturbances on the system. The pathway with regulatory interactions is shown by the diagram.

Preliminary result: By non-dimensionalizing equations the number of parameters was reduced, and the steady state solutions were found analytically. Constraining the concentrations to be non-negative gave us valid parameters region. The stability of the steady state solutions was studied by investigation of eigenvalues of the Jacobian matrix of system of equations looking for bifurcation points.

Further work: We are studying how different choices of the control variables can affect the trajectory of the system and the stability of the steady state solutions.

Acknowledgements: This research was supported by the Research Council of Norway through grant 248840.

ABSTRACT. Recently, a synthetic circuit in E. coli demonstrated that two proteins engineered with LAA tags targeted to the native protease ClpXP are susceptible to crosstalk due to competition for degradation between proteins. To understand proteolytic crosstalk beyond the single protease regime, we investigated in E. coli a set of synthetic circuits designed to probe the dynamics of existing and novel degradation tags fused to fluorescent proteins. These circuits were tested using both microplate reader and single-cell assays. We first quantified the degradation rates of each tag in isolation. We then tested if there was crosstalk between two distinguishable fluorescent proteins engineered with identical or different degradation tags. We demonstrated that proteolytic crosstalk was indeed not limited to the LAA degradation tag, but was also apparent between other diverse tags, supporting the complexity of the E. coli protein degradation system.

ABSTRACT. Constraints-based reconstruction and analysis (COBRA) is widely used to interpret and predict the interplay between genotype and metabolic fluxes. The community-developed cobrapy Python package implements functionality to read, write, edit, and adjust COBRA models and to perform simulations using numerous popular algorithms. Since the first releases in 2012, the cobrapy project has gained considerable attention thanks to its broad feature-set with extensive documentation, and to being free/open source software without any non-free dependencies. The core classes of cobrapy form a basic infrastructure for constraints-based modeling that is easy to reuse in other packages, facilitating the development of new functionality as well as increasing potential interoperability between packages. In order to simplify the implementation of new algorithms, we have drawn on our experiences with the early versions of cobrapy, the development of the strain design package cameo, and the mathematical modeling package optlang, to enhance the cobrapy core classes. Interaction with the software that actually solves the mathematical problem in the COBRA model is now provided by optlang, which greatly facilitates the implementation of new COBRA algorithms and encourages the contribution of new algorithms from the research community. Several simulation algorithms have already been refactored for increased readability and performance. Here, we present the new functionality and outline the way forward for the role of cobrapy as a freely available infrastructure package for efficient constraints-based modeling of metabolic flux in python.

ABSTRACT. Cells undergo an evolutionary process such that cooperative cancerous characteristics are acquired along with the accumulation of somatic mutations during tumorigenesis. There is still a lack of systems biological understanding as to how the cancerous phenotypic characteristics change and how they are synergistically connected in cancer evolution processes. In our previous study, using a genome-wide analysis of the somatic mutations in colorectal cancer patients on the basis of a large-scale molecular interaction network, we found that a giant cluster of mutation-influencing subnetworks undergoes a percolation transition during colorectal tumorigenesis and ultimately results in a giant percolated cluster (GPC) that includes a set of genes closely related to the cancer development. In this study, to further investigate the cooperation of cancerous phenotypic characteristics in tumorigenesis, we projected the genes obtained from GPC onto a phenotypic network composed of 50 hallmark gene sets from MSigDB, to be named a Hallmark Gene Set Network (HGSN). By clustering patient-specific HGSNs according to the occurrence of driver mutations (e.g., APC, KRAS, PIK3CA, SMAD4, and TP53), we found that the HGSN becomes dense as driver mutations are accumulated in a defined order. When we categorized the hallmark gene sets according to the Weinberg’s cancer hallmarks, we further found that a more malignant cancer has a higher density of the subnetwork related to the metastasis category. Taken together, we conclude that cancerous phenotypic characteristics are connected with each other cooperatively during colorectal tumorigenesis and that HGSN enables us to infer optimal anti-cancer targets that can prevent such critical connections among the cancer-related hallmark gene sets.

ABSTRACT. Bullous Pemphigoid (BP) is the most common autoimmune blistering disease in Europe. As both the incidence of the disease and the relative proportion of the elderly population continue to rise, it represents a significant medical burden. Whereas some progress has been achieved in defining genetic risk factors for autoimmune blistering diseases, no environmental agent has been conclusively identified. Emerging evidence suggests that host immunity may influence the skin microbiota while the latter modulates cutaneous immunity. Nevertheless, the relationship between skin microbial communities and autoimmune bullous disease has yet to be studied in humans. Here, we aim to characterize and compare the skin microbiome of BP patients and healthy, age-matched controls at numerous body sites. Similar to what has been shown in healthy controls, the composition of skin microbiota in BP patients appears to be very divergent and site-specific. Microbial phylum abundances differ between perilesional sites of BP patients and the same anatomic locations of control patients. A distinct cutaneous microbiota profile, which correlates with BP, further strengthens the significance of commensal-host interaction on our immune system. Moreover, these results raise the possibility that the cutaneous microbiome may contribute to the pathogenesis of BP, with important implications for the treatment of this disease.

ABSTRACT. Owing to the development of high-throughput measurement technologies and the advancement of systems biology for network inference methods, the network models for molecular regulatory interactions are rapidly increasing and becoming more available. However, network analysis methods are still limited and often require very accurate nonlinear dynamic models. In particular, to predict the effects of any perturbation in the network components such as nodes or links, we need a detailed nonlinear mathematical model based on kinetic parameter estimation with time-series measurements. Such time-series measurements for parameter estimation limits constructing a large-scale dynamic network model. Hence, there is a pressing need to develop a new network analysis method with which we can predict the perturbation effects based only on network topology. In this study, we present such a novel method with a linear system approximation. Intriguingly, we could predict the perturbation effects with an accuracy of more than 60% based on network topology and mutation profiles of various cell types. The proposed method can be used to predict differential drug responses depending on cell types and also to identify new promising drug targets.

ABSTRACT. Predicting anticancer drug responses of cancer patients is a crucial problem to achieve a higher therapeutic efficacy and to implement precision medicine. It is becoming more feasible to construct a computational model for predicting drug responses since large-scale drug screening data have been accumulated. In this study, we have employed a deep neural network that has feedforward bypass connections to predict anticancer drug responses based on gene expression profiles and drug molecular fingerprints. Drug response data from Genomics of Drug Sensitivity in Cancer (GDSC) database were used. Log(IC50) was used as a prediction target, and this value was normalized for each drug such that all drugs have a same baseline. To divide the whole dataset into a training set and a test set, we undertook two different ways. One is randomly dividing cell lines for each drug into training/test sets, and the other is holding out a subset of cell lines. Our drug response prediction model reached an accuracy of 0.74 Pearson correlation coefficient under first training, test set division setting while 0.49 under second setting which is a higher performance compared to other drug response prediction methods. We also tried to predict missing values from the GDSC database against five MEK1/2 inhibitors and could successfully predict different sensitivities of those cell lines with BRAF mutation although the mutation information has not been used.

Acknowledgements: This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea Government, the Ministry of Science, ICT & Future Planning (2015M3A9A7067220, 2014R1A2A1A10052404, and 2013M3A9A7046303). It was also supported by the KAIST Future Systems Healthcare Project from the Ministry of Science, ICT & Future Planning

ABSTRACT. Machine learning techniques for cancer prediction and biomarker discovery can hasten cancer detection and significantly improve prognosis. Recent “OMICS” studies which include a variety of cancer and normal tissue samples along with machine learning approaches have the potential to accelerate such discovery. In this work, 2,175 gene expression samples from nine tissue types were used to identify gene sets whose expression is characteristic of each cancer class. Nine single-tissue, two multi-tissue cancer-versus-normal, and multi-tissue normal classifiers using random forests classification and ten-fold cross-validation were developed. The single-tissue models -- given a sample of a particular tissue type -- classified the sample as cancer or normal with a testing accuracy between 85.29% and 100%. A multi-tissue bi-class model, which classifies a sample as either cancer or normal achieved a testing accuracy of 97.89%, whereas, a multi-tissue multi-class model, which classifies a sample as cancer versus normal and as a specific tissue type achieved a testing accuracy of 97.43%. A multi-tissue normal model which, given a normal sample of any of the nine tissue, classifies the sample as a particular tissue type achieved a testing accuracy of 97.35%. This study demonstrates the feasibility of predicting the tissue origin of a carcinoma in the context of multiple cancer classes. The machine learning classifiers developed in this study identify potential cancer biomarkers with sensitivity and specificity that exceed those of existing biomarkers. Our tissue-specific biomarkers pointed to metabolic and signaling pathways that are critical to both general and tissue-specific tumor development.

ABSTRACT. In traditional herbal medicine (THM), different THMs may be prescribed for similar symptoms of a disease depending on the overall condition of the patient. To determine how to select the proper prescription for a patient, it is important to clarify the pharmacological properties of each THM prescription. In traditional Japanese medicine (Kampo), both Maobushisaishinto (MBST) and maoto (MT) prescriptions are commonly used to treat flu-like symptoms, but MBST is specifically used for patients with frailty. Both prescriptions share ephedra herbs (Ephedrae Herba (EH)) as the major ingredient, but their other constituents vary. The differences in the pharmacological properties of these Kampos remain unclear. In this study, we compared the pharmacological and compound traits of MBST and MT by target/non-targeted metabolomics analysis. Pharmacology: we evaluated the effect of MBST and MT on flu-like symptoms by using a rat model of polyI:C-induced inflammation. Both MBST and MT ameliorated the decrease in locomotor activity and proinflammatory cytokines in model rats, but the effects did not differ between the two Kampos. Targeted metabolomics: we compared changes in endogenous metabolites in rat plasma after MBST and MT treatment. Although both prescriptions affected several metabolic pathways such as amino acid metabolism, we identified metabolites that were affected specifically by only one of the two Kampos. Non-targeted metabolomics: we comprehensively compared the compound and metabolomic profiles of MBST and MT in rat plasma after treatment. Notably, one-quarter of EH-derived compounds were detected only in either Kampo, suggesting that the combination of herbs in each Kampo might affect the extraction and/or absorption of EH compounds. Our extensive pharmacological/compound profiling highlight differences between MBST and MT in chemical composition and in their effects on host metabolism. These analyses will lay the foundation for elucidating the mechanisms of the pharmacological/biological effects of THM.

ABSTRACT. The regulation of cell growth and division has long been considered a central problem in cell biology. Much work has been done both theoretically and experimentally towards expanding our understanding of the processes of the cell division cycle. However, cellular growth, which is regulated by a complex interplay of nutrient signaling and growth factor signaling pathways remains poorly understood in higher eukaryotes. Specifically, growth in unicellular eukaryotes, in contrast to multicellular eukaryotes, depends only on nutrient availability and not on growth factor signaling. However, the complexity of the interactions in the regulatory network for cellular growth precludes an intuitive understanding of nutrient signaling mediated regulation of cell growth.

Here, we propose an ODE-based dynamical model of the regulatory mechanism governing cell growth in the budding yeast Saccharomyces cerevisiae. This model can simulate variations in cellular growth rates as a function of varying macronutrient inputs, namely carbon and nitrogen sources. The model captures the metabolic signaling network composed of components from the TORC1, Snf1, and the Ras-cAMP-PKA pathways, and captures the interactions that impinge on the ribosome biogenesis regulon, which governs the mass growth rate. Previous dynamical models have focused on modeling the responses of pathways involved in specific nutrient responses, or have focused on general stress responses in yeast, involving an overlap with nutritional stress responses. Our model differs from existing efforts in its comprehensiveness in modeling the nutrient response, and its goal of predicting cellular growth rate robustly for a given nutritional input. In the future, we plan on integrating this model of nutrient signaling with the highly successful model of the yeast cell cycle, with the goal of refining the phenotypic predictions made by the cell cycle model.

ABSTRACT. Complex diseases can be modeled as damage to intra-cellular networks that results in abnormal cell behaviors. Network-based dynamic models such as Boolean models have been employed to model a variety of biological systems including those corresponding to disease. Previous work designed compensatory interactions to stabilize an attractor of a Boolean network after single node damage. We generalize this method to a multi-node damage scenario and to the simultaneous stabilization of multiple steady state attractors. We classify the emergent situations, with a special focus on combinatorial effects, and characterize each class through simulation. We explore how the structural and functional properties of the network affect its resilience and its possible repair scenarios. We demonstrate the method's applicability to two intra-cellular network models relevant to cancer. This work has implications in designing prevention strategies for complex disease.

ABSTRACT. In single cell time lapse studies cell segmentation and tracking can be time consuming. Using deep learning algorithms, we developed a pipeline that performs automatic segmentation and tracking after learning from a training set of data. We tested the method with data obtained with an E-cad/Vim double-FP-labeled T47D cell line using the CRISPR technique.

ABSTRACT. The development of drug resistance by acquired mutations is the main reason that melanoma cancers show only a modest response to single-agent targeted drug therapies. A well-designed combination of drugs that targets multiple molecular pathways can improve overall survival. To understand the molecular mechanisms underlying the development of drug resistance, we have created a comprehensive computational model of the signaling network in melanoma cells and how it responds to systematic perturbations with 90 different drug combinations. For this detailed mathematical model, we use a new modeling framework in which all reactions are classified into three basic types: protein synthesis and degradation (→ C →), phosphorylation and de-phosphorylation (C ↔ CP), and binding to activator or inhibitor partners (C+A ↔ C:A). The molecular mechanisms of signaling and cell cycle regulation in melanoma cells have been deduced from protein-protein interaction data obtained by an automated literature extraction method. We used published experimental data on 90 treatments consisting 12 different drugs, each at two different doses, and also different drug-pair combinations. Our initial approach was to train the model to fit experimental data (cell viability or inviability) on 12 drug treatments (either 12 low-dose treatments or 12 high-dose treatments) and compare the model predictions with data on the remaining 78 drug combinations. Trained in this way, the model is successful in explaining ~82% of viability data. We have also tested the model against detailed proteomic data, generating a response map for protein concentrations and an error map for comparison of experimental data with model predictions. We have parameterized the model using both viability and proteomic data, and we plan to use this version of the model to investigate drug resistance in melanoma and to predict multi-agent targeted therapies that might make it harder for melanoma cells to resist treatment. Our broader goal is to make a generic, customizable model of molecular signaling in cancer cells.

ABSTRACT. In recent years, it has become increasingly apparent that antisense transcription play an important role in the regulation of gene expression. Recently, it was reported that the antisense transcript of the mammalian core-clock gene Per2, which we named Per2AS, oscillates with a circadian period and about 12 h phase shift from the peak of expression of Per2 mRNA. In this study, we address the question as to whether Per2AS plays a regulatory role in the mammalian circadian-clock. In particular, we study the potential effects of Per2 and Per2AS interactions on the circadian rhythm in silico, in the context of two hypotheses about how Per2 and Per2AS transcripts mutually interfere with each other's expression. In our 'pre-transcriptional' model, we assume that the process of transcribing Per2AS RNA from the non-coding DNA strand represses the transcription of Per2 mRNA from the coding strand. In our 'post-transcriptional' model, we assume that Per2 and Per2AS transcripts form double stranded RNAs, due to their complementary sequences, and that the duplex RNA is rapidly destroyed. To study these alternative hypotheses, we have modified a mathematical model of the molecular regulatory network of the mammalian circadian clock, originally put forward by Relogio et al [1], by adding new terms describing our proposed hypotheses. Our pre-transcriptional model predicts that Per2-Per2AS interactions can generate alternative modes of circadian oscillations, and we characterize these modes in terms of the amplitude and phase of oscillation of various clock genes. In our post-transcriptional model, antisense overexpression dampens the circadian rhythm. In a model that combines pre- and post-transcriptional controls, the period, amplitude and phase of circadian proteins exhibit non-monotonic dependencies on the rate of expression of Per2AS, presumably as a consequence of the double regulatory functions of Per2AS.

[1]. Relógio A,et al., Tuning the Mammalian Circadian Clock: Robust Synergy of Two Loops. PLOS Computational Biology. 2011;7(12):e1002309.

ABSTRACT. Granulocyte-monocyte progenitor (GMP) cells play a vital role in the immune system as they mature into a variety of white blood cells, including neutrophils and macrophages. In the classical motif of GMP differentiation, GMP cells mature into one of two competing lineages, monocytes or granulocytes, depending on exposure to cytokines such as various types of colony stimulating factors (CSF). Granulocyte-CSF (G-CSF) induces granulopoiesis and macrophage-CSF (M-CSF) induces monopoiesis, while granulocyte macrophage-CSF (GM-CSF) favors monocytic and granulocytic differentiation at low and high concentrations, respectively. Although these differentiation pathways are well documented, the mechanisms behind the diverse behavioral responses of GMP cells to CSFs are not well defined. Using dynamic systems theory we explore the differentiation of GMP cells in response to varying dosages of G-CSF, M-CSF, and GM-CSF. Our model reproduces experimental observations of GM-CSF induced differentiation, and for the first time, we propose a mechanism for this intriguing behavior. Furthermore, we explore the differentiation of a fourth phenotype, monocytic myeloid-derived suppressor cells (M-MDSC), how they fit into the classical motif of GMP differentiation, and how progenitor cells can be primed for M-MDSC differentiation. Finally, we legitimize our model by comparing its results to numerous experiments and make intriguing predictions that should be explored by future experimental studies.

ABSTRACT. Small cell lung cancer (SCLC) is an aggressive tumor type with a strong ability to become resistant to all known treatments and to survive in diverse microenvironments. Proposals to stratify patients based on tumor phenotype have been met with resistance due to unclear clinical relevance, as the “small blue round” SCLC cells are extremely uniform by histopathology, but more recently it has become increasingly understood that SCLC tumors exhibit phenotypic heterogeneity implicated in the aggressiveness of the disease. My central hypothesis is that interactions between these plastic phenotypes form a functional ecosystem that drives growth and controls the overall response to therapy. The discovery of multiple SCLC phenotypes necessitates further analysis to identify the core SCLC phenotypes and eventually map the phenotypic space. Consensus clustering and weighted gene co-expression network analysis (WGCNA) applied to transcriptomics data of 53 cell lines and 81 primary human tumors reveal 4 clusters with several gene modules distinguishing them. The clusters found include the neuroendocrine tumor-propagating cell (NE TPC) and supporting mesenchymal-like (ML) phenotypes previously reported in the literature, as well as two novel hybrid phenotypes. We constructed an expanded panel of SCLC candidate markers for single-cell analysis (e.g., by mass cytometry) optimized to distinguish between and ensure broad coverage of these phenotypes. Characterization of the core SCLC phenotypes will contribute substantially to the goal of a global definition of the multi-dimensional variables that drive cooperation between them and support tumor progression and drug evasion.

ABSTRACT. Cells constantly process extracellular stimuli that can lead to multiple phenotypic outcomes. In this work, we study how L929 cells, stimulated with TNFa, execute death (programmed necroptosis) or survival through the Nuclear Factor kappa-B (NF-B) signaling network. Although the NF-B pro-survival network and the necroptosis network are typically studied in isolation, recent work by the Hoffmann lab has shown evidence that crosstalk between the networks exists. A mechanistic understanding of this crosstalk would thus provide insights about cellular intracellular communication that employs a complex system of interactions to achieve alternate outcomes. To understand the complex interactions that lead to TNF stimulation, we build and calibrate mathematical models of complex biochemical systems, as a tool to probe molecular interactions and mechanisms, outside the reach of current experimental technologies. In this work, we recapitulate the NF-κB model by, Hoffmann et al., a set of ordinary differential equations using mass-action kinetics, in PySB to ensure that our implementation recapitulates key aspects of NF-κB signaling. We then extend this model with a necroptosis execution module in PySB that extends the original NF-κB model to enable the exploration of programmed necroptosis, through the regulation of anti-apoptotic signals. We also intend to link the model with other models of necroptosis to better understand cell fate outcomes in cancers.

ABSTRACT. The cell cycle comprises a chain of events that results in the division of a cell into two daughter cells. It is carefully regulated by a complex network of control mechanisms including cyclin-Cdk interaction, DNA replication and checkpoints. In particular, the G2 checkpoint checks the integrity of the DNA before proceeding to mitosis. In the presence of DNA damage, the cell cycle is arrested in this checkpoint until the damage is repaired. Failure to activate or maintain this checkpoint causes genome instability and, in some cases, cancer cells. This is critically important in radiation therapy since it has been shown that G2 checkpoint activation is compromised for low doses of radiation [1, 2]. Here, we study the effect of radiation on the cell cycle through a mathematical model based on a Minimal Cdk Network [3]. In our modified version of the model, we include a DNA damage pathway and study its effect on the cell cycle (represented by a stable limit cycle). We identify the G2 checkpoint activation in the cell cycle with a saddle-node on an invariant circle (SNIC) bifurcation. For a small dose of radiation below a threshold, we observe that the period of the limit cycle increases (which corresponds to a delay in the progression into M-phase); for higher dose of radiation above the threshold, we observe a loss of the limit cycle and the appearance of a node and a saddle (corresponding to the activation of the checkpoint). We also observe that the G2 checkpoint, determined by the saddle point, is located right before mitosis and depends dynamically on the amount of radiation. Our results provide a foundation for understanding many phenomena observed in low-dose radiation, including hyper-radiosensitivity and increased radioresistance (HRS/IRR) phenomenon observed in the study of survival fraction after radiation.

[1] Brian Marples. Is low-dose hyper-radiosensitivity a measure of G2-phase cell radiosensitivity? Cancer and Metastasis Reviews, 23(3-4):197–207, 2004. [2] Markus Löbrich and Penny A Jeggo. The impact of a negligent G2/M checkpoint on genomic instability and cancer induction. Nature Reviews Cancer, 7(11):861–869, 2007. [3] Claude Gérard, John J Tyson, Damien Coudreuse, and Béla Novák. Cell Cycle Control by a Minimal Cdk Network. PLOS Computational Biology, 11(2):e1004056, 2015.

ABSTRACT. Eukaryotic single cells such as migrating cancer cells, developing neurons, or mating yeast move (chemotaxis) or grow (chemotropism) in a specific direction by using a chemical gradient for directional reference. Existing models for how eukaryotic cells sense chemical gradients assume cells compare concentrations of ligand-bound receptors across the cell surface to infer the gradient. However, our studies, using yeast chemotropic mating as a model, suggest that cells may sense a chemical gradient with a sensitive front that is mobile on the cell cortex. We have modeled the dynamics of this moving sensitive front with stochastically perturbed reaction-diffusion equations.

In this study, we show that an in silico cell based on this computational model is capable of sensing pheromone gradients. However, the model cannot sense shallow gradients as well as yeast in in vivo experiments. When we incorporate the effects of an additional pathway recently implicated in gradient sensing in yeast into our mathematical model, the model's ability to track pheromone gradients improves significantly. Furthermore, we experimentally demonstrate that the dynamics of a sensitive front that emerge from genetic perturbations recapitulate the dynamics generated by the mathematical model. In summary, we propose a mechanistic basis for a novel pathway necessary for tracking chemical gradients accurately.

ABSTRACT. The dermal extracellular matrix (ECM) plays a vital role in providing tensile strength, elasticity and hydration to the skin. The ECM is primarily composed of type 1 and 3 collagens but also consists of an array of various other components, the abundance of which optimizes skin function. With age, this balance is driven away from normal homeostasis resulting in less and disorganized collagen, reduced strength and elasticity and the other characteristics of skin ageing. In this study, time course expression profiles (12 time points over 96h, n=6) were collected from three cell lines (neonatal human dermal fibroblasts [HDFs], irradiation induced senescent HDFs [IR] and adult HDFs [56 years +], n=3 per cell line ) in response to transforming growth factor beta (TGFb) or control using both Affymetrix microarray and high throughput (WaferGen) PCR technology. Here we report on the analysis of a subset of this data (pilot data) which has been used to establish computational methods to extract meaningful information. We apply an integrated top-down and bottom-up systems biology methods using the ‘-omic’ level data to highlight differences between young, IR and adult cell lines with respect to the fibroblast TGFb response and ODE modelling strategies to build dynamic models of these differences. We aim to identify novel means of treating the aged phenotype.

ABSTRACT. Background: Administration of the cytokine interferon alpha (IFN-α) is a standard therapy against chronic hepatitis C infection- the leading cause of liver diseases. However, the success of the treatment is patient dependent. For optimised patient dependent IFN-α treatment design, mechanistic explanation of the underlying physiological processes across different scales of biological organisation is crucial.

Method: Using the case of IFN-α treatment in humans we here present a novel approach for the integration of molecular pathway models at the cellular level into physiology-based pharmacokinetic (PBPK) models at the organism scale.

Results: The multi scale model describes the whole-body distribution of IFN-α and the resulting cellular signalling response in the JAK/STAT pathway. It captures the non-linear pharmacokinetic behaviour of IFN-α within the body shedding light on the changes in signalling behaviour when considered an in-vivo context.

Conclusion: This work is a significant step towards understanding the mutual dependencies of the tissue specific pharmacokinetic availability of IFN-α and the resulting therapeutic response at the cellular signaling level. Moreover, it provides generic workflow for the integration of cellular models based on in-vitro data within an in-vivo context.

ABSTRACT. Visualization data tools help people to understand the dynamics of a mathematical model by establishing it in a visual context, instead of a set of mathematical equations. The purpose of the project is to help develop and test a computer visualization tool to highlight the characteristics of budding yeast cells during synthesis and mitosis. The mathematical model consists of a variety of differential equations and discrete simulations of eukaryotic cells. The budding yeast cell cycle visualization tool allows anyone to manipulate the cell by selecting different species, time, species scale, and tracking mother or daughter cells. In order for a budding yeast cell cycle to achieve its main goal, the cell must grow its size, duplicate its genes and proteins, and divide into two daughter cells, which is done by an alternative process of synthesis and mitosis. Not all budding yeast cells experience the same cycles. The data collected for the visualization tool were originally taken from numerical simulation, using Matlab, of a mathematical model of budding yeast cells. It was learned that in some mutant cases the cell would just continuously grow and die after the model was stimulated with the visualization tool. Some of the budding yeast cells may immerse in specific stages. Various cases would immerse mainly in the G1 phase. There were three different mathematical models for each case; labeled as Deterministic, Chen, and Hybrid. In the deterministic model, the corresponding yield will be produced every time. In the Chen model, we have the average behavior but not the variance. The Hybrid model is the stochastic and deterministic model combined. In conclusion of the current phase of the project, out of 122 mutant cases only 30 were visualized. More mutant cases are expected to be visualized in the next phase of the project.