The community of particular interest (COSI) in systems modeling (SysMod) organizes annual gatherings. This brief talk introduces all speakers, organizers, and the main topics of the 2022 meeting. This year’s meeting comprises three sessions covering various topics, beginning with “Methods for Modeling and their Applications,” followed by “Modeling Signal Transduction, Gene Regulation, and Protein-Protein Interactions.” It concludes with a session on “Applications of Single Cell Technology.” Two outstanding keynote speakers will present their visions on developments in these fields: Reinhard Laubenbacher from the University of Florida Health and Ana Conesa from the Politechnical University of Valencia, Spain. Bhanwar Lal Puniya will close the event by awarding this year’s poster prizes. The event is hosted by Andreas Dräger, Bhanwar Lal Puniya, and Juilee Thakar. This talk will end with an outline of the COSI team’s plans for a change in its governance structure, leading to more involvement of all members and restructuring of the COSI team.

Respiratory diseases typically involve processes across several spatial and temporal scales. The immune system plays an important role in a variety of ways, both in infectious and sterile conditions.  Computational models can help understand the mechanisms that link different scales together, and help discover new and compare existing treatment options. Their usefulness would be greatly enhanced if they could be personalized to characteristics of individual patients. This talk will describe such a model of the innate immune response to fungal and viral pathogens, and show some examples of how the model can be used as a virtual laboratory. The talk will also describe tools helpful for building large-scale models of this type.

MRA (Modular Response Analysis) is a method used to infer biological networks. From a set of independent perturbations applied on network nodes, it is possible to compute the connectivity between every pair of nodes. It is well-known that classical MRA is sensitive to measurement noise and perturbation intensity. One of the most important questions about MRA concerns discovered edges meaning. We have developed a new approach linking MRA and multiple linear regression. Connectivity coefficients estimation is equivalent to reckon regression parameters as confidence intervals. We have confirmed this approach successfully by comparing it with classical MRA, in the case of an “in silico” six nodes Map kinases network. One important MRA’s application, coupled to regression, is to identify null edges, notably concerning gene networks, which are often sparse. Many regression methods have been compared: multiple regression, “Lasso”, “threshold regression”, applied to “in-silico” gene networks stemming from Dream Challenge 4. Results have shown a correct error rate for 10 and 100 genes networks.

Metabolomics, synthetic biology, and microbiome research demand information about organism-scale metabolic networks. The convergence of genome sequencing and computational inference of metabolic networks has enabled great progress toward satisfying that demand by generating metabolic reconstructions from the genomes of thousands of sequenced organisms. Visualization of whole metabolic networks is critical for aiding researchers in understanding, analyzing, and exploiting those reconstructions. We have developed bioinformatics software tools that automatically generate a full metabolic-network diagram for an organism, and that enable searching and analyses of the network. The software generates metabolic-network diagrams for unicellular organisms, for multi-cellular organisms, and for pan-genomes and organism communities. The diagrams are zoomable to enable researchers to study local neighborhoods in detail and to see the big picture. The diagrams also serve as tools for comparison of metabolic networks and for interpreting high-throughput datasets, including transcriptomics and metabolomics data. These data can be overlaid on the metabolic charts to produce animated zoomable displays of gene expression and metabolite abundance. The BioCyc.org website contains whole-network diagrams for more than 20,000 sequenced organisms.

Human Milk Oligosaccharides (HMOs) are abundant carbohydrates fundamental to infant health and development and modulation of the infant microbiome. Although these oligosaccharides were discovered more than half a century ago, their biosynthesis in the mammary gland remains largely uncharacterized. Here, we use a constraint-based modeling framework that integrates glycan and RNA expression data to construct an HMO biosynthetic network and predict glycosyltransferases involved. To accomplish this, we construct models describing the most likely pathways for the synthesis of the oligosaccharides accounting for >95% of the HMO content in human milk. Through our models, we propose candidate genes for elongation, branching, fucosylation, and sialylation of HMOs. Our model aggregation approach recovers 2 of 2 previously known gene-enzyme relations and 2 of 3 empirically confirmed gene-enzyme relations. The top genes we propose for the remaining 5 linkage reactions are consistent with previously published literature. These results provide the molecular basis of HMO biosynthesis necessary to guide progress in HMO research, including the study of maternal genetics on HMO composition and metabolic engineering of microbes to facilitate the manufacturing of these molecules as invaluable nutraceuticals to improve infant health and development.

In recent years, the production of chemicals using biological systems has been gaining traction. Bioproduction is a sustainable alternative to traditional chemical processes, but it is often associated with higher operational costs. The economy of bioprocesses can be improved by co-producing multiple products using the same system. There are multiple algorithms that exist for in silico metabolic engineering of organisms to achieve overproduction of a single product. However, there is a lack of computational tools that can co-optimize a set of metabolites. Here, we present co-FSEOF (co-production using Flux Scanning based on Enforced Objective Flux), an algorithm designed to identify intervention strategies to co-optimize the production of multiple metabolites. Co-FSEOF can identify all pairs of metabolites that can be co-produced in an organism using a single intervention. It can also identify higher-order intervention strategies for a chosen set of metabolites. We have utilized this tool to identify intervention strategies for the co-production of pairs of metabolites in Escherichia coli and Saccharomyces cerevisiae under aerobic and anaerobic conditions. The proposed computational tool provides a systematic approach to study co-production and thereby aids the design of better bioprocesses.

Signaling crosstalk occurs when the stimulation of a signaling pathway’s receptors results in downstream effects on another signaling pathway. While numerous network-based methods identify the presence of signaling crosstalk, they do not distinguish between different signaling events and therefore offer limited mechanistic insights into crosstalk. Given the context-specific and often concurrent types of interactions in cellular signaling, multilayer networks offer the potential to better understand signaling pathways and predict their crosstalk. We built a multilayer network consisting of a gene regulatory layer and a signaling layer, and developed a statistical framework, MuXTalk, that uses high-dimensional edges, or multilinks, to model signaling crosstalk. Using statistically over-represented multilinks as proxies of crosstalk between signaling pathways, we identified potentially crosstalking pathway pairs among 61 KEGG pathways. In our benchmark, MuXTalk had a higher area under the ROC and precision-recall curves compared to all single layer-based methods tested, identifying additions to the current gold-standard. Crosstalk predictions in our “discovery” set of pathway pairs were highly supported in the literature with a precision >80% for the top 50 pairs. Overall, our findings suggest the utility of the multilayer modeling of signaling crosstalk, with possible future applications to extend our approach to tissue- and disease-specific crosstalk.

A holistic understanding of the biological interactions requires modelling integrated networks to elucidate normal and disease phenotypes. Despite signalling and metabolic networks being highly integrated, systems view of these networks allowed modellers to analyse them in isolation. In this work, a mathematical framework for integrating the signalling and metabolic networks is developed by integrating the constraint-based and ordinary differential equation (ODE)-based models. Particularly, metabolic networks are modelled using constraint-based modelling approach, while signalling networks are modelled using a set of ODEs. At each step, the signalling network’s information (activation or inhibition) is fed into the metabolic network as constraints on reaction rates. The changes in the metabolic network are then looped back into the signalling network as initial conditions to simulate the ODEs. This process is repeated till steady-state conditions arrive at the signalling network. Further, the integration approach is corroborated by integrating the insulin-glucagon signalling network with hepatocyte glucose metabolism. The developed framework captures the normal and disease phenotypes of the system. The algorithm developed is generic and can be employed to integrate signalling and metabolic networks for any biological system.

Connexins are the functional transmembranous components of cellular hemichannels and gap junctions, enabling both direct and long range metabolite exchange with other cells and the extracellular space. Connexin 43 is the most ubiquitous connexin in humans, supporting critical physiological functions. Its dysfunctional regulation is associated with a diverse set of pathophysiologies e.g. chronic wounds, psoriasis, cardiomyopathies, tumour progression and metastasis. Connexin mimetic peptides, e.g. Gap27, have been shown experimentally able to modify different connexin-based channel behaviours, offering promising therapeutic potential. In our work, we aim to explore the dynamical behaviour of connexin 43 cycling and consider how the introduction of Gap27 modulates this behaviour for potential experimental applications. We derive a new mathematical model for connexin 43 cycling with Gap27 binding based on mass action kinetics. We find through our numerical steady-state stability analysis, in-keeping with the biological literature, introduction of Gap27 affects only transient behaviour of connexin 43-based cycling species; the effect being Gap27 concentration dependent. This work provides a useful mathematical characterisation of connexin 43 cycling with Gap27 binding, reproducing qualitatively similar dynamical behaviour as reported experimentally. With further quantitative work, this model could be used in future collaboration with experimentalists to generate potential dosing regimen for Gap27.

Alphavirus is a genus of positive-sense single-stranded RNA viruses that cause numerous diseases, posing a significant threat to human welfare. For example, eastern equine encephalitis virus has a human mortality rate of 30% when neuroinvasive. Currently, there is no established cure for any alphavirus infection and only palliative care is available. Although numerous molecular components and interactions of alphaviruses within the host cell are characterized, how these mechanisms determine the dynamics of RNA viral replication and host immune responses remains unclear, limiting our ability to advance therapeutic development. We present a mathematical model to elucidate the mechanisms regulating the precise dynamics of alphavirus replication. Specifically, this model describes attachment, entry, uncoating, replication, assembly and export of both infectious virions and non-infectious virus-like particles within mammalian cells. The model recapitulates known characteristics of alphavirus infection, including the timing and amplitude of virion production, and identifies genome replication as the significant rate-limiting step. We are working to expand the model to incorporate the type I interferon induction signaling pathway. This will provide a comprehensive perspective on the conditions required for maximizing host response efficacy and determine key steps of immune system activation necessary for successful suppression of viral infection.

A vast amount of raw, heterogeneous and unstructured ‘omics’ data has boomed in the past few decades thanks to breakthrough advancements in genome sequencing. With such an influx of information, researchers strive to work together to understand details from single plasmids to entire eukaryotic organisms, necessitating a digital space for collaboration. Genome-Scale Metabolic Models (GEMMs) have emerged as a key player in encouraging laboratory scientists to store, analyze and validate experiments in computational formats. Such models stand out as the de-facto framework is mathematically defining metabolic activities within biological pathways. However, openly available tools for building, visualization and analysis of such networks are yet lacking. Cell Collective is an easy-to-use web-based interface to construct, visualise and simulate large-scale logical and metabolic models in a real-time and interactive fashion. Models can be easily loaded and exported in SBML or COBRA-like formats thereby making them compatible with already available tools. Such models and analysis workflows can easily be shared with other users for rapid knowledge-discovery, scientific reproducibility, conformity and interoperability within the modelling community at large. Cell Collective GEMM module maintains a publicly available model catalogue from widely used centralized repositories such as BioModels and BiGG with over 30,000 metabolic data.

During vertebrate somitogenesis an inherent segmentation clock coordinates the spatiotemporal signaling to generate segmented structures that pattern the body axis. Using our experimental and quantitative approach, we study the cell movements and the genetic oscillations of her1 expression level at single-cell resolution simultaneously and scale up to the entire pre-somitic mesoderm (PSM) tissue. From the experimentally determined phases of PSM cellular oscillators, we deduced an in vivo frequency profile gradient along the anterior-posterior PSM axis and inferred precise mathematical relations between spatial cell–level period and tissue-level somitogenesis period. We also confirmed a gradient in the relative velocities of cellular oscillators along the axis. The phase order parameter within an ensemble of oscillators revealed the degree of synchronization in the tailbud and the posterior PSM being only partial, whereas synchronization can be almost complete in the presumptive somite region but with temporal oscillations. Collectively, the degree of synchronization itself, possibly regulated by cell movement and the synchronized temporal phase of the transiently expressed clock protein Her1, can be an additional control mechanism for making precise somite boundaries.

Polycyclic aromatic hydrocarbons (PAHs) with several aromatic rings have gained growing importance since their presence is a common cause of soil pollution: they are slowly degrading, resistant substances, some of which are carcinogenic, meaning that if an area is contaminated with PAHs, neither agricultural production nor proper nature management is feasible there. This is why the process of disinfecting areas contaminated with PAHs has become an important task for microbiologists and bioinformatics researchers over the past few decades.
In our previous research, new potentially PAH-degrading enzymes were detected and applied to PAH-contaminated areas. It can be assumed that this application may cause changes in the ecosystem of these areas: after the difficult task of PAH degradation by the new enzymes is completed, other bacteria may appear in a more significant amount and accomplish the simpler phases of the degradation process. In this work, applying in silico methods (e.g., the AmphoraVizu webserver), the changes in ratio and numbers of short reads detectable in soil samples from different stages of enzyme treating are compared and interpreted in different phylotype levels.

Computational modeling and simulation become increasingly important to study self-organization and patterning during tissue morphogenesis. Morpheus (1,2) is an extensible open-source software framework (3) based entirely on declarative modeling. It uses the domain-specific language MorpheusML to define multicelllar models through a user-friendly GUI.
We here present how MorpheusML and the open-source framework allow
for rapid model prototyping and advanced scientific work-flows.
MorpheusML describes the various model components by a bio-mathematical language. It represents the spatial and mechanical aspects of interacting cells in terms of the cellular Potts model formalism and enables separation of model from implementation, model sharing, versioning and archiving in a public repository (4). Morpheus is SBML-compliant, it supports simulations based on experimental data and offers a broad set of analysis tools to extract features right during simulation. The associated FitMultiCell toolbox supports parameter estimation for stochastic Morpheus models (5). The rich C++ API allows to extend MorpheusML and the simulator with user-tailored plugins.
We demonstrate the use of Morpheus for liver systems medicine and regenerative processes in tissues and organs.

1. Starruß et al. (2014). Bioinformatics 30, 1331.
2. Morpheus homepage: morpheus.gitlab.io
3. Source code: gitlab.com/morpheus.lab/morpheus
4. Model repository: morpheus.gitlab.io/models
5. FitMultiCell homepage: fitmulticell.gitlab.io

While advances in next-generation repertoire sequencing technologies have revealed the dynamics of antibody diversification in both healthy and antigen-stimulated B cell donors, there is still a limited understanding of how antibodies achieve high specificity towards particular antigens. We investigate different methods to characterize antibodies specificity and their affinity to a particular epitope. Namely, we explore methods based on antibody sequence similarity, paratope identification and comparison, as well as antibody structural analysis as predicted by the sequence. We also exploit phylogenetic tree inference analysis to characterize the evolutionary processes that take place in Germinal centers (GCs), specialized compartments where B cells proliferate, differentiate, and mutate their antibody genes in response to foreign antigens. We leverage this framework to address central questions about GC evolution, such as whether spatially separated GCs generate antibodies that recognize similar antigens. Our study suggests that the re-engagement of B cell clones in different GCs contributes to the appearance of B cell clones in distinct GCs that react to the same antigen. In particular, our analyses indicates that communication between GCs happens regularly and contributes to the long-term maturation of antibodies, which has implication in the production of broadly protective antibodies targeting conserved, non-immunodominant epitopes.

The human nose and the respiratory tract harbour numerous bacteria and are a vital target of virus particles. COVID-19 demonstrated the threat of virus infections. The mutual interactions between bacteria, their host, and viruses influence the risk of an infection outbreak. Altering the microbiome could lead to novel therapeutic applications.

We systematically reconstruct genome-scale metabolic models of bacteria in the respiratory tract and refine existing modeling protocols and procedures. We contributed to standard file formats and software for benchmarking and efficient simulation, reconstructing cell-type-specific human cells, and combining models. So far, we have reconstructed eight species, including manual curation. By simulating the virus reproduction of COVID-19, we identified potential weaknesses.

Computational models provide a valuable basis for integrated research projects. Their predictions guide experimental work and identify potential targets for investigating complicated mutual interactions between microbes, viruses, and their host. By validating model-derived hypotheses, it may become possible to drastically reduce the growth and colonization ability to threaten pathogens from the nose to the lungs. We laid the foundation for systematic model building, quality control, open distribution, validation, and simulation. The outcomes of this project support the search for novel antimicrobial therapies against multi-resistant pathogens or severe viral infections.

To understand mechanisms that underlie the pathological processes of myocardial ischemia in humans, we investigated genetic-transcriptomic relationships in response to ischemia. Left ventricular biopsies from 118 patients undergoing aortic valve replacement surgery were obtained at baseline and after an average of 79 minutes of cold cardioplegic arrest/ischemia. Transcriptomes were quantified by RNA-sequencing.
To understand mechanisms of myocardial ischemia, here we investigated the effect of gene-ischemia interaction on gene expression response (reQTL) using the conventional approaches and a novel causal model.

To assess the role of these reQTLs in myocardial ischemic injury, we investigated the effect of genetic variants on troponin in a separate set of 1751 patient undergoing a coronary artery bypass graft surgery. Troponin is a complex of three regulatory proteins that is integral to muscle contraction in cardiac muscle and elevated levels in blood after cardiac surgery. We performed a genome-wide association study (GWAS) for myocardial ischemic injury during surgery, followed by two-sample multivariate Mendelian randomization and colocalization analysis between state-specific cis-eQTLs (baseline, post-ischemia, and response eQTL) and variants associated with heart diseases. We also applied other enrichment analysis including replication in GTEx, and finally mouse model.

Proteins perform the majority of essential biological processes governing cellular functions, yet the proteome remains largely unexplored at the resolution of single cells, representing crucial gaps in our understanding of cellular complexity. Here, we present a novel deconvolution algorithm that combines single-cell RNA-sequencing (scRNA-Seq) with bulk proteomics to model the global proteome at the single-cell level. Our approach leverages cell profile similarities to overcome the weak correlation between RNA-Seq and proteomics that confounds existing deconvolution strategies. We apply our algorithm to cell differentiation datasets and demonstrate its ability to accurately reconstruct single-cell profiles from bulk-level measurements at both the proteome and transcriptome levels. Furthermore, we show that our algorithm is able to successfully cross the protein-RNA divide by using scRNA-Seq in combination with bulk proteomics to distinguish established canonical markers. Our method provides a generalizable computational framework for charting the relationship between bulk and single-cell molecular layers, and offers researchers the ability to study the proteome at the single-cell level using established bulk-proteomics workflows. This work also lays the foundation for transferring cell-state information between RNA and protein modalities, integrating the under-served layer of proteomics into the single-cell analysis toolkit to enhance the prioritization of cell populations for targeted therapeutics.

Despite recent advances in single-cell omics analysis, studying metabolism at single-cell resolution remains challenging due to the lack of single-cell technologies for metabolites. In this study, we developed a method for predicting metabolic fluxes and metabolic state from single-cell transcriptome (scRNA-seq) data by combining machine learning (ML) modeling with flux balance analysis (FBA). Specifically, we built a ML model that incorporates the gene-protein-reaction (GPR) mapping from gene expression to reaction activity using a perceptron and a graph neural network (GNN) of metabolic reactions. For the GNN, we used the graph of genome-scale metabolic reaction network—each node is a reaction and the edges represent product-substrate relationship. We implemented a message-passing algorithm using reaction activity as the node attributes. These reaction activities were used to set the lower and upper flux bounds in the FBA. We demonstrated our method in an application to budding yeast, specifically to elucidate single-cell metabolic flux alterations caused by a shift in the sugar substrate in the culture media, from glucose to maltose. Our method identified a strong increase in the hexokinase pathway after the switch, consistent with previous reports of hexokinase suppression by glucose. The switch was also predicted to boost ethanol exchange and transport.

Genome-scale metabolic models (GEMs) can grant the key steps toward understanding the principles of microbes. The construction of the GEMs is technically possible despite being time-consuming. Recent studies have developed many automated network reconstruction tools to accelerate the reconstruction process1–6.
One of the more time-consuming steps toward curation the GEMs is mass and charge balancing of metabolites and the reactions involved. The presence of metabolites included several chemical formulas or undefined side groups (e.g., alkyl groups -R), and inconsistencies in mass balances can lead the difficulty in computing the in silico biomass and the synthesis of protons or energy (ATP) without chemical basis7, 8.
To facilitate the mass and charge curation, we created an automated Python module named “MassChargeCuration (MCC)” to reduce this curation manually. This Python module manipulates reconstructions encoded in Systems Biology Markup Language (SBML) and adds information by using multiple existing data resources and consistently consolidating their information. This also generates a visual comparison to track the changes.
Next, we applied this module to the metabolic reconstruction of Corynebacterium tuberculostearicum. This module improves the standardization policy of SBML models to provide a high-quality GEM for a better understanding of the microbe’s future role in health and disease.

Chemotaxis, cell migration in response to chemical gradients, is known to promote self-organization of microbiological populations. Despite the great attention paid by modelists, chemotaxis in heterogeneous environments is still less studied. In this study, we extend a multi-agent phenomenological model, designed for Disctyostelium discoideum colonies, to provide new insights on their self-organization processes in environments with obstacles of various kinds. Consequently, we study dynamical features emerging from complex social interactions among individual cells belonging to amoeba colonies.
An appropriate metric to quantitatively estimate the gathering process is first defined. Then, we find that: a) obstacles play the role of local topological perturbation, as they alter the flux of chemical signals; b) physical obstacles (blocking cellular motions and fluxes of extracellular signals) and purely chemical obstacles (only interfering with chemical fluxes) elicit similar dynamical behaviors; c) a minimal program for robustly gathering simulated cells does not involve mechanisms for obstacle sensing and avoidance; d) random cell movements concur in preventing multiple stable clusters and improve the gathering efficacy. Hence, we speculate that chemotactic cells can avoid obstacles without needing specialized mechanisms. Social interactions are sufficient to guarantee the aggregation of the whole colony past numerous obstacles.

COVID-19 is one of the deadliest respiratory diseases, and its emergence caught the pharmaceutical industry off guard. This study presents a novel workflow to predict robust druggable targets against emerging RNA viruses using metabolic networks and information of the viral structure and its genome sequence. For this purpose, we implemented pymCADRE and PREDICATE to create tissue-specific metabolic models, construct viral biomass functions and predict host-based antiviral targets out of one or more genome sequences. We observed that pymCADRE reduces the computational time of flux variability analysis for internal optimizations. We applied these tools to the bronchial epithelial cells infected with SARS-CoV-2 and identified enzymatic reactions with inhibitory effects. The most promising reported target was the Nucleoside Diphosphate
Kinase (NDPK1), for which the literature reports inhibitors. Finally, we computationally tested the robustness of our targets in all known variants of concern, verifying NDPK1 ’s inhibitory effect. Since our workflow focuses on metabolic fluxes within infected cells, it is applicable for rapid hypothesis-driven identification of potentially exploitable antivirals concerning various viruses and host cell types.

The surfaces of eukaryotic cells are decorated with glycans: branched sugar polymers covalently attached to proteins and lipids. These molecules encode cell identity and play key roles in cell-cell interactions. Eukaryotic glycans are manufactured in the Golgi apparatus, so perturbations to Golgi organisation are expected to impact the glycans made by a cell. In this context, the Golgi acts as an information channel whose input is the perturbation, and whose output is the glycan distribution. We consider an experimental paradigm in which a one-dimensional input perturbation (X) is reflected as a two-dimensional output (Y1,Y2): the levels of two specific surface glycans in single cells. We construct an optimal Bayesian decoder to infer X from the measured values (Y1,Y2). The decoder has two qualitatively distinct regimes: an all-or-nothing regime that decodes to the minimum and maximum values of X, and a continuous regime that decodes to intermediate values of X. We show that this behaviour, analogous to a second-order phase transition, arises due to the non-linearity of the input-output map X -> (Y1,Y2). We estimate a lower bound for the Shannon information capacity of the Golgi channel in this artificial setting.

Analysis of viral load profiles leads to insights into dosage regimens, duration of efficacy, and timing of treatment, aiding the drug development process. A SARs-CoV-2model by Goyal et al. (Fig. 1a) characterized the off-treatment viral dynamics and showed potent and early treatment is important for clearance of the virus. An asymptotic analysis of the short-term viral dynamics of that model in the prophylactic setting was conducted. Fig 1b illustrates a general case with nullclines (red: infected cells, blue: virus) and stability of the fixed points under prophylaxis. The viral entry inhibition under the effect of a hypothetical neutralizing antibody (Fig.9 of Goyal et al)was considered. Analytically, (i)a necessary condition for the clearance of the infection and(ii) minimal inhibition required to impede the entry of the virus into susceptible cells were derived. The viral dynamics will converge to the stable equilibrium (Fig 1c) under no viral inhibition and (Fig 1d) showed an entry inhibition of 80% and more is required to reduce the viral copies to nearly zero levels, thereby achieving successful prophylaxis.

Footprint-based methods have been available for decades in the context of omic data sets such as transcriptomic and phosphoproteomic to estimate transcription factor and kinase activity deregulations. Adapting such methodology for metabolomic data isn’t trivial, because the effect of a metabolic enzyme can be observed not only at the level of its direct reactant and products, but also in the abundance of metabolites located several reaction steps away from the metabolic enzyme of interest. ocEAn is a method that defines global metabolic enzyme footprint from reaction networks and explores coordinated deregulations of metabolite abundances with respect to their relative position in the reaction network in the same manner as TF-targets Enrichment analysis. ocEAn maps the relative position of every metabolite and enzyme and can be used to identify metabolic deregulations upstream/downstream of metabolic enzymes, even if the deregulated metabolites are located multiple reaction steps away. Each metabolic enzyme score represents the imbalance of metabolic abundance deregulations upstream/downstream of the metabolic enzyme. We applied ocEAn with cell line models of kidney cancer metastasis. We identified how nitrogen partitioning downstream of Branched Chain Amino-acid Transaminase BCAT1 through Argininosuccinate synthase ASS1 can be critical to support metastatic transition.

Influenza A viruses (IAV) drive one of the world’s most common viral infections. Although most cases resolve, about 650,000 individuals die due to respiratory complications each year. IAV primarily infects epithelial cells. To shed light on the molecular first-lines of immune defense against IAV infection, we integrated (1) gene regulatory network (GRN) modeling with (2) single-cell (sc)RNA-seq and scATAC-seq profiling of an IAV infection time course in a human airway epithelial (HAE) tissue model. The single-cell measurements enabled resolution of the heterogeneous cell types composing the tissue model, as well as IAV-infected versus uninfected “bystander” cells in the same culture. Both scRNA-seq and scATAC-seq informed the GRN. GRNs describe regulatory interactions between hundreds of transcription factors (TFs) and thousands of genes. To improve inference in this high-dimensional setting, we (1) assume GRN sparsity and (2) integrate TF binding site predicitons from the scATAC-seq, to guide GRN model construction from the scRNA-seq data. Then, we identified “core” subnetworks within the GRN, or regulatory interactions specific to cell types and conditions (e.g., 12h IAV response in basal cells). Our GRN analyses recovered canonical TFs mediating the interferon response such as STAT1, IRF1, IRF7, IRF9, as well as lineage-specific TFs such as MAFF.

The tumour micro-environment (TME) defines a complex system containing multiple cell types interacting through contact and cytokine exchanges. Particularly, immune cells play a major role in cancer development and their characterization allows a better understanding of the TME. In this context, transcriptomics time courses allow studying the gene regulatory networks and interactions between immune and cancer cells. Such computations rely on using these networks to obtain relevant information about the biology behind them, and to identify novel molecular interactions and potential drug targets.
In this project, we aim to characterise the formation of Nurse Like Cells (NLC) macrophages found in Chronic Lymphocytic Leukaemia (CLL) and to investigate the cross-talk between them. We performed GRN inference on both cell types, using advanced inference methods on a unique transcriptomics time-series on purified cell types from a 13-days co-culture. We applied network analysis and identified the main regulators through topological analysis of the inferred directed networks. Furthermore TF enrichment analysis shed light on the processes taking place inside the two cell populations. We are currently exploring how to integrate these entirely data-driven results with information from databases to turn these networks into executable logical models.

Biological systems have been under varying selection pressures through the course of evolution. In microorganisms, these selection pressures span from selection for viability in the face of environmental perturbations to antibiotic resistance. Robust biological systems contain multiple alternate routes for a given metabolic function. These alternate routes, upon inspection, reveal complex inter-dependencies within and across metabolic modules. Higher order synthetic lethals are excellent examples of these alternate routes. A synthetic double lethal comprises of a reaction pair which when deleted simultaneously, result in cell death; however, individually deleting one of the reactions does not give rise to the lethal phenotype. In order to ensure viability when individual reactions belonging to a lethal pair are deleted, there is a dynamic redistribution of fluxes. Analyzing this redistribution would help us to uncover the linkage between the reactions in a synthetic double lethal and also to understand the complexity underlying the reroutings.

We propose a novel method, minRerouting, which performs an $l_0$-norm MOMA-like formulation on the flux distributions obtained after deletion of individual reactions. Using this optimization approach, we simultaneously obtain the minimal rerouting set and the optimal flux distribution for both deletions, that are closest to each other in the flux space.

Microbes exhibit rich metabolic diversity leading to a broad set of behaviors and exposing a large substrate for bioengineering. To predict or design new biological functions, we utilize metabolic models. Building models relies mostly on empirical analysis, and there are few systematic automatic approaches to identify holes or potential errors. This could impact the quality of our predictions. Software engineering has developed methodologies to analyze, test, and predict behavior of large, complex systems which have similar properties to biological networks. In this work we utilize a common technique for inferring behavior in software to identify potential problems in metabolic models by developing a mapping between software and metabolic models. Invariants are logical relations between variables that remain constant across all program inputs. They can be used to verify correctness and infer program specifications. We extracted invariants from two models of E. coli with different fidelity: one of which is the iML1515 model. Some invariants point to potential model problems. For example, the total concentration of fluoride in extracellular space, periplasm, and cytoplasm remains constant in the iML1515 model, contradicting what we expect in nature. Initial results suggest invariants could be useful in refining metabolic models.

All populations display patterns in allele frequencies over time. Some alleles cease to exist, while some grow to become the norm. These frequencies can shift or stay constant based on the conditions the population lives in. If in Hardy-Weinberg equilibrium, the allele frequencies stay constant. Most populations, however, have bias from environmental factors, sexual preferences, other organisms, etc.

​We propose a stochastic Markov chain model to study allele progression across generations. In such a model, the allele frequencies in the next generation depend only on the frequencies in the current one.

We use this model to track a recessive allele through successive generations. Eventually, the allele will be “cancelled out” by the genotype of an organism becoming homozygous dominant. We estimate the number of generations it will take for this allele to be “cancelled out” by computing a hitting time in the Markov chain. This will allow us to efficiently communicate the trends of allele frequencies and estimate the speed of growth or decay of alleles.

The Toll-like receptor (TLR) signaling pathway is crucial for the initiation of effective innate immune responses. In this investigation, experimental and computational techniques are being integrated to generate a strongly data driven model of the TLR pathway. Targeted mass spectrometry was used to measure the absolute abundance of 54 (phospho)proteins (using 136 unmodified peptides and 29 phosphopeptides). The protein abundances ranged from 1,332 to 227,000,000 copies per cell (mouse bone marrow-derived macrophages). They moderately correlated with transcript abundance values (r = 0.699, p = 1.37e-17), and these data were used to make proteome-wide abundance estimates. Hundreds of TLR pathway protein-protein association rates were estimated using protein structures and molecular simulations (TransComp and Simulation of Diffusional Association). Rule-based pathway modeling and simulation is being performed using the Simmune software suite. The obtained values for absolute protein abundances and protein-protein interaction rates are being used as model parameters, and targeted phosphoproteomics is being used for model training, testing, and validation. This work was supported by the Intramural Research Program of NIAID, NIH.

The increasing number of available large RNA-seq datasets, combined with genome-wide association studies (GWAS), differential gene expression (DEG) studies, and gene regulatory networks (GRN) analyses have led to the discovery of many novel therapeutics. Despite this progress, our ability to translate GWAS and DEG analyses into an improved mechanistic understanding of many diseases remains limited, as both analyses disregard information about the cell types where the causative mechanisms driving the disease take place. This is critical because regulatory mechanisms may differ widely across cellular types. We explore several independent approaches to elucidate cell-type specific regulatory information about candidate genes associated with rheumatoid arthritis (RA). We compute sample-specific GRNs, which is a substantial advance compared to cohort-specific GRNs, as it enables the use of statistical techniques to compare network properties between phenotypic groups. Our analysis makes very precise experimental predictions, such as the impact of the knockdown of a specific TF in a specific cell type, and therefore, we expect it to be very useful to both rheumatologists and the broader scientific community interested in identifying cell-specific driver genes in other complex diseases.

Many proteins only function when in complex with other proteins, yet experimental methods for protein-protein interaction (PPI) determination lack the ability to accurately construct interaction networks on a proteome-wide scale. To address this shortcoming, computational methods can be used to complement high throughput experimental datasets as well as motivate small-scale biochemical experiments. We have developed such an approach, PEPPI, which predicts PPIs on a whole-proteome scale through a combination of homology, functional association, and machine learning; we find that PEPPI shows superior performance when compared with current state-of-the-art methods while remaining species-agnostic and computationally efficient. With PEPPI, we have predicted several cross-species interaction networks, including the network between humans and SARS-CoV-2 virus as well as the network between humans and the probiotic bacterial strain E. coli Nissle. Ongoing refinements include incorporation of a deep learning-based approach which leverages interchain co-evolution to make its predictions. Finally, to more faithfully represent the complexity of PPI networks, we have developed a 3D network visualization program, FALCON, which supports visualization of networks in immersive virtual reality space. Through the understanding these programs provide of protein function expressed through PPI networks, we can more effectively develop therapeutics through future drug development and protein engineering studies.

Cells and tissues respond to perturbations in multiple ways that can be sensitively reflected in alterations of gene expression. Current approaches to finding and quantifying the effects of perturbations on cell-level responses over time disregard the temporal consistency of identifiable gene programs. To leverage the occurrence of these patterns for perturbation analyses, we developed CellDrift (https://github.com/KANG-BIOINFO/CellDrift), a generalized linear model-based functional data analysis method capable of identifying covarying temporal patterns of various cell types in response to perturbations. As compared to several other approaches, CellDrift demonstrated superior performance in the identification of temporally varied perturbation patterns and the ability to impute missing time points. We applied CellDrift to multiple longitudinal datasets, including COVID-19 disease progression and gastrointestinal tract development, and demonstrated its ability to identify specific gene programs associated with sequential biological processes, trajectories, and outcomes.\

Single cell RNA-seq assays have dramatically advanced our ability to study and model cellular differentiation and cell fate decision. RNA velocity is an analysis framework that has become fundamental in the toolbox of the single-cell research community. However preprocessing choices and model assumptions of current RNA velocity models can dramatically influence their predictions and lead to misinterpretation of developmental order and cell fate. To this end, we propose Pyro-Velocity a probabilistic and end-to-end inference framework for RNA velocity based on variational inference that is scalable to millions of cells, based on unsmoothed data, and that naturally provides uncertainty estimation of cell fate based on a joint learning of transcriptional processes across genes. In addition, Pyro-Velocity can be used to learn cell fate decisions from one dataset and predict RNA velocity to other datasets with similar biological contexts. Our method outperforms existing methods in predicting cell fate on a compendium of single cell RNA-seq datasets with or without lineage information and for different biological systems. In summary, Pyro-Velocity is a new probabilistic RNA velocity framework and a user-friendly end-to-end software package to study cell fate decision from single cell data.

Despite the abundance of multi-modal data, suitable statistical models that can improve our understanding of diseases with genetic underpinnings are challenging to develop. Here we present SparseGMM, a novel statistical approach for gene regulatory network discovery. SparseGMM uniquely uses latent variable modeling with sparsity constraints regulators to learn gaussian mixtures from multi-omic data. By combining co-expression patterns with a Bayesian framework, sparseGMM quantitatively measures confidence in regulators and uncertainty in target gene assignment by computing gene entropy. We apply SparseGMM to liver cancer and normal liver tissue data and evaluate the discovered gene modules in an independent scRNA-seq dataset. sparseGMM identifies PROCR as a regulator of angiogenesis, and PDCD1LG2 and HNF4A as regulators of immune response and blood coagulation in cancer, respectively. Additionally, we show that more genes have significantly higher entropy in cancer compared to normal liver; among high entropy genes are key multifunctional components shared by critical pathways, such as p53 and estrogen signaling.

Aging profoundly affects immune system function, rendering the elderly more susceptible to pathogens, cancers and chronic inflammation. Single-cell genomics studies have accelerated the discovery of age-dependent immune-cell populations, linking aging phenotypes to changes in diverse immune populations. Here, we used single-cell RNA-seq (scRNA-seq) and chromatin accessibility (scATAC-seq) to deeply profile the CD4+ memory T cell (CD4+TM) compartment over time, enriching for IL10+ cells. We captured many T cell subsets including a population of IL10-producing, T follicular helper-like cells (Tfh10) we previously linked to suppressed vaccine responses in aged mice. From these data, we inferred gene regulatory networks (GRNs) and predicted transcription factor control of gene expression across T cell subsets in youth and old age. Further, we integrated pan-cell sc-genomics studies to identify factors from the microenvironment driving age-dependent changes in CD4+TM. Through computational modeling and broad integration across sc-genomics aging studies, our atlas of finely resolved CD4+TM subsets, GRNs and extracellular-signaling networks opens new opportunities to manipulate IL-10 production and improve immune responses in the elderly.

Multiomics experiments are now wide-spread and a powerful tool to assist in precision genomics. However, the transformation of multiomics data into a multi-layered representation of the biology is still a highly challenging task. I will review the challenges that field faces for multi-scale modeling and present novel statistical methods for interpretative multi-omics integration and for connecting epigenetic modifications to the control of cellular metabolism.

This talk briefly reviews this year’s SysMod community meeting, including speakers, chairpersons, and organizers. The day comprised three sessions on various methods and applications of computational systems modeling in biology. In addition, two keynote talks were given by Reinhard Laubenbacher and Ana Conesa. It is now the time to thank all contributions to scientific posters and bestow the best ones with the annual SysMod poster awards in 2022. Finally, we will conclude the meeting and open the discussion with the audience for feedback on the next SysMod in 2023. In particular, this involves a discussion on the transition of this COSI’s governance structure for improved involvement of the community members.