The community of particular interest (COSI) in systems modeling (SysMod) organizes annual gatherings. In 2021 the meeting comprises six sessions covering various topics, beginning with two sessions on disease and multi-scale modeling, followed by a dedicated session on infectious diseases. It concludes with two sessions on integrative approaches and methodologies and one session on structure-based dynamic modeling. As a highlight, three keynote speakers will outline the trend-setting developments in these fields: Ines Thiele, Ruth E. Baker, and Boris N. Kholodenko. Juilee Thakar will close the event by bestowing this year’s poster awards. The event is hosted by Claudine Chaouiya, Anna Niarakis, Andreas Dräger, Laurence Calzone, and Matteo Barbaris on behalf of the team of all twelve COSI organizers. The first meeting day includes a virtual social event. This brief talk introduces all speakers, organizers, and the main topics of the 2021 meeting.

Precision medicine relies on the availability of realistic, mechanistic models that capture the complexity of the human body. Comprehensive computational models of human metabolism have been assembled by the systems biology community, which summarise known metabolic processes occurring in at least one human cell or organ. However, these models have not yet been expanded to connect with whole-body level processes. To address this shortcoming, we have built whole-body metabolic models of a male (deemed Harvey) and a female (deemed Harvetta) starting from the existing human metabolic models, physiological and anatomic information, comprehensive proteomic and metabolomic data, as well as biochemical data obtained from an extensive manual literature review. We tested the predictive capabilities of the resulting whole-body metabolic models against the current knowledge of organ-specific and inter-organ metabolism. The final models contain 28 organs. Importantly, these whole-body models can be expanded to include the strain-resolved metabolic models of gut microbes. By parameterising the whole-body metabolic models with physiological and metabolomic data, we connected physiology with molecular-level processes through networks of genes, proteins, and biochemical reactions. As a sample application of the whole-body metabolic models, I will demonstrate how different microbial composition leads to differences in host metabolism, such as the capability to produce important neurotransmitters in the brain and flux through liver enzymes, with implications for the gut-brain axis as well as for microbiome-mediated liver toxicity. The predictions were consistent with our current understanding but also highlighted that different microbiota composition can lead to high inter-person variability. I envisage the microbiome-associated whole-body metabolic models will usher in a new era for research into causal host-microbiome relationships and greatly accelerate the development of targeted dietary and microbial intervention strategies.

The human nose harbors diverse microbes that play a fundamental role in the well-being of their host. Still, the contribution of many nasal microorganisms to human health remains undiscovered. Among all microbial species in the human nose, Staphylococcus aureus belongs to the most common human nasal pathogens.  Multiple epidemiological studies identify Dolosigranulum pigrum as a candidate beneficial bacterium based on its positive association with health, including negative associations with S. aureus.In this work, we propose a workflow for understanding the composition of the nasal microbiome community and the intricate interplay from a metabolic modeling perspective. To this end, we first create a basic community model that mimics the human nasal environment. Then we incorporate accurate genome-scale metabolic network models of D. pigrum and S. aureus.Our analysis supports the role of negative microbe-microbe interactions involving D. pigrum examined experimentally in the lab. By this, we identify and characterize metabolic exchange factors involved in a specific interaction between D. pigrum and S. aureus as an in silico candidate factor for a deep insight into the associated species. This method can open developing new ways to inhibit S. aureus and corresponding disease-causing infections through microbe-microbe interactions.

A central goal in systems biology is to understand how all of the molecules and processes in cells interact to generate behavior. While small molecules can be concretely described by molecular formats such as SMILES, macromolecules are frequently described as sequences of canonical residues. However, non-canonical residues, caps, crosslinks, and nicks are important to many functions of DNAs, RNAs, proteins, and complexes. One barrier towards models that explain how networks of such non-canonical macromolecules perform complex functions is our limited formats for concretely describing them. To overcome this barrier, we develop BpForms and BcForms, a toolkit for unambiguously representing the primary structure of macromolecules as combinations of residues, caps, crosslinks, and nicks. The toolkit can help omics researchers perform quality control and exchange information about macromolecules, help systems biologists assemble global models of cells that encompass processes such as post-translational modification, and help bioengineers design cells.

Motivation: Mathematical models of biological processes are often represented as complex networks of signaling pathways, describing intracellular behaviors of specific cell types (epithelial, T cells, etc.). However, this representation prevents us from describing spatial information or cell-cell interactions, that plays an important role in the dynamics and gives a more complete view of a biological process. To test the effectiveness that comes from merging models of signalling pathways with spatial models of cell populations, we present a model of cell invasion made with PhysiBoSS, a multiscale framework which combines agent-based simulation and continuous time Markov processes applied on Boolean network.Methods: The purpose of this model is to study the different modes of cell migration through an extracellular matrix by combining the spatial information obtained from the agent simulation and the intracellular information obtained from the possible stable states of the transcription factors signalling network. This includes different pathways involved in cell fates processes that can lead to death, proliferation, quiescence and invasion.Results: The model allows to simulate different initial conditions and mutations, as well as monitoring each cell behaviors using 2D and 3D representations, successfully reproducing single, collective and trail migration processes.

Psoriasis arises from complex interactions between keratinocytes and infiltrating immune cells in the skin. Chronic psoriatic inflammation is perpetuated by a Th17-dependent intercellular signaling loop, and pro-inflammatory eicosanoids are hypothesized to play a role in this process. We aimed to emulate the regulatory network of psoriatic keratinocytes with a logical model representing current knowledge about disease mechanisms and use this to study modes of action of cPLA2 inhibitors, alone or in combination with other drugs, and explore the cytokine-mediated signaling that takes place in psoriatic lesions. Through the integration of in vitro and in silico experimentation, we describe the PLA2-dependent release of PGE2 in response to Th17 cytokines. Further analyses of the computational model revealed the immunomodulatory role of Th1 cytokines, the modulation of the physiological states of keratinocytes by Th17 cytokines, and how together they promote the development of psoriasis. The response to treatment with cPLA2 inhibitor and/or Calcipotriol revealed a distinct mode of action of the two drugs. Lastly, novel entities were identified as potential drug targets that could restore a normal phenotype. In addition to contributing to the knowledge about psoriasis, this work showcases how the study of complex diseases can benefit from integrated systems approaches.

Type 2 Diabetes (T2D) is a major epidemic characterized by an increased blood glucose resulting from a defect of insulin secretion and insulin sensitivity. Available drugs targeting these mechanisms do not cure T2D. The intestinal absorption of ingested glucose is an underestimated but significant contributor to T2D. Since the 1960, hundreds of mathematical models have been proposed with the overall objective of predicting the dynamics of glucose homeostasis with increasing accuracy. Despite decades of research and massive funding, the underlying mechanisms of T2D remain poorly understood.In this preliminary work, we consider one of the most cited model of postprandial glucose response and explore its limits. First, we evaluate this model to predict an original dataset of obese diabetic patients. Second, based on partial parameter estimation, we investigate the capability of this model to predict which physiological compartments (intestine, pancreas, liver, etc.) are most likely able to restore a normal glycemia from a pathological one. Considering the parameter values given by the original authors of the model, we show that the model predictions are as expected. However, considering our own experimental dataset, the model fails to predict how the glucose homeostasis is altered after a physiological modification of intestinal absorption.

CD4+ T cells provide cell-mediated protection against diseases. When dysregulated, CD4+ T cells are associated with autoimmune and other immune-mediated diseases. Metabolism of CD4+ T cells regulates their function, therefore offer an opportunity to explore as a drug target against autoimmune diseases. In this study, we developed constraint-based models of naive and T helper 1, 2, and 17 subtypes. We mapped existing drugs and compounds and simulated metabolic behaviors under drug-induced inhibitions of metabolic genes. We integrated these metabolic behaviors with gene expression data of three autoimmune diseases, rheumatoid arthritis (RA), multiple sclerosis (MS), and primary biliary cholangitis (PBC). We identified and prioritized drugs and their targets that reversed the directions of differentially expressed genes in diseases. We identified 68 metabolic drug targets for the three studied diseases. We performed in vitro experiments and mined experiments available in the literature to validate results. The experimental results showed that 50% of the drug targets suppressed CD4+ T cell proliferation. In the end, we developed an integrated pipeline to explore metabolic models to identify drug targets and repurposable drugs.

We are closing the session on “Disease and multi-scale modeling” with a brief intervention from experts of the field to launch a discussion on current limitations and future challenges for the effective modeling of diseases and multi-scale biological processes. Finally, the panel will conclude with a brief review of the first two SysMod sessions and prepare the audience for the upcoming infectious diseases sessions.

Respiratory viruses present major health challenges, as evidenced by the 2009 influenza pandemic and the ongoing SARS-CoV-2 pandemic. Severe virus respiratory infections often correlate with high viral load and excessive inflammation. Understanding the dynamics of the innate immune response and its manifestation at the cell and tissue levels are vital to understanding the mechanisms of immunopathology and developing strain independent treatments. Here, we present a multicellular spatial model of two principal components: RNA virus replication and type-I interferon mediated antiviral response to infection within lung epithelial cells. The model exhibits either linear radial growth of viral plaques or arrested plaque growth depending on the local concentration of type I interferons. Modulating the phosphorylation of STAT or altering the ratio of the diffusion constants of interferon and virus in the cell culture could lead to plaque growth arrest. The dependence of arrest on diffusion constants highlights the importance of developing validated spatial models of cytokine signaling and the need for in vitro experiments to measure these diffusion constants. Findings suggest that plaque growth and cytokine assay measurements should be collected during arrested plaque growth, as the model parameters are significantly more sensitive and more likely to be identifiable.

Typical analysis of large-scale data utilizes identification of differentially expressed genes and enrichment analysis for mechanistic information. Enrichment analysis uses predefined sets of genes such as gene-sets from literature. Gene-sets do not have network topology or regulatory information. Enrichment analysis is limited in single-cell transcriptomic data, fewer genes are measured compared to bulk sequencing. To enhance regulatory information from gene-sets we are developing an algorithm optimizing gene-set subnetworks, using differential network analysis across experimental groups. Specifically, we use a genetic algorithm(GA) to optimize subnetworks of co-expression networks informed by gene-sets. The optimization function utilizes prior information, such as transcriptional regulation. We investigate novel regulators of atherosclerosis(AS) in people living with HIV. In a set of cells expressing naive T-cell markers we identified associations among 360(AS-) and 304(AS+) genes using mutual information(MI)(MI>(MIµ±2σ)). These contain genes related to an inflammatory response in HIV infection. Gene-sets from MSigDB were selected based on overlap with co-expression networks. After optimization, gene-sets with differential degree centrality were identified which reveal cell type-specific dysregulation of gene-sets and can be used to determine novel targets for interventions and future experiments. STREGA-NONA integrates knowledge with data to identify gene-set specific topology and differential activity using single-cell RNAseq data.

Mucins are present in mucosal membranes throughout the body and play a key role in microbe clearance and infection prevention. Thus, understanding the metabolic responses of pathogens to mucins will further enable the development of protective approaches against infections. We update the genome-scale metabolic network reconstruction (GENRE) of one such pathogen, Pseudomonas aeruginosa PA14, through metabolic coverage expansion, format update, extensive annotation addition, and literature-based curation to produce iPau21. We then validate iPau21 through Memote, growth rate, carbon source utilization, and gene essentiality testing to demonstrate its improved quality and predictive capabilities. We then integrate the GENRE with transcriptomic data in order to generate context-specific models of P. aeruginosa metabolism. The contextualized models recapitulated known phenotypes of unaltered growth and differential utilization of fumarate metabolism while also providing a novel insight into increased propionate metabolism upon MUC5B exposure. This work serves to validate iPau21 and demonstrate its utility for providing novel biological insights.

A highly curated genome-scale metabolic model of cystic fibrosis (CF) pathogen, Pseudomonas aeruginosa PA14 was developed in this study. During the reconstruction process, a set of substrate utilization and gene-essentiality data was used to improve the predictive abilities of the model. Furthermore, strain-specific processes (e.g., phenazine transport and redox metabolism, cofactor metabolism, carnitine metabolism, oxalate metabolism, etc.) were added to the reconstruction after a thorough literature review. Through this extensive process, a three-compartment, BiGG model, iSD1511 was created which is a significant improvement over the previous modeling effort for this strain. iSD1511 was assessed using another set of gene essentiality and substrate utilization data, and the prediction accuracies as high as 92.7% and 93.5%, respectively, were computed. The model can simulate growth in both aerobic and anaerobic conditions. iSD1511 can also be used to simulate condition-dependent production of phenazines and their effect on the growth of PA14 strain. Finally, the model was used to provide mechanistic explanations for two case studies: a) effect of gene deletion on gluconate production, and b) metabolic influence on drug tolerance of P. aeruginosa. Overall, this work provides a computational framework to aid in the development of effective intervention strategies against P. aeruginosa.

Antigen-presenting cells (APC) are immune cells mediating immune response by processing and presenting antigens to lymphocytes such as T cells. Dendritic cells (DC) are considered critical APCs that play an essential role in immune responses. Developments in Omics-technologies revolutionized the research in biological fields, which significantly increased the demand for integrative approaches to address novel predictions and analysis. One of these approaches is computational modeling that provides an executable, dynamic network revealing system-level behaviors through in silico perturbations. The in silico simulations provide a platform for researchers to generate new or rank existing hypotheses before starting the wet-lab experiments, saving time, energy, and materials.Our DC model is a manual curated, entirely annotated, and expert validated dynamic network comprising different cellular compartments to facilitate tracking flow of information in the system (e.g. extracellular stimuli, surface markers, signaling pathways, transcription factors, and secreted cytokines). The model is built in the Cell Collective platform utilizing logical formalism and contains more than two hundred nodes and interactions illustrating the APC function in DC.
Taking together, our DC cell-specific model of APC functions will address the links between the innate and adaptive immunity crosstalk.

Typically, clinical metabolome analysis is performed on blood samples. However, drawing blood is not only a cumbersome procedure for patients but requires qualified personnel which impairs measurement during real life settings. A promising alternative is the analysis of the metabolome from finger sweat, which benefits from simple sampling procedures. However, a major obstacle is the inability to control the amount of sweat produced by the sweat glands on the fingertips at any given time. Not addressing this problem prevents a reliable quantification of metabolites.Here we present a computational method based on the identification of metabolically linked species in the sweat metabolome that allows us to estimate sweat volumes and enables an individualized, accurate, quantitative finger sweat analysis for clinical applications. In a proof-of-principal application we used short interval sampling of sweat from fingertips to monitor the dynamic response of 43 individuals after caffeine consumption. We not only identified corresponding xenobiotics concentration time-series but extracted individualized kinetic parameters of caffeine metabolites from sweat and show their long-time stability.In conclusion, this work highlights the feasibility of individualized and reliable biomonitoring using sweat samples from fingertips which may have far reaching implications for personalised medical diagnostics and biomarker discovery.

Interconnected feedback loops in gene regulatory networks can give rise to sophisticated dynamics such as stepwise cell lineage commitment or low-rate irreversible differentiation. Despite their importance in biology, there is a lack of tools for extracting and visualizing these interconnected, “high-feedback” subnetworks from complex networks. Though many algorithms can enumerate or detect enrichment of specific motifs in networks, they cannot search for motifs defined by cycle interconnection. We developed a software toolbox HiLoop that can enumerate and intuitively display instances of high-feedback motifs even from large networks. With HiLoop, dynamics of selected subnetworks can be modeled with ordinary differential equations under varied parameter sets to characterize the landscape of possible attractors. HiLoop can compare a network to numerous permutations thereof to quantify the enrichment of high-feedback motifs. We used HiLoop to show that a network of genes involved in epithelial-mesenchymal transition is strongly enriched in a particular interconnection of three positive feedback loops. Modeling output of several selected subnetworks compactly yet clearly communicated the types of multistable and/or oscillatory systems that they are expected to produce. We contribute HiLoop as both a hypothesis generator and quantification method for exploring dynamics of complex feedback-rich regulatory networks.

Gene regulatory networks form an attractive avenue for understanding biological systems. Among gene regulatory networks, the lysis/lysogeny switch in bacteriophage lambda has received particular attention as the classical example of a genetic switch. We have previously found success in analyzing bacteriophage lamba using a hybrid system model formalism which reveals that the lytic and lysogenic states form stable attractors in our model’s state space. In our current work we further refine our hybrid system based modeling approach by expanding our phage model with additional data, analyzing the model’s overall robustness across a variety of conditions and network layouts. We apply and adapt our insights into the regulatory patterns of phage lambda to additional temperate phages, notably the HK022 lambdoid phage and Mu.

The Chemical Master Equation (CME) is a standard approach to model biochemical reaction networks. Alternatively, sample paths from the underlying stochastic process can be obtained by the Stochastic Simulation Algorithm (SSA). Both approaches are only applicable for small systems due to a large number of system states or high propensities.
In our study, we focus particularly on reactions that are driven by catalytic molecules, including for example epigenetic changes such as DNA methylation. We show that these complex systems can analogously be described by a concatenation of simpler systems. This simplification is achieved by simulating one catalyst molecule at a time instead of all molecules together. The advantages of our approach are:
(i) Intractable state transition graphs are replaced by a concatenation of much simpler graphs, resulting in lower dimensional linear differential equations for the CME approach.
(ii) The implementation of the SSA is considerably simplified.
We show the validity of our approach by applying it to two test-bed reaction systems, degradation of a molecule modelled with a simple first order kinetics, and methyltransferase-mediated DNA methylation. We furthermore discuss conditions under which our approach is applicable to a larger bandwidth of systems.

Defining differential equations at certain biological levels may be hampered by the lack of mechanistic details. Phenomenological equations are a rule of thumb to fill these knowledge gaps and get a model, but they come at the cost of being arbitrary choices. In this context, deep learning and hierarchical modeling can complement each other to automatically decipher trajectory patterns and variable interplay from data and underlying physics. Such integrated technologies may bridge low-level molecular descriptions and clinical manifestations as needed in quantitative systems pharmacology (QSP) applications.We developed a PyTorch pipeline that integrates a Long Short-Term Memory (LSTM) neural network and known ODEs to reproduce the dynamics of model variables with missing equations and to forecast future model states. Our approach is designed to generalize outside the training domain and does not require prior intuition on suitable regression classes nor dynamical laws. In our benchmarks, the variable interplay within the AI-ODE system and the adopted learning protocol improved the training procedure, leading to accurate predictions as a result. Further tests on relevant QSP applications and a modeling-oriented extension of the pipeline to account for model rate calibration are currently under investigation.

The reconstruction of genome-scale metabolic models (GEMs) is a growing field within systems biology. GEMs allow predicting varying phenotypic states. Thus, they have proven useful to guide metabolic engineering strategies of relevant biotechnological organisms. One microbe of particular biotechnological relevance is the Corynebacterium glutamicum, ideal for producing ʟ-amino acids at an industrial scale. We here present an updated high-quality GEM of the Corynebacterium glutamicum ATCC 13032.
The high reconstruction quality of this GEM is obtained by applying a high level of annotations and cross-references of metabolites, reactions, and genes. The GEM comprises 1,045 metabolites, 1,545 reactions, and 807 genes. Our model reproduces experimentally validated data and had realistic growth rates when tested on different media under aerobic and anaerobic conditions. Moreover, the model produces all canonical amino acids.
This new in-silico model is useful in filling knowledge gaps in C. glutamicum: ʟ-glutamate could still be produced even when the formerly believed relevant enzyme pyruvate carboxylase was knocked out. Our study shows that by integrating high reconstruction standards, GEMs can prove fruitful in consistently reproducing experimentally validated knowledge, filling knowledge gaps, and supporting metabolic engineers.

To enable earlier detection of compound toxicity in drug development, it is necessary to better understand the biological mechanisms. We aimed to identify early changes in pathway and transcription factor (TF) activity relevant for Drug-Induced Liver Injury which take place before adverse histological changes are observed. To do so, we used data from the TG-GATEs database, which comprises time-resolved transcriptomics and histopathology from repeat-dose studies in rats across eight timepoints, ranging from 3 hours to 4 weeks. We were able to recover known events preceding DILI, with some being more frequent, e.g. mitophagy, and others more confident, e.g. bile acid recycling. Furthermore, we prioritized additional pathways and transcription factors (TFs) based on time concordance. Among TFs, we were able to separate induced TFs, such as Cebpa, from post-transcriptionally activated ones, e.g. Srebf2, based on whether differential TF expression is observed before changes in regulon activity. Furthermore, we identified interactions between TFs which are supported by functional interactions as well as time concordance providing hypotheses on toxicity-related gene-regulatory mechanisms, e.g. Hnf4a-dependent Cebpa expression. Overall, we demonstrate how time-resolved transcriptomics can derive mechanistic hypothesis and how this can be combined with other streams of causal evidence.

Staphylococcus aureus is a high-priority pathogen causing severe infections with high morbidity and mortality worldwide. It is one of the most successful and prominent modern pathogens. An effective fight against S. aureus infections requires novel antimicrobial targets. Recent advances in whole-genome sequencing and high-throughput techniques facilitate the generation of genome-scale metabolic models (GEMs). Among the multiple applications of GEMs is drug-targeting in pathogens. Hence, comprehensive and predictive metabolic reconstructions of S. aureus could facilitate the identification of novel targets for antimicrobial therapies.This work aims at reviewing all 114 available GEMs of multiple S. aureus strains. We updated each model to a current version of SBML and evaluated its scope, including the number of reactions, metabolites, and genes. Furthermore, all models were quality-controlled using Mᴇᴍᴏᴛᴇ. Growth capabilities and model similarities were examined.This work should lead as a guide for choosing the appropriate GEM for a given research question. With the information about the availability, the format, and the strengths and potentials of each model, one can either choose an existing model or combine several models to create models with even higher predictive values. This facilitates model-driven discoveries of novel antimicrobial targets to fight multi-drug resistant S. aureus strains.

CD4+ T cells provide cell-mediated protection against diseases. When dysregulated, CD4+ T cells are associated with autoimmune and other immune-mediated diseases. Metabolism of CD4+ T cells regulates their function, therefore offer an opportunity to explore as a drug target against autoimmune diseases. In this study, we developed constraint-based models of naive and T helper 1, 2, and 17 subtypes. We mapped existing drugs and compounds and simulated metabolic behaviors under drug-induced inhibitions of metabolic genes. We integrated these metabolic behaviors with gene expression data of three autoimmune diseases, rheumatoid arthritis (RA), multiple sclerosis (MS), and primary biliary cholangitis (PBC). We identified and prioritized drugs and their targets that reversed the directions of differentially expressed genes in diseases. We identified 68 metabolic drug targets for the three studied diseases. We performed in vitro experiments and mined experiments available in the literature to validate results. The experimental results showed that 50% of the drug targets suppressed CD4+ T cell proliferation. In the end, we developed an integrated pipeline to explore metabolic models to identify drug targets and repurposable drugs.

We present a statistical learning framework to infer deterministic mean-field models from stochastic dynamics of mechanically active agents, such as cytoskeletal filaments and motor proteins, cells, or animals in a swarm. Physical theories of self-organized, non-equilibrium active systems exist both at the microscopic scale of individual agents (e.g., Brownian dynamics, Langevin equations), and at the continuum mean-field scale (e.g., active polar gel equtions). Naturally, the question arises how interactions between microscopic components lead to mean-field dynamics and which continuum description is sufficient to describe a given microscopic collective. Deriving such a coarse-graining by hand using classic methods like homogenization or volume averaging is tedious and difficult, due to the non-equilibrium physics, nonlinearities, and inherent fluctuations of the microscopic system. Data-driven inference of coarse-grained models from microscopic data has therefore emerged as a complementary approach, enabled by advances in machine learning (Rudy et al., 2017).

Single-cell RNA-Sequencing has made possible to infer high-resolution gene regulatory networks (GRNs), providing deep biological insights by revealing regulatory interactions at single-cell resolution. Current single-cell GRN inference methods produce a single GRN per input dataset, limiting their capability to uncover complex regulatory relationships across related cell phenotypes and potentially missing relationships between cells from different phenotypes.
We present SimiC, a single-cell GRN inference framework that enables uncovering complex regulatory dynamics at single-cell resolution across a range of systems, both model and non-model alike, that would have otherwise been missed by previously proposed methods. SimiC produces one GRN per phenotype while imposing a similarity constraint that allows for a smooth transition between GRNs, enabling a direct comparison between different states, treatments, or conditions.
We tested SimiC on simulated and real datasets, showing improved performance with respect to previous methods. Specifically, SimiC was able to: i) recapitulate CAR T cell dynamics after tumor recognition; ii) unravel well-known regulatory patterns on a regenerating liver; and iii) implicating glial cells in the generation of distinct behavioral states in honeybees.
Thus, SimiC establishes a new approach to quantitating regulatory architectures between GRNs of distinct cellular phenotypes, with far-reaching implications for systems biology.

Summary: Studying biological systems generally relies on computational modeling and simulation, e.g., for model-driven discovery and hypothesis testing. Progress in standardization efforts led to the development of interrelated file formats to exchange and reuse models in systems biology, such as SBML, the Simulation Experiment Description Markup Language (SED-ML), or the Open Modeling EXchange format (OMEX). Conducting simulation experiments based on these formats requires efficient and reusable implementations to make them accessible to the broader scientific community and to ensure the reproducibility of the results. The Systems Biology Simulation Core Library (SBSCL) provides interpreters and solvers for these standards as a versatile open-source API in Java™. The library simulates even complex biomodels and supports deterministic Ordinary Differential Equations (ODEs); Stochastic Differential Equations (SDEs); constraint-based analyses; recent SBML and SED-ML versions; exchange of results, and visualization of in silico experiments; open modeling exchange formats (COMBINE archives); hierarchically structured models; and compatibility with standard testing systems, including the Systems Biology Test Suite and published models from the BioModels and BiGG databases.Availability: SBSCL is freely available at https://draeger-lab.github.io/SBSCL/ and via Maven.

Effector-Triggered Immunity (ETI), which is induced by some pathogen effectors, is one major mode of plant innate immunity and is associated with massive transcriptional reprogramming. Plant cells directly receiving such effectors (cell population 1) and those indirectly responding (cell population 2) probably have different transcriptome dynamics but cannot be characterized separately with bulk RNA-seq. We analyzed a public time-series transcriptome dataset profiled during ETI in Arabidopsis after inoculation of Pseudomonas syringae strain expressing the effector AvrRpt2. We found that double-peak transcript patterns were prevalent among 2329 ETI-inducible genes. To quantitatively interpret the double-peak patterns, we developed a novel computational approach based on Multi-compartment models (MCMs) and fit the model to the transcript data of each gene. The model enables a linear decomposition of a double-peak pattern into two single-peak patterns and therefore allows us to characterize the peak patterns parametrically, such as the peak time and peak amplitude. We demonstrated that the first- and second-peaks of most genes highly likely represent the responses in cell populations 1 and 2, respectively and the timing and amplitude of the second peak response for each gene echoes the dynamics of the first peak response for the majority of genes.

Abstract: Clonal haematopoiesis (CH) describes the clonal expansion of mutant haematopoietic stem cells in healthy and mostly elderly individuals. Yet CH is also considered an important step in the progression to blood cancers, such as acute myeloid leukaemia and myelodysplastic syndrome. Therefore, understanding CH evolution is important to better understand the progression from seemingly inoffensive somatic mutation to cancer. Here, we use longitudinal deep targeted sequencing data for 385 elderly individuals to estimate clonal dynamics and age at onset for all clones through hierarchical Bayesian modelling. We observe distinct annual growth rates for each gene, with the lowest occurring in DNMT3A (4%) and the highest in U2AF1 (21%). The extrapolated age at onset of clonal growth shows similar patterns for most genes, with clones arising uniformly through life. Against this trend, mutations such as those in U2AF1 or SRSF2-P95H appear exclusively later in life. With single cell colonies for 3 individuals we are able to verify our estimates and observe clonal expansions with no known drivers, suggesting that a portion of CH goes by unnoticed. This work offers a comprehensive view on the dynamics and evolution of CH, disentangling the effects that driver genetics and other effects have on CH progression.

The phylum Apicomplexa comprises various medically-relevant parasites with diverse life cycle strategies, driven at least partly by distinct metabolic capacities. While the malaria-causing Plasmodium falciparum is mosquito-borne and infects liver cells and erythrocytes, Toxoplasma gondii is able to infect any nucleated cell of warm-blooded animals. Through constraint-based reconstruction and analysis, we seek to understand how each parasite has adapted its metabolism to best exploit its corresponding host environment. Our central hypothesis is that each parasite has reduced its metabolic potential in response to its ability to scavenge its hosts’ nutrients. We begin by reconstructing the metabolic models of various Apicomplexan parasites together with their free-living relatives (Vitrella brassicaformis and Chromera velia). We then apply mixed-integer linear programming to predict minimal metabolic models from the free-living models—functional models with the smallest number of reactions. Following in silico knock-out experiments in the free-living models, we identify essential reactions and compute their degree of conservation in the parasitic models. Next, we perform simulations to investigate which metabolic pathways are favoured in minimal models given different profiles of nutrient availabilities (representing different host environments). Overall, our analyses suggest that Apicomplexan metabolism has evolved reductively and so as to most efficiently exploit host environments.

The underlying dynamics of biological signalling networks is often unknown. Previous work on unravelling these complex networks includes heuristic and exact methods that aim to uncover the Boolean functions behind each interaction in an annotated network by explaining experimental observations given as stimuli/response sets.Our main goal is to infer the logical topology of a network, without the need for a prior annotated network. We consider a generalized framework for inferring Boolean functions of increasing logical complexity, including rudimentary Boolean functions, threshold gates and completely unrestricted truth tables. As a secondary objective, we compare integer linear programming (ILP) to satisfiability solving (SAT) solutions. For each complexity level, we provide two alternative formulations, one based on ILP, which extends previous work to the generalized setting, and a novel formulation based on SAT. Preliminary results of experiments on a well-known EGFR dataset as well as on synthetic data indicate that the SAT implementation significantly outperforms the ILP implementation in terms of running time and can infer an optimal model in networks of up to 205 nodes and 204 edges given 18 experiments in less than 5 minutes on a standard laptop computer.

Overstimulation of catabolic events compared to biosynthetic activity of cartilage chondrocytes (CC) leads to articular cartilage degradation when osteoarthritis (OA) appears. Mechanical and biochemical signals are involved in this dysregulation. To allow the semi-quantitative interpretation of these proceedings, we propose a network-based model (NBM) at the CC level that incorporates the actions of catabolic and anabolic factors on the expression of structural proteins and proteases. Understanding better such mechanisms might leverage the identification of new therapeutic targets that could stop OA progression and improve current conservative pain-killers’ treatments. Our NBM appears from a combination of knowledge-based and a data-driven approaches. First, an OA specific protein interactome was developed based on reported actions in specialized literature. Then, a targeted network enrichment is developed using STRING. Finally, it is successfully calibrated against experimental data through a genetic algorithm. Independent validation with multiplexed phosphoproteomics measurements shows that the optimized NBM has a relative error of 3.5%. It also captures CC’s reported behaviors with 95% of accuracy, and it correctly predicts the main outcomes of an OA treatment with biologics. Therefore, the proposed methodology allows us to model an optimal NBM that controls CC metabolism. Further research should target the incorporation of mechanical signals.

Through mechanotransduction (MT) mechanical loads cause changes on chondrocyte (CC) phenotype. Systematic representations of these procedures can inform early strategies for osteoarthritis management, characterized by an increase of catabolic events compared to biosynthetic activity. To this end, a network model of CC activity that incorporates important mechanosensors (MS) and the downstream signalling cascades is proposed. Key mechanoregulatory molecular interactions is mapped in an interactome, based on 72 journal articles. It is translated into a semi-quantitative mathematical model based on a system of differential equations that converges to steady states (SS), which can be extrapolated to systematic descriptions of CC metabolism. An anabolic SS is reached when a MS related to physio-osmotic conditions is activated: TRPV4. A pro-inflammatory stimulus and an injurious load perturbation (PIEZO1/2 activation) can induce a catabolic shift (t-test, α=0.05) compared to the physio-osmotic SS: pro-inflammatory cytokines and proteases have higher expression rates than structural proteins. Remarkably, Sox9 (a healthy marker), is highly expressed in the physio-osmotic SS; but in an injurious environment NFkB, HIF2a or Runx2 increase (related to inflammatory and hypertrophic events). An intracellular network-based model of a CC is developed that could predict expected MT and inflammation effects revealing the potential of exploitation in OA.

Our group developed a whole-cell protein complex prediction tool called Cytocast. With that tool, we can simulate how the proteins bind and unbind and form different kinds of protein complexes, which can help us to understand how the cell behaves in different conditions. Based on the input data, our model can simulate any cell in any state. The input file describes the cell’s geometrical structure, the proteins and complexes, and localization and abundances. The interactions are described on the level of the proteins’ binding sites. These binding sites do not have any structure, and they are referred to as point-like objects. The cell is divided into tens of thousands of sub-volumes. Thus the reaction simulation can run separately, and the particles are mixed through the diffusions.
We created a pipeline that provides data from our database and analysis tools to evaluate the simulation results of the different conditions. We can use this to compare treated and not-treated cells for different tissues and predicting the main and side effects of drugs.

The prediction of emergent variant characteristics remains extremely challenging. Faced with a new pathogen, it is desirable to be able to forecast the case fatality rate (CFR) into the future. Leveraging a compartment model, we reveal general constraints among human pathogenic respiratory viruses where the variation of multiple parameters in concert leads to decreased virulence and increased pathogen fitness but not independent variation. Highly virulent viruses are constrained by host behavior, whereas moderately virulent viruses are constrained by the relationship between the duration of immunity and CFR. When the immune population is rapidly expanded through vaccination, specific predictions can be made in the face of dramatically altered selective pressures. The potential emergence of vaccine resistance constrains optimal vaccine distribution. Analogous to low-dose antibiotic exposure, recently vaccinated, partially immunized individuals play an outsized role in the emergence of resistance. When an escape variant is modestly less infectious than the originating strain, there exists an optimal rate of vaccine distribution. Exceeding this rate increases the cumulative number of infections due to vaccine escape. Modulating the rate of host-host contact for the recently vaccinated population by less than an order of magnitude can alter the cumulative number of infections by more than 20%.

Kinetic modeling is essential in understanding the dynamic behavior of biochemical networks, such as metabolic and signal transduction pathways. However, parameter estimation remains a major bottleneck in the development of kinetic models.We present RCGAToolbox, software for real-coded genetic algorithms (RCGAs), which accelerates the parameter estimation of kinetic models. RCGAToolbox provides two RCGAs: the unimodal normal distribution crossover with minimal generation gap (UNDX/MGG) and real-coded ensemble crossover star with just generation gap (REXstar/JGG), using the stochastic ranking method. The RCGAToolbox also provides user-friendly graphical user interfaces.We tested 33 problems: 27 mathematical benchmarks and six parameter estimations. The performance of RCGAToolbox was comparable to that of an existing parameter estimation tool. We also found that REXstar/JGG outperformed a widely-used algorithm in one of the parameter estimation problems.RCGAToolbox is available from https://github.com/kmaeda16/RCGAToolbox under GNU GPLv3, with application examples. RCGAToolbox runs on MATLAB in Windows, Linux, and Mac.

In the tumor microenvironment, tumor-associated macrophages are known to play a critical role in the survival and chemoresistance of cancer cells. In the case of chronic lymphocytic leukemia (CLL), these tumor-associated macrophages are called Nurse-Like Cells (NLCs) and reside mainly in the lymph nodes, where they are able to protect leukemic B cells (B-CLL) from spontaneous apoptosis and contribute to their chemoresistance. NLCs are differentiated from monocytes through cytokines signalling and physical contact with the cancer cells [1], however, the precise mechanisms by which B-CLL cells influence this differentiation are still unknown.We propose here an agent-based model (ABM) of monocyte-to-macrophage differentiation in an in-vitro co-culture of monocytes and cancer B-CLL cells. This model is a first step to a better understanding of the spatio-temporal dynamics of the tumor microenvironment and of the mechanisms used by cancer cells to influence monocyte differentiation into pro-tumoral macrophages. A more complex model taking into account other immune cell types and integrating Boolean modeling of relevant signaling pathways inside each agent is the next step.[1] F Boissard et al. Nurse like cells: Chronic lymphocytic leukemia associated macrophages. Leuk. Lymphoma, 56(5):1570- 1572, 2015.

Recent rapid advances in genome synthesis have increased the significance of genome design techniques. Whole-cell modeling is one of the effective methods for genome design. However, whole-cell models are often constructed manually, and it is difficult to keep them up-to-date against the huge amount of new knowledge. In addition, it is now required to build models on-demand for single nucleotide-level design. Here, we developed a technique to automatically construct and simulate whole-cell models of a Gammaproteobacteria represented by Escherichia coli, based on a genome sequence and multiple omics data as its input. This workflow includes predictions of annotations, such as CDSs, operons, and replication initiation sites, on the given sequence, kinetic modeling of metabolic pathways, and parameter determination against multi-omics data (transcriptome, proteome, and metabolome). The workflow is constructed on Snakemake and confers the reproducibility of modeling. The whole-cell model consists of gene expression and replication system described in an agent-based manner with more than 1 million agents, and a kinetic metabolic model governed by ordinary differential equations of hundreds of variables. The automated modeling can facilitate the construction and maintenance of large complex models and allow us to design a genome from scratch toward an era of synthetic cells.

This brief talk will welcome the audience and speakers to the second day of the annual SysMod meeting in 2021 and announce the schedule of today’s meeting. There will be two sessions on integrative approaches and methodologies, followed by the final keynote talk by Boris N. Kholodenko about structure-based dynamic modeling, as well as Juilee Thakar’s announcement of this year’s poster awards.

Simple mathematical models have had remarkable successes in biology, framing how we understand a host of mechanisms and processes. However, with the advent of a host of new experimental technologies, the last ten years has seen an explosion in the amount and types of quantitative data now being generated. This sets a new challenge for the field – to develop, calibrate and analyse new, biologically realistic models to interpret these data. In this talk I will showcase how quantitative comparisons between models and data can help tease apart subtle details of biological mechanisms, as well as present some steps we have taken to tackle the mathematical challenges in developing models that are both identifiable and can be efficiently calibrated to quantitative data.

COVID-19 has been characterized as one of the deadliest respiratory diseases. In this regard, scientists globally try to understand the host’s immunopathological response, how the novel coronavirus (SARS-CoV-2) adapts, and how it spreads. Currently, great efforts are made to detect effective therapies against the coronaviruses. Identifying potential antiviral targets is of great interest, and one way to detect them is by analyzing metabolic changes in infected cells. In 2012, Wang et al. published mCADRE, aiming to reconstruct tissue-specific models using gene expression data and network topology information. The algorithm is implemented in MATLAB, and its functionality is based solely on the first version of the human model resulting in its limited usability.We present pymCADRE, a re-implementation of mCADRE in Python 3.8. Its functionality was tested using all three currently available versions of the human metabolic network. Internal optimizations done with fastFVA resulted in context-specific models closer to the ground truth. Additionally, host-virus models were created to help to identify potential antiviral targets against SARS-CoV-2. With those models, the recently identified potential antiviral target enzyme guanylate kinase was further investigated. Further improvements could be done to make it feasible with more complex models, like Recon2.2 and Recon3D.

Cell classifiers are synthetic bio-devices performing type-specific, in vivo classification of the cell’s molecular fingerprint. In particular, they can recognize cancerous cells and trigger their apoptosis, shaping novel therapies for cancer patients. Although a single circuit’s processing logic is usually described using a Boolean function, other architectures have also been considered, e.g. multi-circuit designs. Distributed classifiers consist of a group of single-circuit classifiers deciding collectively according to a predefined threshold function whether a cell is cancerous. Such architecture has shown the potential to predict the cell condition with high accuracy. However, the lack of far-reaching machinery to design and evaluate the classifiers, in particular, assessing their robustness to noise and novel information, makes their application limited. Here, we present a framework for designing miRNA-based distributed cell classifiers combining genetic algorithms and logic programming. We develop optimization criteria comprising the accuracy and robustness of the circuits that allow achieving high performance as shown in multiple simulated data studies. The evaluation performed on real cancer data demonstrates that distributed classifiers outperform single-circuit designs. Our classifiers include relevant miRNAs previously described in the literature, as well as more complex regulation patterns included in the data.

The Gillespie algorithm is commonly applied for simulating memoryless processes that follow an exponential waiting-time. However, stochastic processes governing biological interactions, such as cell apoptosis and epidemic spreading, are empirically known to exhibit properties of memory, an inherently non-Markovian feature. The presence of such non-Markovian processes can significantly influence the outcome of a simulation, While several extensions to the Gillespie algorithm have been proposed, most of them suffer from either a high computational cost or are only applicable to a narrow selection of probability distributions that do not match the experimentally observed biological data distributions. To tackle the aforementioned issues, we developed a Rejection Gillespie for non-Markovian Reactions (REGINR) that is capable of generating simulations with non-exponential waiting-times, while remaining an order of magnitude faster than alternative approaches. REGINR uses the Weibull distribution, which interpolates between the exponential, normal, and heavy-tailed distributions. We applied our algorithm to a mouse stem cell dataset with known non-Markovian dynamics and found it to faithfully recapitulate the underlying biological processes. We conclude that our algorithm is suitable for gaining insight into the role of molecular memory in stochastic models, as well as for accurately simulating real-world biological processes.

Glioblastoma multiforme (GBM) is a highly aggressive form of brain cancer that has a 5-year survival rate of about 5%. The low survival rate can be at least partially attributed to its heterogeneity. The presence of multiple genetically distinct clones within the solid tumor, as well as its stem cell nature, leads to treatment resistance. Predicting treatment resistance or susceptibility is difficult, even when the tumor has been classified according to genetic subtype. These classifications rely on few biomarkers and are informed by limited mechanistic details. Here we present a discrete computational model of GBM stem cell dynamics, which has been informed by expert knowledge, biomedical literature, and biological data. We show that the interactions within the model are well-supported by both literature and database sources. Our approach for parameterizing the model includes fitting to biological data such as gene expression to better predict individual cell line responses to treatment. We show that our GBM stem cell model is capable of predicting the success of a high percentage of kinase inhibitor treatments, and that these predictions are cell line specific. Our computational model offers a rapid, individualized approach for predicting drug treatment efficacy across GBM stem cell lines.

Understanding the mechanical strategies by which cancer cells migrate within confined spaces is essential to raise novel insights regarding metastasis prevention. We explored a phase-field model where cells are described as droplets with surface tension that interact with the extracellular matrix by adhesion and excluded volume. With the purpose of verifying that our results are model independent, we produced an equivalent system using a Dissipative Particle Dynamics (DPD) approach. We observed that adhesion regulates cell deformability and enhances migration under extreme confinement conditions for both models. Furthermore, we were able to reproduce different experimental conditions, both with and without matrix metalloproteinases (MMPs) inhibition, solely by changing the adhesion coefficient in our models. This might indicate that MMPs act not only on degradation but also play a key role in cell migration and adhesion.

The study of metabolism and its regulation is finding increasing application in various fields, including health, wellness, and biotransformations. Complete characterization of regulatory mechanisms controlling metabolism requires knowledge of metabolic fluxes, whose direct determination lags behind other omic technologies, such as metabolomics and transcriptomics. In isolation, these methodologies do not allow accurate characterization of metabolic regulation. Hence, there is a need for integrated methodologies to disassemble the interdependence between different regulatory layers controlling metabolism.To this aim, we propose a computational pipeline to characterize the landscape of metabolic regulation in different biological samples. The method integrates intracellular and extracellular metabolomics, and transcriptomics, using constraint-based stoichiometric metabolic models as a scaffold. We compute differential reaction expression from transcriptomic data and use constraint-based modeling to predict if the differential expression of metabolic enzymes directly originates differences in metabolic fluxes. In parallel, using metabolomic data, we predict how differences in substrate availability translate into differences in metabolic fluxes. By intersecting these two output datasets, we discriminate fluxes regulated at the metabolic and/or transcriptional level. This information is valuable to better inform targeted action planning in different fields, including personalized prescriptions in multifactorial diseases, such as cancer, and metabolic engineering.

The eukaryotic cell cycle is driven by waves of cyclin-dependent kinase (cyclin/Cdk) activities that rise and fall with a timely pattern called “waves of cyclins”. This pattern guarantees coordination and alternation of DNA synthesis with cell division, and its failure results in altered cyclin/Cdk dynamics and abnormal cell proliferation. Although details about transcription of cyclins are available, the network motifs responsible for this timely pattern are currently unknown. Here we reveal a novel principle of design that ensures cell cycle time keeping through interlocking transcription with cyclin/Cdk dynamics in budding yeast. Deterministic and stochastic analyses of 1024 kinetic models of the cyclin/Cdk network are verified against quantitative data of Clb dynamics. A novel regulatory design is unravelled, which involves the evolutionarily conserved Forkhead (Fkh) transcription factors. The Fkh-mediated cascade among Clb cyclins and Clb/Cdk1-mediated positive feedback loops are pivotal for synchronizing Clb/Cdk1 waves. Furthermore, our model predicts a definite Fkh activation pattern underlying this design, with a progressive Clb/Cdk1-mediated Fkh phosphorylation. Experimental validation confirms the computational prediction, highlighting the Clb/Cdk–Fkh axis being pivotal for timely cell cycle dynamics. This work rationalizes the quantitative model of Cdk control for budding yeast, identifying the regulatory motifs underlying cell proliferation dynamics.

This is a brief discussion on the two sessions on Integrative approaches and methodologies and an outlook to the final session about structure-based dynamic modeling and SysMod poster award 2021.

Major problem encountered using small molecule cancer therapeutics in clinic is that even in susceptible cancers, these drugs rarely give durable responses, almost inevitably being hampered by signaling reactivation and development of resistance. Studying the causes of this resistance has revealed severe limitations in our understanding of the network properties and molecular mechanisms that control drug responses. We show that contrary to a common opinion, feedback loops by themselves cannot restore or overshoot steady state signaling. De novo synthesized negative feedback regulators can lead to a transient overshoot but still cannot fully restore output signaling. These findings can rationalize recent scientific and clinical disappointments that were based on the hypothesis that negative feedback loops can fully explain drug resistance. We demonstrate that there are two major means of complete, steady state revival of signaling, enabled by (1) the network topology or (2) molecular mechanisms rendering the primary drug target active again. Network topology analysis shows that at least two, activating and inhibitory, connection routes from a primary drug target to the output, must exist for complete reactivation or overshoot of steady-state output activity that existed before the inhibition.
Irrespective of the network topology, drug-induced overexpression of the primary drug target or drug-induced increase in its dimerization or oligomerization can restore the pathway output activity. The formation of kinase homo- or heterodimers is a major course of resistance. In this constellation one protomer is drug-bound and allosterically activates the other, drug-free protomer thereby conferring resistance. The emergence of different drug affinities between protomers in a dimer has been enigmatic, but can be explained by thermodynamics (https://www.ncbi.nlm.nih.gov/pubmed/26344764). A striking example is so-called paradoxical activation of the extracellular regulated kinase (ERK) pathway by RAF inhibitors, which is caused by RAF homo- or heterodimerization. This dimerization is promoted by RAF inhibitors and amplified by mutant RAS and negative feedback regulations, but if an inhibitor does not facilitate dimerization, negative feedback can only result in a transient overshoot of the pathway activity. Exciting and counterintuitive discoveries of ways to overcome resistance were made using next generation modelling, which combines aspects of protein structure, posttranslational modifications, thermodynamics, network architecture, mutation data and dynamic reaction mechanisms (https://www.ncbi.nlm.nih.gov/pubmed/30007540). As a specific example, we show that a treatment with Type I½ and Type II RAF inhibitors can counterbalance ERK pathway reactivation and concomitant drug resistance.

This talk briefly reviews the two days of SysMod, including all the speakers, chairpersons, and organizers. The first day comprised three sessions on disease and multi-scale modeling, with one session on infectious disease modeling in particular. In addition, two keynote talks were given by Ines Thile and Ruth E. Baker. The second day included two sessions on integrative approaches and methodologies and one session on Structure-based dynamic modeling, with the final keynote just given by Boris N. Kholodenko. Next, it is time to thank all contributions to scientific posters and bestow the best ones with the annual SysMod poster awards in 2021. Finally, we will conclude the meeting and open the discussion with the audience for feedback to the next SysMod in 2022.