Rheumatoid Arthritis (RA) is a chronic autoimmune inflammatory disease that affects about 0.1% to 2% of the population worldwide. Despite the development of several novel therapies, there is only limited benefit for many patients. Thus, there is room for new approaches to improve response to therapy, including designing better trials e.g., by identifying subpopulations that can benefit from specific classes of therapy and enabling reverse translation by analyzing completed clinical trials. We have developed an open-source, mechanistic multi-scale model of RA, which captures the interactions of key immune cells and mediators in an inflamed joint. The model consists of a treatment-naive Virtual Population (Vpop) that responds appropriately (i.e. as reported in clinical trials) to standard-of-care treatment options-Methotrexate (MTX) and Adalimumab (ADA, anti-TNF-α) and an MTX inadequate responder sub-population that responds appropriately to Tocilizumab (TCZ, anti-IL-6R) therapy. The clinical read-outs of interest are the American College of Rheumatology score (ACR score) and Disease Activity Score (DAS28-CRP), which is modeled to be dependent on the physiological variables in the model. Further, we have validated the Vpop by predicting the therapy response of TCZ on ADA Non-responders. This paper aims to share our approach, equations, and code to enable community evaluation and greater adoption of mechanistic models in drug development for autoimmune diseases.

Understanding the dynamic states and transitions of heterogeneous cell populations is crucial for addressing fundamental biological questions. High-content imaging provides rich datasets, but it remains increasingly difficult to integrate and annotate high-dimensional and time-resolved datasets to profile heterogeneous cell population dynamics in different microenvironments. Using hepatic stellate cells (HSCs) LX-2 as model, we proposed a novel analytical strategy for image-based integration and annotation to profile dynamics of heterogeneous cell populations in 2D/3D microenvironments. High-dimensional features were extracted from extensive image datasets, and cellular states were identified based on feature profiles. Time-series clustering revealed distinct temporal patterns of cell shape and actin cytoskeleton reorganization. We found LX-2 showed more complex membrane dynamics and contractile systems with an M-shaped actin compactness trend in 3D culture, while they displayed rapid spreading in early 2D culture. This image-based integration and annotation strategy enhances our understanding of HSCs heterogeneity and dynamics in complex extracellular microenvironments.

Amyotrophic Lateral Sclerosis (ALS) is a devastating, immensely complex neurodegenerative disease by lack of effective treatments. We developed a network medicine methodology via integrating human brain multi-omics data to prioritize drug targets and repurposable treatments for ALS. We leveraged non-coding ALS loci effects from genome-wide associated studies (GWAS) on human brain expression quantitative trait loci (QTL) (eQTL), protein QTL (pQTL), splicing QTL (sQTL), methylation QTL (meQTL), and histone acetylation QTL (haQTL). Using a network-based deep learning framework, we identified 105 putative ALS-associated genes enriched in known ALS pathobiological pathways. Applying network proximity analysis of predicted ALS-associated genes and drug-target networks under the human protein-protein interactome (PPI) model, we identified potential repurposable drugs (i.e., Diazoxide and Gefitinib) for ALS. Subsequent validation established preclinical evidence for top-prioritized drugs. In summary, we presented a network-based multi-omics framework to identify drug targets and repurposable treatments for ALS and other neurodegenerative disease if broadly applied.

Petri nets are commonly applied in modeling biological systems. However, construction of a Petri net model for complex biological systems is often time consuming, and requires expertise in the research area, limiting their application. To address this challenge, we developed GINtoSPN, an R package that automates the conversion of multi-omics molecular interaction network extracted from the Global Integrative Network (GIN) into Petri nets in GraphML format. These GraphML files can be directly used for Signaling Petri Net (SPN) simulation. To demonstrate the utility of this tool, we built a Petri net model for neurofibromatosis type I. Simulation of NF1 gene knockout, compared to normal skin fibroblast cells, revealed persistent accumulation of Ras-GTPs as expected. Additionally, we identified several other genes substantially affected by the loss of NF1's function, exhibiting individual-specific variability. These results highlight the effectiveness of GINtoSPN in streamlining the modeling and simulation of complex biological systems.

Since nail psoriasis restricts the patient's daily activities, therapeutic intervention based on reliable and reproducible evaluation is critical. The Nail Psoriasis Severity Index (NAPSI) is a validated scoring tool, but its usefulness is limited by interobserver variability. This study aimed to develop a reliable and accurate NAPSI scoring tool using deep learning. The tool "NAPSI calculator" includes two parts: nail detection from images and NAPSI scoring. NAPSI was annotated by nine nail experts who are board-certified dermatologists with sufficient experience in a specialized clinic for nail diseases. In the final test set, the "NAPSI calculator" correctly located 137/138 nails and scored NAPSI with higher accuracy than the compared six non-board-certified residents: 83.9% vs 65.7%; P = 0.008 and four board-certified non-nail expert dermatologists: 83.9% vs 73.0%; P = 0.005. The "NAPSI calculator" can be readily used in a clinical situation, contributing to raising the medical practice level for nail psoriasis.

Normal epithelial cells receive proper signals for growth and survival from attachment to the underlying extracellular matrix (ECM). They perceive detachment from the ECM as a stress and die - a phenomenon termed as 'anoikis'. However, metastatic cancer cells acquire anoikis-resistance and circulate through the blood and lymphatics to seed metastasis. Under normal (adherent) growth conditions, the serine-threonine protein kinase Akt stimulates protein synthesis and cell growth, maintaining an anabolic state in the cancer cell. In contrast, previously we showed that the stress due to matrix deprivation is sensed by yet another serine-threonine kinase, AMP-activated protein kinase (AMPK), that inhibits anabolic pathways while promoting catabolic processes. We illustrated a switch from Akt/AMPK in adherent condition to AMPK/Akt in matrix-detached condition, with consequent metabolic switching from an anabolic to a catabolic state, which aids cancer cell stress-survival. In this study, we utilized these experimental data and developed a deterministic ordinary differential equation (ODE)-based mechanistic mathematical model to mimic attachment-detachment signaling network. To do so, we used the framework of insulin-glucagon signaling with consequent metabolic shifts to capture the pathophysiology of matrix-deprived state in breast cancer cells. Using the developed metastatic breast cancer signaling (MBCS) model, we identified perturbation of several signaling proteins such as IRS, PI3K, PKC, GLUT1, IP3, DAG, PKA, cAMP, and PDE3 upon matrix deprivation. Further, in silico molecular perturbations revealed that several feedback/crosstalks like DAG to PKC, PKC to IRS, S6K1 to IRS, cAMP to PKA, and AMPK to Akt are essential for the metabolic switching in matrix-deprived cancer cells. AMPK knockdown simulations identified a crucial role for AMPK in maintaining these adaptive changes. Thus, this mathematical framework provides insights on attachment-detachment signaling with metabolic adaptations that promote cancer metastasis.

Intensive malaria control and elimination efforts have led to substantial reductions in malaria incidence over the past two decades. However, the reduction in Plasmodium falciparum malaria cases has led to a species shift in some geographic areas, with P. vivax predominating in many areas outside of Africa. Despite its wide geographic distribution, P. vivax vaccine development has lagged far behind that for P. falciparum, in part due to the inability to cultivate P. vivax in vitro, hindering traditional approaches for antigen identification. In a prior study, we have used a positive-unlabeled random forest (PURF) machine learning approach to identify P. falciparum antigens based on features of known antigens for consideration in vaccine development efforts. Here we integrate systems data from P. falciparum (the better-studied species) to improve PURF models to predict potential P. vivax vaccine antigen candidates. We further show that inclusion of known antigens from the other species is critical for model performance, but the inclusion of only the unlabeled proteins from the other species can result in misdirection of the model toward predictors of species classification, rather than antigen identification. Beyond malaria, incorporating antigens from a closely related species may aid in vaccine development for emerging pathogens having few or no known antigens.

Phosphosignaling networks control cellular processes. We built kinase-mediated regulatory networks elicited by thrombin stimulation of brain endothelial cells using two computational strategies: Temporal Pathway Synthesizer (TPS), which uses phosphoproteomics data as input, and Temporally REsolved KInase Network Generation (TREKING), which uses kinase inhibitor screens. TPS and TREKING predicted overlapping barrier-regulatory kinases connected with unique network topology. Each strategy effectively describes regulatory signaling networks and is broadly applicable across biological systems.

Cellular networks realize their functions by integrating intricate information embedded within local structures such as regulatory paths and feedback loops. However, the precise mechanisms of how local topologies determine global network dynamics and induce bifurcations remain unidentified. A critical step in unraveling the integration is to identify the governing principles, which underlie the mechanisms of information flow. Here, we develop the cumulative linearized approximation (CLA) algorithm to address this issue. Based on perturbation analysis and network decomposition, we theoretically demonstrate how perturbations affect the equilibrium variations through the integration of all regulatory paths and how stability of the equilibria is determined by distinct feedback loops. Two illustrative examples, i.e., a three-variable bistable system and a more intricate epithelial-mesenchymal transition (EMT) network, are chosen to validate the feasibility of this approach. These results establish a solid foundation for understanding information flow across cellular networks, highlighting the critical roles of local topologies in determining global network dynamics and the emergence of bifurcations within these networks. This work introduces a novel framework for investigating the general relationship between local topologies and global dynamics of cellular networks under perturbations.

Biomarkers associated with the progression from gastric intestinal metaplasia (GIM) to gastric adenocarcinoma (GA), i.e., GA-related GIM, could provide valuable insights into identifying patients with increased risk for GA. The aim of this study was to utilize multi-bioinformatics to reveal potential biomarkers for the GA-related GIM and predict potential drug repurposing for GA prevention in patients. The multi-bioinformatics included gene expression matrix (GEM) by microarray gene expression (MGE), ScType (a fully automated and ultra-fast cell-type identification based solely on a given scRNA-seq data), Ingenuity Pathway Analysis, PageRank centrality, GO and MSigDB enrichments, Cytoscape, Human Protein Atlas and molecular docking analysis in combination with immunohistochemistry. To identify GA-related GIM, paired surgical biopsies were collected from 16 GIM-GA patients who underwent gastrectomy, yielding 64 samples (4 biopsies per stomach x 16 patients) for MGE. Co-analysis was performed by including scRNAseq and immunohistochemistry datasets of endoscopic biopsies of 37 patients. The results of the present study showed potential biomarkers for GA-related GIM, including GEM of individual patients, individual genes (such as RBP2 and CD44), signaling pathways, network of molecules, and network of signaling pathways with key topological nodes. Accordingly, potential treatment targets with repurposed drugs were identified including epidermal growth factor receptor, proto-oncogene tyrosine-protein kinase Src, paxillin, transcription factor Jun, breast cancer type 1 susceptibility protein, cellular tumor antigen p53, mouse double minute 2, and CD44.

Intratumoral delivery and localized chemotherapy have demonstrated promise in tumor treatment; however, the rapid drainage of therapeutic agents from well-vascularized tumors limits their ability to achieve maximum therapeutic efficacy. Therefore, innovative approaches are needed to enhance treatment efficacy in such tumors. This study utilizes a mathematical modeling platform to assess the efficacy of combination therapy using anti-angiogenic drugs and drug-loaded nanoparticles. Anti-angiogenic drugs are included to reduce blood microvascular density and facilitate drug retention in the extracellular space. In addition, incorporating negatively charged nanoparticles aims to enhance diffusion and distribution of therapeutic agents within well-vascularized tumors. The findings indicate that, in the case of direct injection of free drugs, using compounds with lower drainage rates and higher diffusion coefficients is beneficial for achieving broader diffusion. Otherwise, drugs tend to accumulate primarily around the injection site. For instance, the drug doxorubicin, known for its rapid drainage, requires the prior direct injection of an anti-angiogenic drug with a high diffusion rate to reduce microvascular density and facilitate broader distribution, enhancing penetration depth by 200%. Moreover, the results demonstrate that negatively charged nanoparticles effectively disperse throughout the tissue due to their high diffusion coefficient. In addition, a faster drug release rate from nanoparticles further enhance treatment efficacy, achieving the necessary concentration for complete eradication of tumor compared to slower drug release rates. This study demonstrates the potential of utilizing negatively charged nanoparticles loaded with chemotherapy drugs exhibiting high release rates for localized chemotherapy through intratumoral injection in well-vascularized tumors.

Endothelial cells are key players in the cardiovascular system. Among other things, they are responsible for sprouting angiogenesis, the process of new blood vessel formation essential for both health and disease. Endothelial cells are strongly regulated by the juxtacrine signaling pathway Notch. Recent studies have shown that both Notch and angiogenesis are influenced by extracellular matrix stiffness; however, the underlying mechanisms are poorly understood. Here, we addressed this challenge by combining computational models of Notch signaling and YAP/TAZ, stiffness- and cytoskeleton-regulated mechanotransducers whose activity inhibits both Dll4 (Notch ligand) and LFng (Notch-Dll4 binding modulator). Our simulations successfully mimicked previous experiments, indicating that this YAP/TAZ-Notch crosstalk elucidates the Notch and angiogenesis mechanoresponse to stiffness. Additional simulations also identified possible strategies to control Notch activity and sprouting angiogenesis via cytoskeletal manipulations or spatial patterns of alternating stiffnesses. Our study thus inspires new experimental avenues and provides a promising modeling framework for further investigations into the role of Notch, YAP/TAZ, and mechanics in determining endothelial cell behavior during angiogenesis and similar processes.

Cytokines operate in concert to maintain immune homeostasis and coordinate immune responses. In cases of ER breast cancer, peripheral immune cells exhibit altered responses to several cytokines, and these alterations are correlated strongly with patient outcomes. To develop a systems-level understanding of this dysregulation, we measured a panel of cytokine responses and receptor abundances in the peripheral blood of healthy controls and ER breast cancer patients across immune cell types. Using tensor factorization to model this multidimensional data, we found that breast cancer patients exhibited widespread alterations in response, including drastically reduced response to IL-10 and heightened basal levels of pSmad2/3 and pSTAT4. ER patients also featured upregulation of PD-L1, IL6Rα, and IL2Rα, among other receptors. Despite this, alterations in response to cytokines were not explained by changes in receptor abundances. Thus, tensor factorization helped to reveal a coordinated reprogramming of the immune system that was consistent across our cohort.

T cells are dynamically regulated immune cells that are implicated in a variety of diseases ranging from infection, cancer and autoimmunity. Recent advancements in sequencing methods have provided valuable insights in the transcriptional and epigenetic regulation of T cells in various disease settings. In this review, we identify the key sequencing-based methods that have been applied to understand the transcriptomic and epigenomic regulation of T cells in diseases.

Intratumor heterogeneity reduces treatment efficacy and complicates our understanding of tumor progression and there is a pressing need to understand the functions of heterogeneous tumor cell subpopulations within a tumor, yet systems to study these processes in vitro are limited. Single-cell RNA sequencing (scRNA-seq) has revealed that some cancer cell lines include distinct subpopulations. Here, we present clusterCleaver, a computational package that uses metrics of statistical distance to identify candidate surface markers maximally unique to transcriptomic subpopulations in scRNA-seq which may be used for FACS isolation. With clusterCleaver, ESAM and BST2/tetherin were experimentally validated as surface markers which identify and separate major transcriptomic subpopulations within MDA-MB-231 and MDA-MB-436 cells, respectively. clusterCleaver is a computationally efficient and experimentally validated workflow for identification of surface markers for tracking and isolating transcriptomically distinct subpopulations within cell lines. This tool paves the way for studies on coexisting cancer cell subpopulations in well-defined in vitro systems.

Antisense transcripts are a unique group of non-coding RNAs and play regulatory roles in a variety of biological processes, including circadian rhythms. Per2AS is an antisense transcript to the sense core clock gene Period2 (Per2) in mouse and its expression is rhythmic and antiphasic to Per2. To understand the impact of Per2AS-Per2 interaction, we developed a new mathematical model that mechanistically described the mutually repressive relationship between Per2 and Per2AS. This mutual repression can regulate both amplitude and period of circadian oscillation by affecting a negative feedback regulation of Per2. Simulations from this model also fit with experimental observations that could not be fully explained by our previous model. Our revised model can not only serve as a foundation to build more detailed models to better understand the impact of Per2AS-Per2 interaction in the future, but also be used to analyze other sense-antisense RNA pairs that mutually repress each other.

The cell cycle of budding yeast is governed by an intricate protein regulatory network whose dysregulation can lead to lethal mistakes or aberrant cell division cycles. In this work, we model this network in a Boolean framework for stochastic simulations. Our model is sufficiently detailed to account for the phenotypes of 40 mutant yeast strains (83% of the experimentally characterized strains that we simulated) and also to simulate an endoreplicating strain (multiple rounds of DNA synthesis without mitosis) and a strain that exhibits 'Cdc14 endocycles' (periodic transitions between metaphase and anaphase). Because our model successfully replicates the observed properties of both wild-type yeast cells and many mutant strains, it provides a reasonable, validated starting point for more comprehensive stochastic-Boolean models of cell cycle controls. Such models may provide a better understanding of cell cycle anomalies in budding yeast and ultimately in mammalian cells.

Tumor hypoxia is a dynamic phenomenon marked by fluctuations in oxygen levels across both rapid (seconds to minutes) and slow (hours to days) time scales. While short hypoxia cycles are relatively well understood, the mechanisms behind longer cycles remain largely unclear. In this paper, we present a novel mechanistic mathematical model that explains slow hypoxia cycles through feedback loops involving vascular expansion and regression, oxygen-regulated tumor growth, and toxic cytokine production. Our study reveals that, for the emergence of slow hypoxia cycles, endothelial cells must adapt by decreasing receptor activation as ligand concentration increases. Additionally, the interaction between tumor cells and toxic cytokines influences frequency and intensity of these cycles. By examining the effects of pharmacological interventions, specifically poly (ADP-ribose) polymerase inhibitors, we also demonstrate how targeting cell proliferation can help regulate oxygen levels. Our findings enhance the understanding of hypoxia regulation and suggest PARP proteins as promising therapeutic targets.

Elucidating the emergent dynamics of cellular differentiation networks is crucial to understanding cell-fate decisions. Toggle switch - a network of mutually repressive lineage-specific transcription factors A and B - enables two phenotypes from a common progenitor: (high A, low B) and (low A, high B). However, the dynamics of networks enabling differentiation of more than two phenotypes from a progenitor cell has not been well-studied. Here, we investigate the dynamics of a four-node network A, B, C, and D inhibiting each other, forming a toggle tetrahedron. Our simulations show that this network is multistable and predominantly allows for the co-existence of six hybrid phenotypes where two of the nodes are expressed relatively high as compared to the remaining two, for instance (high A, high B, low C, low D). Finally, we apply our results to understand naïve CD4 T cell differentiation into Th1, Th2, Th17 and Treg subsets, suggesting Th1/Th2/Th17/Treg decision-making to be a two-step process.

Genome-scale metabolic models (GEMs) cover the entire list of metabolic genes in an organism and associated reactions, in a tissue/condition non-specific manner. RNA-seq provides crucial information to make the GEMs condition-specific. Integrative Metabolic Analysis Tool (iMAT) and Integrative Network Inference for Tissues (INIT) are the two most popular algorithms to create condition-specific GEMs from human transcriptome data. The normalization method of choice for raw RNA-seq count data affects the model content produced by these algorithms and their predictive accuracy. However, a benchmark of the RNA-seq normalization methods on the performance of iMAT and INIT algorithms is missing in the literature. Another important phenomenon is covariates such as age and gender in a dataset, and they can affect the predictivity of analysis. In this study, we aimed to compare five different RNA-seq data normalization methods (TPM, FPKM, TMM, GeTMM, and RLE) and covariate adjusted versions of the normalized data by mapping them on a human GEM using the iMAT and INIT algorithms to generate personalized metabolic models. We used RNA-seq data for Alzheimer's disease (AD) and lung adenocarcinoma (LUAD) patients. The results demonstrated that RNA-seq data normalized by the RLE, TMM, or GeTMM methods enabled the production of condition-specific metabolic models with considerably low variability in terms of the number of active reactions compared to the within-sample normalization methods (FPKM, TPM). Using these models, we could more accurately capture the disease-associated genes (average accuracy of ~0.80 for AD and ~0.67 for LUAD) for the RLE, TMM, and GeTMM normalization methods. An increase in the accuracies was observed for all the methods when covariate adjustment was applied. We found a similar accuracy trend when we compared the metabolites of perturbed reactions to metabolome data for AD. Together, our benchmark study shows that the between-sample RNA-seq normalization methods reduce false positive predictions at the expense of missing some true positive genes when mapped on GEMs.

Chronic obstructive pulmonary disease (COPD) is an etiologically complex disease characterized by acute exacerbations and stable phases. We aimed to identify biological functions modulated in specific COPD conditions, using whole blood samples collected in the AERIS clinical study (NCT01360398). Considered conditions were exacerbation onset, severity of airway obstruction, and presence of respiratory pathogens in sputum samples. With an integrative multi-network gene community detection (MNGCD) approach, we analyzed expression profiles to identify communities of correlated genes. The approach combined different layers of gene interactions for each explored condition/subset of samples: gene expression similarity, protein-protein interactions, transcription factors, and microRNAs validated regulons. Heme metabolism, interferon-alpha, and interferon-gamma pathways were modulated in patients at both exacerbation and stable-state visits, but with the involvement of distinct sets of genes. An important gene community was enriched with G2M checkpoint, E2F targets, and mitotic spindle pathways during exacerbation. Targets of TAL1 regulator and hsa-let-7b - 5p microRNA were modulated with increasing severity of airway obstruction. Bacterial infections with Moraxella catarrhalis and, particularly, Haemophilus influenzae triggered a specific cellular and inflammatory response in acute exacerbations, indicating an active reaction of the host to infections. In conclusion, COPD is a complex multifactorial disease that requires in-depth investigations of its causes and features during its evolution and whole blood transcriptome profiling can contribute to capturing some relevant regulatory mechanisms associated with this disease. In this work, we explored multi-network modeling that integrated diverse layers of regulatory gene networks and enhanced our comprehension of the biological functions implicated in the COPD pathogenesis.