Longitudinal intervention studies with repeated measurements over time are an important type of experimental design in biomedical research. Due to the advent of "omics"-sciences (genomics, transcriptomics, proteomics, metabolomics), longitudinal studies generate increasingly multivariate outcome data. Analysis of such data must take both the longitudinal intervention structure and multivariate nature of the data into account. The ASCA+-framework combines general linear models with principal component analysis and can be used to separate and visualize the multivariate effect of different experimental factors. However, this methodology has not yet been developed for the more complex designs often found in longitudinal intervention studies, which may be unbalanced, involve randomized interventions, and have substantial missing data. Here we describe a new methodology, repeated measures ASCA+ (RM-ASCA+), and show how it can be used to model metabolic changes over time, and compare metabolic changes between groups, in both randomized and non-randomized intervention studies. Tools for both visualization and model validation are discussed. This approach can facilitate easier interpretation of data from longitudinal clinical trials with multivariate outcomes.

Controllability perception significantly influences motivated behavior and emotion and requires an estimation of one's influence on an environment. Previous studies have shown that an agent can infer controllability by observing contingency between one's own action and outcome if there are no other outcome-relevant agents in an environment. However, if there are multiple agents who can influence the outcome, estimation of one's genuine controllability requires exclusion of other agents' possible influence. Here, we first investigated a computational and neural mechanism of controllability inference in a multi-agent setting. Our novel multi-agent Bayesian controllability inference model showed that other people's action-outcome contingency information is integrated with one's own action-outcome contingency to infer controllability, which can be explained as a Bayesian inference. Model-based functional MRI analyses showed that multi-agent Bayesian controllability inference recruits the temporoparietal junction (TPJ) and striatum. Then, this inferred controllability information was leveraged to increase motivated behavior in the vmPFC. These results generalize the previously known role of the striatum and vmPFC in single-agent controllability to multi-agent controllability, and this generalized role requires the TPJ in addition to the striatum of single-agent controllability to integrate both self- and other-related information. Finally, we identified an innate positive bias toward the self during the multi-agent controllability inference, which facilitated behavioral adaptation under volatile controllability. Furthermore, low positive bias and high negative bias were associated with increased daily feelings of guilt. Our results provide a mechanism of how our sense of controllability fluctuates due to other people in our lives, which might be related to social learned helplessness and depression.

Nervous systems extract and process information from the environment to alter animal behavior and physiology. Despite progress in understanding how different stimuli are represented by changes in neuronal activity, less is known about how they affect broader neural network properties. We developed a framework for using graph-theoretic features of neural network activity to predict ecologically relevant stimulus properties, in particular stimulus identity. We used the transparent nematode, Caenorhabditis elegans, with its small nervous system to define neural network features associated with various chemosensory stimuli. We first immobilized animals using a microfluidic device and exposed their noses to chemical stimuli while monitoring changes in neural activity of more than 50 neurons in the head region. We found that graph-theoretic features, which capture patterns of interactions between neurons, are modulated by stimulus identity. Further, we show that a simple machine learning classifier trained using graph-theoretic features alone, or in combination with neural activity features, can accurately predict salt stimulus. Moreover, by focusing on putative causal interactions between neurons, the graph-theoretic features were almost twice as predictive as the neural activity features. These results reveal that stimulus identity modulates the broad, network-level organization of the nervous system, and that graph theory can be used to characterize these changes.

Cancer reflects the dysregulation of the underlying gene network, which is strongly related to the 3D genome organization. Numerous efforts have been spent on experimental characterizations of the structural alterations in cancer genomes. However, there is still a lack of genomic structural-level understanding of the temporal dynamics for cancer initiation and progression. Here, we use a landscape-switching model to investigate the chromosome structural transition during the cancerization and reversion processes. We find that the chromosome undergoes a non-monotonic structural shape-changing pathway with initial expansion followed by compaction during both of these processes. Furthermore, our analysis reveals that the chromosome with a more expanding structure than those at both the normal and cancer cell during cancerization exhibits a sparse contact pattern, which shows significant structural similarity to the one at the embryonic stem cell in many aspects, including the trend of contact probability declining with the genomic distance, the global structural shape geometry and the spatial distribution of loci on the chromosome. In light of the intimate structure-function relationship at the chromosomal level, we further describe the cell state transition processes by the chromosome structural changes, suggesting an elevated cell stemness during the formation of the cancer cells. We show that cell cancerization and reversion are highly irreversible processes in terms of the chromosome structural transition pathways, spatial repositioning of chromosomal loci and hysteresis loop of contact evolution analysis. Our model draws a molecular-scale picture of cell cancerization from the chromosome structural perspective. The process contains initial reprogramming towards the stem cell followed by the differentiation towards the cancer cell, accompanied by an initial increase and subsequent decrease of the cell stemness.

Robustness and plasticity are essential features that allow biological systems to cope with complex and variable environments. In a constant environment, robustness, i.e., insensitivity of phenotypes, is expected to increase, whereas plasticity, i.e., the changeability of phenotypes, tends to diminish. Under a variable environment, existence of plasticity will be relevant. The robustness and plasticity, on the other hand, are related to phenotypic variances. As phenotypic variances decrease with the increase in robustness to perturbations, they are expected to decrease through the evolution. However, in nature, phenotypic fluctuation is preserved to a certain degree. One possible cause for this is environmental variation, where one of the most important "environmental" factors will be inter-species interactions. As a first step toward investigating phenotypic fluctuation in response to an inter-species interaction, we present the study of a simple two-species system that comprises hosts and parasites. Hosts are expected to evolve to achieve a phenotype that optimizes fitness. Then, the robustness of the corresponding phenotype will be increased by reducing phenotypic fluctuations. Conversely, plasticity tends to evolve to avoid certain phenotypes that are attacked by parasites. By using a dynamic model of gene expression for the host, we investigate the evolution of the genotype-phenotype map and of phenotypic variances. If the host-parasite interaction is weak, the fittest phenotype of the host evolves to reduce phenotypic variances. In contrast, if there exists a sufficient degree of interaction, the phenotypic variances of hosts increase to escape parasite attacks. For the latter case, we found two strategies: if the noise in the stochastic gene expression is below a certain threshold, the phenotypic variance increases via genetic diversification, whereas above this threshold, it is increased mediated by noise-induced phenotypic fluctuation. We examine how the increase in the phenotypic variances caused by parasite interactions influences the growth rate of a single host, and observed a trade-off between the two. Our results help elucidate the roles played by noise and genetic mutations in the evolution of phenotypic fluctuation and robustness in response to host-parasite interactions.

With the rapid accumulation of biological omics datasets, decoding the underlying relationships of cross-dataset genes becomes an important issue. Previous studies have attempted to identify differentially expressed genes across datasets. However, it is hard for them to detect interrelated ones. Moreover, existing correlation-based algorithms can only measure the relationship between genes within a single dataset or two multi-modal datasets from the same samples. It is still unclear how to quantify the strength of association of the same gene across two biological datasets with different samples. To this end, we propose Approximate Distance Correlation (ADC) to select interrelated genes with statistical significance across two different biological datasets. ADC first obtains the k most correlated genes for each target gene as its approximate observations, and then calculates the distance correlation (DC) for the target gene across two datasets. ADC repeats this process for all genes and then performs the Benjamini-Hochberg adjustment to control the false discovery rate. We demonstrate the effectiveness of ADC with simulation data and four real applications to select highly interrelated genes across two datasets. These four applications including 21 cancer RNA-seq datasets of different tissues; six single-cell RNA-seq (scRNA-seq) datasets of mouse hematopoietic cells across six different cell types along the hematopoietic cell lineage; five scRNA-seq datasets of pancreatic islet cells across five different technologies; coupled single-cell ATAC-seq (scATAC-seq) and scRNA-seq data of peripheral blood mononuclear cells (PBMC). Extensive results demonstrate that ADC is a powerful tool to uncover interrelated genes with strong biological implications and is scalable to large-scale datasets. Moreover, the number of such genes can serve as a metric to measure the similarity between two datasets, which could characterize the relative difference of diverse cell types and technologies.

In recent literature, the human microbiome has been shown to have a major influence on human health. To investigate this impact, scientists study the composition and abundance of bacterial species, commonly using 16S rRNA gene sequencing, among patients with and without a disease or condition. Methods for such investigations to date have focused on the association between individual bacterium and an outcome, and higher-order pairwise relationships or interactions among bacteria are often avoided due to the substantial increase in dimension and the potential for spurious correlations. However, overlooking such relationships ignores the environment of the microbiome, where there is dynamic cooperation and competition among bacteria. We present a method for identifying and ranking pairs of bacteria that have a differential dichotomized relationship across outcomes. Our approach, implemented in an R package PairSeek, uses the stability selection framework with data-driven dichotomized forms of the pairwise relationships. We illustrate the properties of the proposed method using a published oral cancer data set and a simulation study.

Advances in genetic engineering technologies have allowed the construction of artificial genetic circuits, which have been used to generate spatial patterns of differential gene expression. However, the question of how cells can be programmed, and how complex the rules need to be, to achieve a desired tissue morphology has received less attention. Here, we address these questions by developing a mathematical model to study how cells can collectively grow into clusters with different structural morphologies by secreting diffusible signals that can influence cellular growth rates. We formulate how growth regulators can be used to control the formation of cellular protrusions and how the range of achievable structures scales with the number of distinct signals. We show that a single growth inhibitor is insufficient for the formation of multiple protrusions but may be achieved with multiple growth inhibitors, and that other types of signals can regulate the shape of protrusion tips. These examples illustrate how our approach could potentially be used to guide the design of regulatory circuits for achieving a desired target structure.

In agricultural landscapes, the amount and organization of crops and semi-natural habitats (SNH) have the potential to promote a bundle of ecosystem services due to their influence on ecological community at multiple spatio-temporal scales. SNH are relatively undisturbed and are often source of complementary resources and refuges, therefore supporting more diverse and abundant natural pest enemies. However, the nexus of SNH proportion and organization with pest suppression is not trivial. It is thus crucial to understand how the behavior of pest and natural enemy species, the underlying landscape structure, and their interaction, may influence conservation biological control (CBC). Here, we develop a generative stochastic landscape model to simulate realistic agricultural landscape compositions and configurations of fields and linear elements. Generated landscapes are used as spatial support over which we simulate a spatially explicit predator-prey dynamic model. We find that increased SNH presence boosts predator populations by sustaining high predator density that regulates and keeps pest density below the pesticide application threshold. However, predator presence over all the landscape helps to stabilize the pest population by keeping it under this threshold, which tends to increase pest density at the landscape scale. In addition, the joint effect of SNH presence and predator dispersal ability among hedge and field interface results in a stronger pest regulation, which also limits pest growth. Considering properties of both fields and linear elements, such as local structure and geometric features, provides deeper insights for pest regulation; for example, hedge presence at crop field boundaries clearly strengthens CBC. Our results highlight that the integration of species behaviors and traits with landscape structure at multiple scales is necessary to provide useful insights for CBC.

Nucleotide sequence and taxonomy reference databases are critical resources for widespread applications including marker-gene and metagenome sequencing for microbiome analysis, diet metabarcoding, and environmental DNA (eDNA) surveys. Reproducibly generating, managing, using, and evaluating nucleotide sequence and taxonomy reference databases creates a significant bottleneck for researchers aiming to generate custom sequence databases. Furthermore, database composition drastically influences results, and lack of standardization limits cross-study comparisons. To address these challenges, we developed RESCRIPt, a Python 3 software package and QIIME 2 plugin for reproducible generation and management of reference sequence taxonomy databases, including dedicated functions that streamline creating databases from popular sources, and functions for evaluating, comparing, and interactively exploring qualitative and quantitative characteristics across reference databases. To highlight the breadth and capabilities of RESCRIPt, we provide several examples for working with popular databases for microbiome profiling (SILVA, Greengenes, NCBI-RefSeq, GTDB), eDNA and diet metabarcoding surveys (BOLD, GenBank), as well as for genome comparison. We show that bigger is not always better, and reference databases with standardized taxonomies and those that focus on type strains have quantitative advantages, though may not be appropriate for all use cases. Most databases appear to benefit from some curation (quality filtering), though sequence clustering appears detrimental to database quality. Finally, we demonstrate the breadth and extensibility of RESCRIPt for reproducible workflows with a comparison of global hepatitis genomes. RESCRIPt provides tools to democratize the process of reference database acquisition and management, enabling researchers to reproducibly and transparently create reference materials for diverse research applications. RESCRIPt is released under a permissive BSD-3 license at https://github.com/bokulich-lab/RESCRIPt.

Genome-scale metabolic models (GEMs) are comprehensive knowledge bases of cellular metabolism and serve as mathematical tools for studying biological phenotypes and metabolic states or conditions in various organisms and cell types. Given the sheer size and complexity of human metabolism, selecting parameters for existing analysis methods such as metabolic objective functions and model constraints is not straightforward in human GEMs. In particular, comparing several conditions in large GEMs to identify condition- or disease-specific metabolic features is challenging. In this study, we showcase a scalable, model-driven approach for an in-depth investigation and comparison of metabolic states in large GEMs which enables identifying the underlying functional differences. Using a combination of flux space sampling and network analysis, our approach enables extraction and visualisation of metabolically distinct network modules. Importantly, it does not rely on known or assumed objective functions. We apply this novel approach to extract the biochemical differences in adipocytes arising due to unlimited vs blocked uptake of branched-chain amino acids (BCAAs, considered as biomarkers in obesity) using a human adipocyte GEM (iAdipocytes1809). The biological significance of our approach is corroborated by literature reports confirming our identified metabolic processes (TCA cycle and Fatty acid metabolism) to be functionally related to BCAA metabolism. Additionally, our analysis predicts a specific altered uptake and secretion profile indicating a compensation for the unavailability of BCAAs. Taken together, our approach facilitates determining functional differences between any metabolic conditions of interest by offering a versatile platform for analysing and comparing flux spaces of large metabolic networks.

New microbial communities often arise through the mixing of two or more separately assembled parent communities, a phenomenon that has been termed "community coalescence". Understanding how the interaction structures of complex parent communities determine the outcomes of coalescence events is an important challenge. While recent work has begun to elucidate the role of competition in coalescence, that of cooperation, a key interaction type commonly seen in microbial communities, is still largely unknown. Here, using a general consumer-resource model, we study the combined effects of competitive and cooperative interactions on the outcomes of coalescence events. To do so, we simulate coalescence events between pairs of communities with different degrees of competition for shared carbon resources and cooperation through cross-feeding on leaked metabolic by-products (facilitation). We also study how structural and functional properties of post-coalescence communities evolve when they are subjected to repeated coalescence events. We find that in coalescence events, the less competitive and more cooperative parent communities contribute a higher proportion of species to the new community because of their superior ability to deplete resources and resist invasions. Consequently, when a community is subjected to repeated coalescence events, it gradually evolves towards being less competitive and more cooperative, as well as more speciose, robust and efficient in resource use. Encounters between microbial communities are becoming increasingly frequent as a result of anthropogenic environmental change, and there is great interest in how the coalescence of microbial communities affects environmental and human health. Our study provides new insights into the mechanisms behind microbial community coalescence, and a framework to predict outcomes based on the interaction structures of parent communities.

Metabolic network models are increasingly being used in health care and industry. As a consequence, many tools have been released to automate their reconstruction process de novo. In order to enable gene deletion simulations and integration of gene expression data, these networks must include gene-protein-reaction (GPR) rules, which describe with a Boolean logic relationships between the gene products (e.g., enzyme isoforms or subunits) associated with the catalysis of a given reaction. Nevertheless, the reconstruction of GPRs still remains a largely manual and time consuming process. Aiming at fully automating the reconstruction process of GPRs for any organism, we propose the open-source python-based framework GPRuler. By mining text and data from 9 different biological databases, GPRuler can reconstruct GPRs starting either from just the name of the target organism or from an existing metabolic model. The performance of the developed tool is evaluated at small-scale level for a manually curated metabolic model, and at genome-scale level for three metabolic models related to Homo sapiens and Saccharomyces cerevisiae organisms. By exploiting these models as benchmarks, the proposed tool shown its ability to reproduce the original GPR rules with a high level of accuracy. In all the tested scenarios, after a manual investigation of the mismatches between the rules proposed by GPRuler and the original ones, the proposed approach revealed to be in many cases more accurate than the original models. By complementing existing tools for metabolic network reconstruction with the possibility to reconstruct GPRs quickly and with a few resources, GPRuler paves the way to the study of context-specific metabolic networks, representing the active portion of the complete network in given conditions, for organisms of industrial or biomedical interest that have not been characterized metabolically yet.

Cortical pyramidal cells (PCs) have a specialized dendritic mechanism for the generation of bursts, suggesting that these events play a special role in cortical information processing. In vivo, bursts occur at a low, but consistent rate. Theory suggests that this network state increases the amount of information they convey. However, because burst activity relies on a threshold mechanism, it is rather sensitive to dendritic input levels. In spiking network models, network states in which bursts occur rarely are therefore typically not robust, but require fine-tuning. Here, we show that this issue can be solved by a homeostatic inhibitory plasticity rule in dendrite-targeting interneurons that is consistent with experimental data. The suggested learning rule can be combined with other forms of inhibitory plasticity to self-organize a network state in which both spikes and bursts occur asynchronously and irregularly at low rate. Finally, we show that this network state creates the network conditions for a recently suggested multiplexed code and thereby indeed increases the amount of information encoded in bursts.

The use of unnatural fluorogenic molecules widely expands the pallet of available genetically encoded fluorescent imaging tools through the design of fluorogen activating proteins (FAPs). While there is already a handful of such probes available, each of them went through laborious cycles of in vitro screening and selection. Computational modeling approaches are evolving incredibly fast right now and are demonstrating great results in many applications, including de novo protein design. It suggests that the easier task of fine-tuning the fluorogen-binding properties of an already functional protein in silico should be readily achievable. To test this hypothesis, we used Rosetta for computational ligand docking followed by protein binding pocket redesign to further improve the previously described FAP DiB1 that is capable of binding to a BODIPY-like dye M739. Despite an inaccurate initial docking of the chromophore, the incorporated mutations nevertheless improved multiple photophysical parameters as well as the overall performance of the tag. The designed protein, DiB-RM, shows higher brightness, localization precision, and apparent photostability in protein-PAINT super-resolution imaging compared to its parental variant DiB1. Moreover, DiB-RM can be cleaved to obtain an efficient split system with enhanced performance compared to a parental DiB-split system. The possible reasons for the inaccurate ligand binding pose prediction and its consequence on the outcome of the design experiment are further discussed.

Working memory (WM) deficits have been widely documented in schizophrenia (SZ), and almost all existing studies attributed the deficits to decreased capacity as compared to healthy control (HC) subjects. Recent developments in WM research suggest that other components, such as precision, also mediate behavioral performance. It remains unclear how different WM components jointly contribute to deficits in schizophrenia. We measured the performance of 60 SZ (31 females) and 61 HC (29 females) in a classical delay-estimation visual working memory (VWM) task and evaluated several influential computational models proposed in basic science of VWM to disentangle the effect of various memory components. We show that the model assuming variable precision (VP) across items and trials is the best model to explain the performance of both groups. According to the VP model, SZ exhibited abnormally larger variability of allocating memory resources rather than resources or capacity per se. Finally, individual differences in the resource allocation variability predicted variation of symptom severity in SZ, highlighting its functional relevance to schizophrenic pathology. This finding was further verified using distinct visual features and subject cohorts. These results provide an alternative view instead of the widely accepted decreased-capacity theory and highlight the key role of elevated resource allocation variability in generating atypical VWM behavior in schizophrenia. Our findings also shed new light on the utility of Bayesian observer models to characterize mechanisms of mental deficits in clinical neuroscience.

We present a comprehensive, experimental and theoretical study of the impact of 5-hydroxymethylation of DNA cytosine. Using molecular dynamics, biophysical experiments and NMR spectroscopy, we found that Ten-Eleven translocation (TET) dioxygenases generate an epigenetic variant with structural and physical properties similar to those of 5-methylcytosine. Experiments and simulations demonstrate that 5-methylcytosine (mC) and 5-hydroxymethylcytosine (hmC) generally lead to stiffer DNA than normal cytosine, with poorer circularization efficiencies and lower ability to form nucleosomes. In particular, we can rule out the hypothesis that hydroxymethylation reverts to unmodified cytosine physical properties, as hmC is even more rigid than mC. Thus, we do not expect dramatic changes in the chromatin structure induced by differences in physical properties between d(mCpG) and d(hmCpG). Conversely, our simulations suggest that methylated-DNA binding domains (MBDs), associated with repression activities, are sensitive to the substitution d(mCpG) ➔ d(hmCpG), while MBD3 which has a dual activation/repression activity is not sensitive to the d(mCpG) d(hmCpG) change. Overall, while gene activity changes due to cytosine methylation are the result of the combination of stiffness-related chromatin reorganization and MBD binding, those associated to 5-hydroxylation of methylcytosine could be explained by a change in the balance of repression/activation pathways related to differential MBD binding.

To help cells cope with protein misfolding and aggregation, Hsp70 molecular chaperones selectively bind a variety of sequences ("selective promiscuity"). Statistical analyses from substrate-derived peptide arrays reveal that DnaK, the E. coli Hsp70, binds to sequences containing three to five branched hydrophobic residues, although otherwise the specific amino acids can vary considerably. Several high-resolution structures of the substrate -binding domain (SBD) of DnaK bound to peptides reveal a highly conserved configuration of the bound substrate and further suggest that the substrate-binding cleft consists of five largely independent sites for interaction with five consecutive substrate residues. Importantly, both substrate backbone orientations (N- to C- and C- to N-) allow essentially the same backbone hydrogen-bonding and side-chain interactions with the chaperone. In order to rationalize these observations, we performed atomistic molecular dynamics simulations to sample the interactions of all 20 amino acid side chains in each of the five sites of the chaperone in the context of the conserved substrate backbone configurations. The resulting interaction energetics provide the basis set for deriving a predictive model that we call Paladin (Physics-based model of DnaK-Substrate Binding). Trained using available peptide array data, Paladin can distinguish binders and nonbinders of DnaK with accuracy comparable to existing predictors and further predicts the detailed configuration of the bound sequence. Tested using existing DnaK-peptide structures, Paladin correctly predicted the binding register in 10 out of 13 substrate sequences that bind in the N- to C- orientation, and the binding orientation in 16 out of 22 sequences. The physical basis of the Paladin model provides insight into the origins of how Hsp70s bind substrates with a balance of selectivity and promiscuity. The approach described here can be extended to other Hsp70s where extensive peptide array data is not available.

Very high risk neuroblastoma is characterised by increased MAPK signalling, and targeting MAPK signalling is a promising therapeutic strategy. We used a deeply characterised panel of neuroblastoma cell lines and found that the sensitivity to MEK inhibitors varied drastically between these cell lines. By generating quantitative perturbation data and mathematical modelling, we determined potential resistance mechanisms. We found that negative feedbacks within MAPK signalling and via the IGF receptor mediate re-activation of MAPK signalling upon treatment in resistant cell lines. By using cell-line specific models, we predict that combinations of MEK inhibitors with RAF or IGFR inhibitors can overcome resistance, and tested these predictions experimentally. In addition, phospho-proteomic profiling confirmed the cell-specific feedback effects and synergy of MEK and IGFR targeted treatment. Our study shows that a quantitative understanding of signalling and feedback mechanisms facilitated by models can help to develop and optimise therapeutic strategies. Our findings should be considered for the planning of future clinical trials introducing MEKi in the treatment of neuroblastoma.

The ability to treat gonorrhoea with current first-line drugs is threatened by the global spread of extensively drug resistant (XDR) Neisseria gonorrhoeae (NG) strains. In Australia, urban transmission is high among men who have sex with men (MSM) and importation of an XDR NG strain in this population could result in an epidemic that would be difficult and costly to control. An individual-based, anatomical site-specific mathematical model of NG transmission among Australian MSM was developed and used to evaluate the potential for elimination of an imported NG strain under a range of case-based and population-based test-and-treat strategies. When initiated upon detection of the imported strain, these strategies enhance the probability of elimination and reduce the outbreak size compared with current practice (current testing levels and no contact tracing). The most effective strategies combine testing targeted at regular and casual partners with increased rates of population testing. However, even with the most effective strategies, outbreaks can persist for up to 2 years post-detection. Our simulations suggest that local elimination of imported NG strains can be achieved with high probability using combined case-based and population-based test-and-treat strategies. These strategies may be an effective means of preserving current treatments in the event of wider XDR NG emergence.

It is generally accepted that the number of neurons in a given brain area far exceeds the number of neurons needed to carry any specific function controlled by that area. For example, motor areas of the human brain contain tens of millions of neurons that control the activation of tens or at most hundreds of muscles. This massive redundancy implies the covariation of many neurons, which constrains the population activity to a low-dimensional manifold within the space of all possible patterns of neural activity. To gain a conceptual understanding of the complexity of the neural activity within a manifold, it is useful to estimate its dimensionality, which quantifies the number of degrees of freedom required to describe the observed population activity without significant information loss. While there are many algorithms for dimensionality estimation, we do not know which are well suited for analyzing neural activity. The objective of this study was to evaluate the efficacy of several representative algorithms for estimating the dimensionality of linearly and nonlinearly embedded data. We generated synthetic neural recordings with known intrinsic dimensionality and used them to test the algorithms' accuracy and robustness. We emulated some of the important challenges associated with experimental data by adding noise, altering the nature of the embedding of the low-dimensional manifold within the high-dimensional recordings, varying the dimensionality of the manifold, and limiting the amount of available data. We demonstrated that linear algorithms overestimate the dimensionality of nonlinear, noise-free data. In cases of high noise, most algorithms overestimated the dimensionality. We thus developed a denoising algorithm based on deep learning, the "Joint Autoencoder", which significantly improved subsequent dimensionality estimation. Critically, we found that all algorithms failed when the intrinsic dimensionality was high (above 20) or when the amount of data used for estimation was low. Based on the challenges we observed, we formulated a pipeline for estimating the dimensionality of experimental neural data.

Finding out the physical structure of neuronal circuits that governs neuronal responses is an important goal for brain research. With fast advances for large-scale recording techniques, identification of a neuronal circuit with multiple neurons and stages or layers becomes possible and highly demanding. Although methods for mapping the connection structure of circuits have been greatly developed in recent years, they are mostly limited to simple scenarios of a few neurons in a pairwise fashion; and dissecting dynamical circuits, particularly mapping out a complete functional circuit that converges to a single neuron, is still a challenging question. Here, we show that a recent method, termed spike-triggered non-negative matrix factorization (STNMF), can address these issues. By simulating different scenarios of spiking neural networks with various connections between neurons and stages, we demonstrate that STNMF is a persuasive method to dissect functional connections within a circuit. Using spiking activities recorded at neurons of the output layer, STNMF can obtain a complete circuit consisting of all cascade computational components of presynaptic neurons, as well as their spiking activities. For simulated simple and complex cells of the primary visual cortex, STNMF allows us to dissect the pathway of visual computation. Taken together, these results suggest that STNMF could provide a useful approach for investigating neuronal systems leveraging recorded functional neuronal activity.

Predictive approaches such as virtual screening have been used in drug discovery with the objective of reducing developmental time and costs. Current machine learning and network-based approaches have issues related to generalization, usability, or model interpretability, especially due to the complexity of target proteins' structure/function, and bias in system training datasets. Here, we propose a new method "DRUIDom" (DRUg Interacting Domain prediction) to identify bio-interactions between drug candidate compounds and targets by utilizing the domain modularity of proteins, to overcome problems associated with current approaches. DRUIDom is composed of two methodological steps. First, ligands/compounds are statistically mapped to structural domains of their target proteins, with the aim of identifying their interactions. As such, other proteins containing the same mapped domain or domain pair become new candidate targets for the corresponding compounds. Next, a million-scale dataset of small molecule compounds, including those mapped to domains in the previous step, are clustered based on their molecular similarities, and their domain associations are propagated to other compounds within the same clusters. Experimentally verified bioactivity data points, obtained from public databases, are meticulously filtered to construct datasets of active/interacting and inactive/non-interacting drug/compound-target pairs (~2.9M data points), and used as training data for calculating parameters of compound-domain mappings, which led to 27,032 high-confidence associations between 250 domains and 8,165 compounds, and a finalized output of ~5 million new compound-protein interactions. DRUIDom is experimentally validated by syntheses and bioactivity analyses of compounds predicted to target LIM-kinase proteins, which play critical roles in the regulation of cell motility, cell cycle progression, and differentiation through actin filament dynamics. We showed that LIMK-inhibitor-2 and its derivatives significantly block the cancer cell migration through inhibition of LIMK phosphorylation and the downstream protein cofilin. One of the derivative compounds (LIMKi-2d) was identified as a promising candidate due to its action on resistant Mahlavu liver cancer cells. The results demonstrated that DRUIDom can be exploited to identify drug candidate compounds for intended targets and to predict new target proteins based on the defined compound-domain relationships. Datasets, results, and the source code of DRUIDom are fully-available at: https://github.com/cansyl/DRUIDom.

Sickle cell disease, a genetic disorder affecting a sizeable global demographic, manifests in sickle red blood cells (sRBCs) with altered shape and biomechanics. sRBCs show heightened adhesive interactions with inflamed endothelium, triggering painful vascular occlusion events. Numerous studies employ microfluidic-assay-based monitoring tools to quantify characteristics of adhered sRBCs from high resolution channel images. The current image analysis workflow relies on detailed morphological characterization and cell counting by a specially trained worker. This is time and labor intensive, and prone to user bias artifacts. Here we establish a morphology based classification scheme to identify two naturally arising sRBC subpopulations-deformable and non-deformable sRBCs-utilizing novel visual markers that link to underlying cell biomechanical properties and hold promise for clinically relevant insights. We then set up a standardized, reproducible, and fully automated image analysis workflow designed to carry out this classification. This relies on a two part deep neural network architecture that works in tandem for segmentation of channel images and classification of adhered cells into subtypes. Network training utilized an extensive data set of images generated by the SCD BioChip, a microfluidic assay which injects clinical whole blood samples into protein-functionalized microchannels, mimicking physiological conditions in the microvasculature. Here we carried out the assay with the sub-endothelial protein laminin. The machine learning approach segmented the resulting channel images with 99.1±0.3% mean IoU on the validation set across 5 k-folds, classified detected sRBCs with 96.0±0.3% mean accuracy on the validation set across 5 k-folds, and matched trained personnel in overall characterization of whole channel images with R2 = 0.992, 0.987 and 0.834 for total, deformable and non-deformable sRBC counts respectively. Average analysis time per channel image was also improved by two orders of magnitude (∼ 2 minutes vs ∼ 2-3 hours) over manual characterization. Finally, the network results show an order of magnitude less variance in counts on repeat trials than humans. This kind of standardization is a prerequisite for the viability of any diagnostic technology, making our system suitable for affordable and high throughput disease monitoring.

Argonaute proteins in combination with short microRNA (miRNAs) can target mRNA molecules for translation inhibition or degradation and play a key role in many regulatory processes. The miRNAs act as guide RNAs that associate with Argonaute and the complementary mRNA target region. The complex formation results in activation of Argonaute and specific cleavage of the target mRNA. Both the binding and activation processes involve essential domain rearrangements of functional importance. For the Thermus Thermophilus Argonaute (TtAgo) system guide-bound (binary) and guide/target-bound (ternary) complexes are known but how the binding of guide and target mediate domain movements is still not understood. We have studied the Argonaute domain motion in apo and guide/target bound states using Molecular Dynamics simulations and a Hamiltonian replica exchange (H-REMD) method that employs a specific biasing potential to accelerate domain motions. The H-REMD technique indicates sampling of a much broader distribution of domain arrangements both in the apo as well as binary and ternary complexes compared to regular MD simulations. In the apo state domain arrangements corresponding to more compact (closed) states are mainly sampled which undergo an opening upon guide and guide/target binding. Whereas only limited overlap in domain geometry between apo and bound states was found, a larger similarity in the domain distribution is observed for the simulations of binary and ternary complexes. Comparative simulations on ternary complexes with 15 or 16 base pairs (bp) formed between guide and target strands (instead of 14) resulted in dissociation of the 3'-guide strand from the PAZ domain and domain rearrangement. This agrees with the experimental observation that guide-target pairing beyond 14 bps is required for activation and gives a mechanistic explanation for the experimentally observed activation process.

Hybrid multiscale agent-based models (ABMs) are unique in their ability to simulate individual cell interactions and microenvironmental dynamics. Unfortunately, the high computational cost of modeling individual cells, the inherent stochasticity of cell dynamics, and numerous model parameters are fundamental limitations of applying such models to predict tumor dynamics. To overcome these challenges, we have developed a coarse-grained two-scale ABM (cgABM) with a reduced parameter space that allows for an accurate and efficient calibration using a set of time-resolved microscopy measurements of cancer cells grown with different initial conditions. The multiscale model consists of a reaction-diffusion type model capturing the spatio-temporal evolution of glucose and growth factors in the tumor microenvironment (at tissue scale), coupled with a lattice-free ABM to simulate individual cell dynamics (at cellular scale). The experimental data consists of BT474 human breast carcinoma cells initialized with different glucose concentrations and tumor cell confluences. The confluence of live and dead cells was measured every three hours over four days. Given this model, we perform a time-dependent global sensitivity analysis to identify the relative importance of the model parameters. The subsequent cgABM is calibrated within a Bayesian framework to the experimental data to estimate model parameters, which are then used to predict the temporal evolution of the living and dead cell populations. To this end, a moment-based Bayesian inference is proposed to account for the stochasticity of the cgABM while quantifying uncertainties due to limited temporal observational data. The cgABM reduces the computational time of ABM simulations by 93% to 97% while staying within a 3% difference in prediction compared to ABM. Additionally, the cgABM can reliably predict the temporal evolution of breast cancer cells observed by the microscopy data with an average error and standard deviation for live and dead cells being 7.61±2.01 and 5.78±1.13, respectively.

Visual stimuli are represented by a highly efficient code in the primary visual cortex, but the development of this code is still unclear. Two distinct factors control coding efficiency: Representational efficiency, which is determined by neuronal tuning diversity, and metabolic efficiency, which is influenced by neuronal gain. How these determinants of coding efficiency are shaped during development, supported by excitatory and inhibitory plasticity, is only partially understood. We investigate a fully plastic spiking network of the primary visual cortex, building on phenomenological plasticity rules. Our results suggest that inhibitory plasticity is key to the emergence of tuning diversity and accurate input encoding. We show that inhibitory feedback (random and specific) increases the metabolic efficiency by implementing a gain control mechanism. Interestingly, this led to the spontaneous emergence of contrast-invariant tuning curves. Our findings highlight that (1) interneuron plasticity is key to the development of tuning diversity and (2) that efficient sensory representations are an emergent property of the resulting network.

Making good decisions requires updating beliefs according to new evidence. This is a dynamical process that is prone to biases: in some cases, beliefs become entrenched and resistant to new evidence (leading to primacy effects), while in other cases, beliefs fade over time and rely primarily on later evidence (leading to recency effects). How and why either type of bias dominates in a given context is an important open question. Here, we study this question in classic perceptual decision-making tasks, where, puzzlingly, previous empirical studies differ in the kinds of biases they observe, ranging from primacy to recency, despite seemingly equivalent tasks. We present a new model, based on hierarchical approximate inference and derived from normative principles, that not only explains both primacy and recency effects in existing studies, but also predicts how the type of bias should depend on the statistics of stimuli in a given task. We verify this prediction in a novel visual discrimination task with human observers, finding that each observer's temporal bias changed as the result of changing the key stimulus statistics identified by our model. The key dynamic that leads to a primacy bias in our model is an overweighting of new sensory information that agrees with the observer's existing belief-a type of 'confirmation bias'. By fitting an extended drift-diffusion model to our data we rule out an alternative explanation for primacy effects due to bounded integration. Taken together, our results resolve a major discrepancy among existing perceptual decision-making studies, and suggest that a key source of bias in human decision-making is approximate hierarchical inference.

Cognitive biases are widespread in humans and animals alike, and can sometimes be reinforced by social interactions. One prime bias in judgment and decision-making is the human tendency to underestimate large quantities. Previous research on social influence in estimation tasks has generally focused on the impact of single estimates on individual and collective accuracy, showing that randomly sharing estimates does not reduce the underestimation bias. Here, we test a method of social information sharing that exploits the known relationship between the true value and the level of underestimation, and study if it can counteract the underestimation bias. We performed estimation experiments in which participants had to estimate a series of quantities twice, before and after receiving estimates from one or several group members. Our purpose was threefold: to study (i) whether restructuring the sharing of social information can reduce the underestimation bias, (ii) how the number of estimates received affects the sensitivity to social influence and estimation accuracy, and (iii) the mechanisms underlying the integration of multiple estimates. Our restructuring of social interactions successfully countered the underestimation bias. Moreover, we find that sharing more than one estimate also reduces the underestimation bias. Underlying our results are a human tendency to herd, to trust larger estimates than one's own more than smaller estimates, and to follow disparate social information less. Using a computational modeling approach, we demonstrate that these effects are indeed key to explain the experimental results. Overall, our results show that existing knowledge on biases can be used to dampen their negative effects and boost judgment accuracy, paving the way for combating other cognitive biases threatening collective systems.

Signaling networks mediate many aspects of cellular function. The conventional, mechanistically motivated approach to modeling such networks is through mass-action chemistry, which maps directly to biological entities and facilitates experimental tests and predictions. However such models are complex, need many parameters, and are computationally costly. Here we introduce the HillTau form for signaling models. HillTau retains the direct mapping to biological observables, but it uses far fewer parameters, and is 100 to over 1000 times faster than ODE-based methods. In the HillTau formalism, the steady-state concentration of signaling molecules is approximated by the Hill equation, and the dynamics by a time-course tau. We demonstrate its use in implementing several biochemical motifs, including association, inhibition, feedforward and feedback inhibition, bistability, oscillations, and a synaptic switch obeying the BCM rule. The major use-cases for HillTau are system abstraction, model reduction, scaffolds for data-driven optimization, and fast approximations to complex cellular signaling.

Patients with coronavirus disease 2019 (COVID-19) often exhibit diverse disease progressions associated with various infectious ability, symptoms, and clinical treatments. To systematically and thoroughly understand the heterogeneous progression of COVID-19, we developed a multi-scale computational model to quantitatively understand the heterogeneous progression of COVID-19 patients infected with severe acute respiratory syndrome (SARS)-like coronavirus (SARS-CoV-2). The model consists of intracellular viral dynamics, multicellular infection process, and immune responses, and was formulated using a combination of differential equations and stochastic modeling. By integrating multi-source clinical data with model analysis, we quantified individual heterogeneity using two indexes, i.e., the ratio of infected cells and incubation period. Specifically, our simulations revealed that increasing the host antiviral state or virus induced type I interferon (IFN) production rate can prolong the incubation period and postpone the transition from asymptomatic to symptomatic outcomes. We further identified the threshold dynamics of T cell exhaustion in the transition between mild-moderate and severe symptoms, and that patients with severe symptoms exhibited a lack of naïve T cells at a late stage. In addition, we quantified the efficacy of treating COVID-19 patients and investigated the effects of various therapeutic strategies. Simulations results suggested that single antiviral therapy is sufficient for moderate patients, while combination therapies and prevention of T cell exhaustion are needed for severe patients. These results highlight the critical roles of IFN and T cell responses in regulating the stage transition during COVID-19 progression. Our study reveals a quantitative relationship underpinning the heterogeneity of transition stage during COVID-19 progression and can provide a potential guidance for personalized therapy in COVID-19 patients.

Transmembranal G Protein-Coupled Receptors (GPCRs) transduce extracellular chemical signals to the cell, via conformational change from a resting (inactive) to an active (canonically bound to a G-protein) conformation. Receptor activation is normally modulated by extracellular ligand binding, but mutations in the receptor can also shift this equilibrium by stabilizing different conformational states. In this work, we built structure-energetic relationships of receptor activation based on original thermodynamic cycles that represent the conformational equilibrium of the prototypical A2A adenosine receptor (AR). These cycles were solved with efficient free energy perturbation (FEP) protocols, allowing to distinguish the pharmacological profile of different series of A2AAR agonists with different efficacies. The modulatory effects of point mutations on the basal activity of the receptor or on ligand efficacies could also be detected. This methodology can guide GPCR ligand design with tailored pharmacological properties, or allow the identification of mutations that modulate receptor activation with potential clinical implications.

The energetic economy of running benefits from tendon and other tissues that store and return elastic energy, thus saving muscles from costly mechanical work. The classic "Spring-mass" computational model successfully explains the forces, displacements and mechanical power of running, as the outcome of dynamical interactions between the body center of mass and a purely elastic spring for the leg. However, the Spring-mass model does not include active muscles and cannot explain the metabolic energy cost of running, whether on level ground or on a slope. Here we add explicit actuation and dissipation to the Spring-mass model, and show how they explain substantial active (and thus costly) work during human running, and much of the associated energetic cost. Dissipation is modeled as modest energy losses (5% of total mechanical energy for running at 3 m s-1) from hysteresis and foot-ground collisions, that must be restored by active work each step. Even with substantial elastic energy return (59% of positive work, comparable to empirical observations), the active work could account for most of the metabolic cost of human running (about 68%, assuming human-like muscle efficiency). We also introduce a previously unappreciated energetic cost for rapid production of force, that helps explain the relatively smooth ground reaction forces of running, and why muscles might also actively perform negative work. With both work and rapid force costs, the model reproduces the energetics of human running at a range of speeds on level ground and on slopes. Although elastic return is key to energy savings, there are still losses that require restorative muscle work, which can cost substantial energy during running.

With the exponential growth of sequence information stored over the last decade, including that of de novo assembled contigs from RNA-Seq experiments, quantification of chimeric sequences has become essential when assembling read data. In transcriptomics, de novo assembled chimeras can closely resemble underlying transcripts, but patterns such as those seen between co-evolving sites, or mapped read counts, become obscured. We have created a de Bruijn based de novo assembler for RNA-Seq data that utilizes a classification system to describe the complexity of underlying graphs from which contigs are created. Each contig is labelled with one of three levels, indicating whether or not ambiguous paths exist. A by-product of this is information on the range of complexity of the underlying gene families present. As a demonstration of CStones ability to assemble high-quality contigs, and to label them in this manner, both simulated and real data were used. For simulated data, ten million read pairs were generated from cDNA libraries representing four species, Drosophila melanogaster, Panthera pardus, Rattus norvegicus and Serinus canaria. These were assembled using CStone, Trinity and rnaSPAdes; the latter two being high-quality, well established, de novo assembers. For real data, two RNA-Seq datasets, each consisting of ≈30 million read pairs, representing two adult D. melanogaster whole-body samples were used. The contigs that CStone produced were comparable in quality to those of Trinity and rnaSPAdes in terms of length, sequence identity of aligned regions and the range of cDNA transcripts represented, whilst providing additional information on chimerism. Here we describe the details of CStones assembly and classification process, and propose that similar classification systems can be incorporated into other de novo assembly tools. Within a related side study, we explore the effects that chimera's within reference sets have on the identification of differentially expression genes. CStone is available at: https://sourceforge.net/projects/cstone/.

Mosquitoes vector harmful pathogens that infect millions of people every year, and developing approaches to effectively control mosquitoes is a topic of great interest. However, the success of many control measures is highly dependent upon ecological, physiological, and life history traits of mosquito species. The behavior of mosquitoes and their potential to vector pathogens can also be impacted by these traits. One trait of interest is mosquito body mass, which depends upon many factors associated with the environment in which juvenile mosquitoes develop. Our experiments examined the impact of larval density on the body mass of Aedes aegypti mosquitoes, which are important vectors of dengue, Zika, yellow fever, and other pathogens. To investigate the interactions between the larval environment and mosquito body mass, we built a discrete time mathematical model that incorporates body mass, larval density, and food availability and fit the model to our experimental data. We considered three categories of model complexity informed by data, and selected the best model within each category using Akaike's Information Criterion. We found that the larval environment is an important determinant of the body mass of mosquitoes upon emergence. Furthermore, we found that larval density has greater impact on body mass of adults at emergence than on development time, and that inclusion of density dependence in the survival of female aquatic stages in models is important. We discuss the implications of our results for the control of Aedes mosquitoes and on their potential to spread disease.

A major challenge in designing proteins de novo to bind user-defined ligands with high affinity is finding backbones structures into which a new binding site geometry can be engineered with high precision. Recent advances in methods to generate protein fold families de novo have expanded the space of accessible protein structures, but it is not clear to what extend de novo proteins with diverse geometries also expand the space of designable ligand binding functions. We constructed a library of 25,806 high-quality ligand binding sites and developed a fast protocol to place ("match") these binding sites into both naturally occurring and de novo protein families with two fold topologies: Rossman and NTF2. Each matching step involves engineering new binding site residues into each protein "scaffold", which is distinct from the problem of comparing already existing binding pockets. 5,896 and 7,475 binding sites could be matched to the Rossmann and NTF2 fold families, respectively. De novo designed Rossman and NTF2 protein families can support 1,791 and 678 binding sites that cannot be matched to naturally existing structures with the same topologies, respectively. While the number of protein residues in ligand binding sites is the major determinant of matching success, ligand size and primary sequence separation of binding site residues also play important roles. The number of matched binding sites are power law functions of the number of members in a fold family. Our results suggest that de novo sampling of geometric variations on diverse fold topologies can significantly expand the space of designable ligand binding sites for a wealth of possible new protein functions.

Pronounced activity is observed in both hemispheres of the motor cortex during preparation and execution of unimanual movements. The organizational principles of bi-hemispheric signals and the functions they serve throughout motor planning remain unclear. Using an instructed-delay reaching task in monkeys, we identified two components in population responses spanning PMd and M1. A "dedicated" component, which segregated activity at the level of individual units, emerged in PMd during preparation. It was most prominent following movement when M1 became strongly engaged, and principally involved the contralateral hemisphere. In contrast to recent reports, these dedicated signals solely accounted for divergence of arm-specific neural subspaces. The other "distributed" component mixed signals for each arm within units, and the subspace containing it did not discriminate between arms at any stage. The statistics of the population response suggest two functional aspects of the cortical network: one that spans both hemispheres for supporting preparatory and ongoing processes, and another that is predominantly housed in the contralateral hemisphere and specifies unilateral output.

The output of neocortical layer 5 pyramidal cells (L5PCs) is expressed by a train of single spikes with intermittent bursts of multiple spikes at high frequencies. The bursts are the result of nonlinear dendritic properties, including Na+, Ca2+, and NMDA spikes, that interact with the ~10,000 synapses impinging on the neuron's dendrites. Output spike bursts are thought to implement key dendritic computations, such as coincidence detection of bottom-up inputs (arriving mostly at the basal tree) and top-down inputs (arriving mostly at the apical tree). In this study we used a detailed nonlinear model of L5PC receiving excitatory and inhibitory synaptic inputs to explore the conditions for generating bursts and for modulating their properties. We established the excitatory input conditions on the basal versus the apical tree that favor burst and show that there are two distinct types of bursts. Bursts consisting of 3 or more spikes firing at < 200 Hz, which are generated by stronger excitatory input to the basal versus the apical tree, and bursts of ~2-spikes at ~250 Hz, generated by prominent apical tuft excitation. Localized and well-timed dendritic inhibition on the apical tree differentially modulates Na+, Ca2+, and NMDA spikes and, consequently, finely controls the burst output. Finally, we explored the implications of different burst classes and respective dendritic inhibition for regulating synaptic plasticity.

The collective behavior of the nuclear array in Drosophila embryos during nuclear cycle (NC) 11 to NC14 is crucial in controlling cell size, establishing developmental patterns, and coordinating morphogenesis. After live imaging on Drosophila embryos with light sheet microscopy, we extract the nuclear trajectory, speed, and internuclear distance with an automatic nuclear tracing method. We find that the nuclear speed shows a period of standing waves along the anterior-posterior (AP) axis after each metaphase as the nuclei collectively migrate towards the embryo poles and partially move back. And the maximum nuclear speed dampens by 28-45% in the second half of the standing wave. Moreover, the nuclear density is 22-42% lower in the pole region than the middle of the embryo during the interphase of NC12-14. To find mechanical rules controlling the collective motion and packing patterns of the nuclear array, we use a deep neural network (DNN) to learn the underlying force field from data. We apply the learned spatiotemporal attractive force field in the simulations with a particle-based model. And the simulations recapitulate nearly all the observed characteristic collective behaviors of nuclear arrays in Drosophila embryos.

Organisms are composed of various cell types with specific states. To obtain a comprehensive understanding of the functions of organs and tissues, cell types have been classified and defined by identifying specific marker genes. Statistical tests are critical for identifying marker genes, which often involve evaluating differences in the mean expression levels of genes. Differentially expressed gene (DEG)-based analysis has been the most frequently used method of this kind. However, in association with increases in sample size such as in single-cell analysis, DEG-based analysis has faced difficulties associated with the inflation of P-values. Here, we propose the concept of discriminative feature of cells (DFC), an alternative to using DEG-based approaches. We implemented DFC using logistic regression with an adaptive LASSO penalty to perform binary classification for discriminating a population of interest and variable selection to obtain a small subset of defining genes. We demonstrated that DFC prioritized gene pairs with non-independent expression using artificial data and that DFC enabled characterization of the muscle satellite/progenitor cell population. The results revealed that DFC well captured cell-type-specific markers, specific gene expression patterns, and subcategories of this cell population. DFC may complement DEG-based methods for interpreting large data sets. DEG-based analysis uses lists of genes with differences in expression between groups, while DFC, which can be termed a discriminative approach, has potential applications in the task of cell characterization. Upon recent advances in the high-throughput analysis of single cells, methods of cell characterization such as scRNA-seq can be effectively subjected to the discriminative methods.

We present artificial neural networks as a feasible replacement for a mechanistic model of mosquito abundance. We develop a feed-forward neural network, a long short-term memory recurrent neural network, and a gated recurrent unit network. We evaluate the networks in their ability to replicate the spatiotemporal features of mosquito populations predicted by the mechanistic model, and discuss how augmenting the training data with time series that emphasize specific dynamical behaviors affects model performance. We conclude with an outlook on how such equation-free models may facilitate vector control or the estimation of disease risk at arbitrary spatial scales.

The Hippo-YAP/TAZ signaling pathway plays a critical role in tissue homeostasis, tumorigenesis, and degeneration disorders. The regulation of YAP/TAZ levels is controlled by a complex regulatory network, where several feedback loops have been identified. However, it remains elusive how these feedback loops contain the YAP/TAZ levels and maintain the system in a healthy physiological state or trap the system in pathological conditions. Here, a mathematical model was developed to represent the YAP/TAZ regulatory network. Through theoretical analyses, three distinct states that designate the one physiological and two pathological outcomes were found. The transition from the physiological state to the two pathological states is mechanistically controlled by coupled bidirectional bistable switches, which are robust to parametric variation and stochastic fluctuations at the molecular level. This work provides a mechanistic understanding of the regulation and dysregulation of YAP/TAZ levels in tissue state transitions.

Over the past decade, biomarker discovery has become a key goal in psychiatry to aid in the more reliable diagnosis and prognosis of heterogeneous psychiatric conditions and the development of tailored therapies. Nevertheless, the prevailing statistical approach is still the mean group comparison between "cases" and "controls," which tends to ignore within-group variability. In this educational article, we used empirical data simulations to investigate how effect size, sample size, and the shape of distributions impact the interpretation of mean group differences for biomarker discovery. We then applied these statistical criteria to evaluate biomarker discovery in one area of psychiatric research-autism research. Across the most influential areas of autism research, effect size estimates ranged from small (d = 0.21, anatomical structure) to medium (d = 0.36 electrophysiology, d = 0.5, eye-tracking) to large (d = 1.1 theory of mind). We show that in normal distributions, this translates to approximately 45% to 63% of cases performing within 1 standard deviation (SD) of the typical range, i.e., they do not have a deficit/atypicality in a statistical sense. For a measure to have diagnostic utility as defined by 80% sensitivity and 80% specificity, Cohen's d of 1.66 is required, with still 40% of cases falling within 1 SD. However, in both normal and nonnormal distributions, 1 (skewness) or 2 (platykurtic, bimodal) biologically plausible subgroups may exist despite small or even nonsignificant mean group differences. This conclusion drastically contrasts the way mean group differences are frequently reported. Over 95% of studies omitted the "on average" when summarising their findings in their abstracts ("autistic people have deficits in X"), which can be misleading as it implies that the group-level difference applies to all individuals in that group. We outline practical approaches and steps for researchers to explore mean group comparisons for the discovery of stratification biomarkers.

Expression QTL (eQTL) analyses have suggested many genes mediating genome-wide association study (GWAS) signals but most GWAS signals still lack compelling explanatory genes. We have leveraged an adipose-specific gene regulatory network to infer expression regulator activities and phenotypic master regulators (MRs), which were used to detect activity QTLs (aQTLs) at cardiometabolic trait GWAS loci. Regulator activities were inferred with the VIPER algorithm that integrates enrichment of expected expression changes among a regulator's target genes with confidence in their regulator-target network interactions and target overlap between different regulators (i.e., pleiotropy). Phenotypic MRs were identified as those regulators whose activities were most important in predicting their respective phenotypes using random forest modeling. While eQTLs were typically more significant than aQTLs in cis, the opposite was true among candidate MRs in trans. Several GWAS loci colocalized with MR trans-eQTLs/aQTLs in the absence of colocalized cis-QTLs. Intriguingly, at the 1p36.1 BMI GWAS locus the EPHB2 cis-aQTL was stronger than its cis-eQTL and colocalized with the GWAS signal and 35 BMI MR trans-aQTLs, suggesting the GWAS signal may be mediated by effects on EPHB2 activity and its downstream effects on a network of BMI MRs. These MR and aQTL analyses represent systems genetic methods that may be broadly applied to supplement standard eQTL analyses for suggesting molecular effects mediating GWAS signals.

Construction of chromosomes 3D models based on single cell Hi-C data constitute an important challenge. We present a reconstruction approach, DPDchrom, that incorporates basic knowledge whether the reconstructed conformation should be coil-like or globular and spring relaxation at contact sites. In contrast to previously published protocols, DPDchrom can naturally form globular conformation due to the presence of explicit solvent. Benchmarking of this and several other methods on artificial polymer models reveals similar reconstruction accuracy at high contact density and DPDchrom advantage at low contact density. To compare 3D structures insensitively to spatial orientation and scale, we propose the Modified Jaccard Index. We analyzed two sources of the contact dropout: contact radius change and random contact sampling. We found that the reconstruction accuracy exponentially depends on the number of contacts per genomic bin allowing to estimate the reconstruction accuracy in advance. We applied DPDchrom to model chromosome configurations based on single-cell Hi-C data of mouse oocytes and found that these configurations differ significantly from a random one, that is consistent with other studies.

Severe acute respiratory coronavirus 2 (SARS-CoV-2), the causative agent of COVID-19, is of zoonotic origin. Evolutionary analyses assessing whether coronaviruses similar to SARS-CoV-2 infected ancestral species of modern-day animal hosts could be useful in identifying additional reservoirs of potentially dangerous coronaviruses. We reasoned that if a clade of species has been repeatedly exposed to a virus, then their proteins relevant for viral entry may exhibit adaptations that affect host susceptibility or response. We perform comparative analyses across the mammalian phylogeny of angiotensin-converting enzyme 2 (ACE2), the cellular receptor for SARS-CoV-2, in order to uncover evidence for selection acting at its binding interface with the SARS-CoV-2 spike protein. We uncover that in rodents there is evidence for adaptive amino acid substitutions at positions comprising the ACE2-spike interaction interface, whereas the variation within ACE2 proteins in primates and some other mammalian clades is not consistent with evolutionary adaptations. We also analyze aminopeptidase N (APN), the receptor for the human coronavirus 229E, a virus that causes the common cold, and find evidence for adaptation in primates. Altogether, our results suggest that the rodent and primate lineages may have had ancient exposures to viruses similar to SARS-CoV-2 and HCoV-229E, respectively.

Gene expression analysis is becoming increasingly utilized in neuro-immunology research, and there is a growing need for non-programming scientists to be able to analyze their own genomic data. MGEnrichment is a web application developed both to disseminate to the community our curated database of microglia-relevant gene lists, and to allow non-programming scientists to easily conduct statistical enrichment analysis on their gene expression data. Users can upload their own gene IDs to assess the relevance of their expression data against gene lists from other studies. We include example datasets of differentially expressed genes (DEGs) from human postmortem brain samples from Autism Spectrum Disorder (ASD) and matched controls. We demonstrate how MGEnrichment can be used to expand the interpretations of these DEG lists in terms of regulation of microglial gene expression and provide novel insights into how ASD DEGs may be implicated specifically in microglial development, microbiome responses and relationships to other neuropsychiatric disorders. This tool will be particularly useful for those working in microglia, autism spectrum disorders, and neuro-immune activation research. MGEnrichment is available at https://ciernialab.shinyapps.io/MGEnrichmentApp/ and further online documentation and datasets can be found at https://github.com/ciernialab/MGEnrichmentApp. The app is released under the GNU GPLv3 open source license.

Because local field potentials (LFPs) arise from multiple sources in different spatial locations, they do not easily reveal coordinated activity across neural populations on a trial-to-trial basis. As we show here, however, once disparate source signals are decoupled, their trial-to-trial fluctuations become more accessible, and cross-population correlations become more apparent. To decouple sources we introduce a general framework for estimation of current source densities (CSDs). In this framework, the set of LFPs result from noise being added to the transform of the CSD by a biophysical forward model, while the CSD is considered to be the sum of a zero-mean, stationary, spatiotemporal Gaussian process, having fast and slow components, and a mean function, which is the sum of multiple time-varying functions distributed across space, each varying across trials. We derived biophysical forward models relevant to the data we analyzed. In simulation studies this approach improved identification of source signals compared to existing CSD estimation methods. Using data recorded from primate auditory cortex, we analyzed trial-to-trial fluctuations in both steady-state and task-evoked signals. We found cortical layer-specific phase coupling between two probes and showed that the same analysis applied directly to LFPs did not recover these patterns. We also found task-evoked CSDs to be correlated across probes, at specific cortical depths. Using data from Neuropixels probes in mouse visual areas, we again found evidence for depth-specific phase coupling of primary visual cortex and lateromedial area based on the CSDs.

It is challenging to associate features such as human health outcomes, diet, environmental conditions, or other metadata to microbial community measurements, due in part to their quantitative properties. Microbiome multi-omics are typically noisy, sparse (zero-inflated), high-dimensional, extremely non-normal, and often in the form of count or compositional measurements. Here we introduce an optimized combination of novel and established methodology to assess multivariable association of microbial community features with complex metadata in population-scale observational studies. Our approach, MaAsLin 2 (Microbiome Multivariable Associations with Linear Models), uses generalized linear and mixed models to accommodate a wide variety of modern epidemiological studies, including cross-sectional and longitudinal designs, as well as a variety of data types (e.g., counts and relative abundances) with or without covariates and repeated measurements. To construct this method, we conducted a large-scale evaluation of a broad range of scenarios under which straightforward identification of meta-omics associations can be challenging. These simulation studies reveal that MaAsLin 2's linear model preserves statistical power in the presence of repeated measures and multiple covariates, while accounting for the nuances of meta-omics features and controlling false discovery. We also applied MaAsLin 2 to a microbial multi-omics dataset from the Integrative Human Microbiome (HMP2) project which, in addition to reproducing established results, revealed a unique, integrated landscape of inflammatory bowel diseases (IBD) across multiple time points and omics profiles.

Time lags in reporting to national surveillance systems represent a major barrier for the control of infectious diseases, preventing timely decision making and resource allocation. This issue is particularly acute for infectious diseases like malaria, which often impact rural and remote communities the hardest. In Guyana, a country located in South America, poor connectivity among remote malaria-endemic regions hampers surveillance efforts, making reporting delays a key challenge for elimination. Here, we analyze 13 years of malaria surveillance data, identifying key correlates of time lags between clinical cases occurring and being added to the central data system. We develop nowcasting methods that use historical patterns of reporting delays to estimate occurred-but-not-reported monthly malaria cases. To assess their performance, we implemented them retrospectively, using only information that would have been available at the time of estimation, and found that they substantially enhanced the estimates of malaria cases. Specifically, we found that the best performing models achieved up to two-fold improvements in accuracy (or error reduction) over known cases in selected regions. Our approach provides a simple, generalizable tool to improve malaria surveillance in endemic countries and is currently being implemented to help guide existing resource allocation and elimination efforts.

The cost of sequencing the genome is dropping at a much faster rate compared to assembling and finishing the genome. The use of lightly sampled genomes (genome-skims) could be transformative for genomic ecology, and results using k-mers have shown the advantage of this approach in identification and phylogenetic placement of eukaryotic species. Here, we revisit the basic question of estimating genomic parameters such as genome length, coverage, and repeat structure, focusing specifically on estimating the k-mer repeat spectrum. We show using a mix of theoretical and empirical analysis that there are fundamental limitations to estimating the k-mer spectra due to ill-conditioned systems, and that has implications for other genomic parameters. We get around this problem using a novel constrained optimization approach (Spline Linear Programming), where the constraints are learned empirically. On reads simulated at 1X coverage from 66 genomes, our method, REPeat SPECTra Estimation (RESPECT), had 2.2% error in length estimation compared to 27% error previously achieved. In shotgun sequenced read samples with contaminants, RESPECT length estimates had median error 4%, in contrast to other methods that had median error 80%. Together, the results suggest that low-pass genomic sequencing can yield reliable estimates of the length and repeat content of the genome. The RESPECT software will be publicly available at https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_shahab-2Dsarmashghi_RESPECT.git&d=DwIGAw&c=-35OiAkTchMrZOngvJPOeA&r=ZozViWvD1E8PorCkfwYKYQMVKFoEcqLFm4Tg49XnPcA&m=f-xS8GMHKckknkc7Xpp8FJYw_ltUwz5frOw1a5pJ81EpdTOK8xhbYmrN4ZxniM96&s=717o8hLR1JmHFpRPSWG6xdUQTikyUjicjkipjFsKG4w&e=.

To obtain a coherent perception of the world, our senses need to be in alignment. When we encounter misaligned cues from two sensory modalities, the brain must infer which cue is faulty and recalibrate the corresponding sense. We examined whether and how the brain uses cue reliability to identify the miscalibrated sense by measuring the audiovisual ventriloquism aftereffect for stimuli of varying visual reliability. To adjust for modality-specific biases, visual stimulus locations were chosen based on perceived alignment with auditory stimulus locations for each participant. During an audiovisual recalibration phase, participants were presented with bimodal stimuli with a fixed perceptual spatial discrepancy; they localized one modality, cued after stimulus presentation. Unimodal auditory and visual localization was measured before and after the audiovisual recalibration phase. We compared participants' behavior to the predictions of three models of recalibration: (a) Reliability-based: each modality is recalibrated based on its relative reliability-less reliable cues are recalibrated more; (b) Fixed-ratio: the degree of recalibration for each modality is fixed; (c) Causal-inference: recalibration is directly determined by the discrepancy between a cue and its estimate, which in turn depends on the reliability of both cues, and inference about how likely the two cues derive from a common source. Vision was hardly recalibrated by audition. Auditory recalibration by vision changed idiosyncratically as visual reliability decreased: the extent of auditory recalibration either decreased monotonically, peaked at medium visual reliability, or increased monotonically. The latter two patterns cannot be explained by either the reliability-based or fixed-ratio models. Only the causal-inference model of recalibration captures the idiosyncratic influences of cue reliability on recalibration. We conclude that cue reliability, causal inference, and modality-specific biases guide cross-modal recalibration indirectly by determining the perception of audiovisual stimuli.

In many models of evolving populations, genetic drift has an outsized role relative to natural selection, or vice versa. While there are many scenarios in which one of these two assumptions is reasonable, intermediate balances between these forces are also biologically relevant. In this study, we consider some natural axioms for modeling intermediate selection intensities, and we explore how to quantify the long-term evolutionary dynamics of such a process. To illustrate the sensitivity of evolutionary dynamics to drift and selection, we show that there can be a "sweet spot" for the balance of these two forces, with sufficient noise for rare mutants to become established and sufficient selection to spread. This balance allows prosocial traits to evolve in evolutionary models that were previously thought to be unconducive to the emergence and spread of altruistic behaviors. Furthermore, the effects of selection intensity on long-run evolutionary outcomes in these settings, such as when there is global competition for reproduction, can be highly non-monotonic. Although intermediate selection intensities (neither weak nor strong) are notoriously difficult to study analytically, they are often biologically relevant; and the results we report suggest that they can elicit novel and rich dynamics in the evolution of prosocial behaviors.

Riboswitches are naturally occurring RNA elements that control bacterial gene expression by binding to specific small molecules. They serve as important models for RNA-small molecule recognition and have also become a novel class of targets for developing antibiotics. Here, we carried out conventional and enhanced-sampling molecular dynamics (MD) simulations, totaling 153.5 μs, to characterize the determinants of binding free energies and unbinding paths for the cognate and synthetic ligands of a PreQ1 riboswitch. Binding free energy analysis showed that two triplets of nucleotides, U6-C15-A29 and G5-G11-C16, contribute the most to the binding of the cognate ligands, by hydrogen bonding and by base stacking, respectively. Mg2+ ions are essential in stabilizing the binding pocket. For the synthetic ligands, the hydrogen-bonding contributions of the U6-C15-A29 triplet are significantly compromised, and the bound state resembles the apo state in several respects, including the disengagement of the C15-A14-A13 and A32-G33 base stacks. The bulkier synthetic ligands lead to significantly loosening of the binding pocket, including extrusion of the C15 nucleobase and a widening of the C15-C30 groove. Enhanced-sampling simulations further revealed that the cognate and synthetic ligands unbind in almost opposite directions. Our work offers new insight for designing riboswitch ligands.

GABA (γ-amino butyric acid) is an inhibitory neurotransmitter in the adult brain that can mediate depolarizing responses during development or after neuropathological insults. Under which conditions GABAergic membrane depolarizations are sufficient to impose excitatory effects is hard to predict, as shunting inhibition and GABAergic effects on spatiotemporal filtering of excitatory inputs must be considered. To evaluate at which reversal potential a net excitatory effect was imposed by GABA (EGABAThr), we performed a detailed in-silico study using simple neuronal topologies and distinct spatiotemporal relations between GABAergic and glutamatergic inputs. These simulations revealed for GABAergic synapses located at the soma an EGABAThr close to action potential threshold (EAPThr), while with increasing dendritic distance EGABAThr shifted to positive values. The impact of GABA on AMPA-mediated inputs revealed a complex temporal and spatial dependency. EGABAThr depends on the temporal relation between GABA and AMPA inputs, with a striking negative shift in EGABAThr for AMPA inputs appearing after the GABA input. The spatial dependency between GABA and AMPA inputs revealed a complex profile, with EGABAThr being shifted to values negative to EAPThr for AMPA synapses located proximally to the GABA input, while for distally located AMPA synapses the dendritic distance had only a minor effect on EGABAThr. For tonic GABAergic conductances EGABAThr was negative to EAPThr over a wide range of gGABAtonic values. In summary, these results demonstrate that for several physiologically relevant situations EGABAThr is negative to EAPThr, suggesting that depolarizing GABAergic responses can mediate excitatory effects even if EGABA did not reach EAPThr.

Sudden changes in visual scenes often indicate important events for behavior. For their quick and reliable detection, the brain must be capable to process these changes as independently as possible from its current activation state. In motion-selective area MT, neurons respond to instantaneous speed changes with pronounced transients, often far exceeding the expected response as derived from their speed tuning profile. We here show that this complex, non-linear behavior emerges from the combined temporal dynamics of excitation and divisive inhibition, and provide a comprehensive mathematical analysis. A central prediction derived from this investigation is that attention increases the steepness of the transient response irrespective of the activation state prior to a stimulus change, and irrespective of the sign of the change (i.e. irrespective of whether the stimulus is accelerating or decelerating). Extracellular recordings of attention-dependent representation of both speed increments and decrements confirmed this prediction and suggest that improved change detection derives from basic computations in a canonical cortical circuitry.

Emergent response properties of sensory neurons depend on circuit connectivity and somatodendritic processing. Neurons of the barn owl's external nucleus of the inferior colliculus (ICx) display emergence of spatial selectivity. These neurons use interaural time difference (ITD) as a cue for the horizontal direction of sound sources. ITD is detected by upstream brainstem neurons with narrow frequency tuning, resulting in spatially ambiguous responses. This spatial ambiguity is resolved by ICx neurons integrating inputs over frequency, a relevant processing in sound localization across species. Previous models have predicted that ICx neurons function as point neurons that linearly integrate inputs across frequency. However, the complex dendritic trees and spines of ICx neurons raises the question of whether this prediction is accurate. Data from in vivo intracellular recordings of ICx neurons were used to address this question. Results revealed diverse frequency integration properties, where some ICx neurons showed responses consistent with the point neuron hypothesis and others with nonlinear dendritic integration. Modeling showed that varied connectivity patterns and forms of dendritic processing may underlie observed ICx neurons' frequency integration processing. These results corroborate the ability of neurons with complex dendritic trees to implement diverse linear and nonlinear integration of synaptic inputs, of relevance for adaptive coding and learning, and supporting a fundamental mechanism in sound localization.

Network propagation refers to a class of algorithms that integrate information from input data across connected nodes in a given network. These algorithms have wide applications in systems biology, protein function prediction, inferring condition-specifically altered sub-networks, and prioritizing disease genes. Despite the popularity of network propagation, there is a lack of comparative analyses of different algorithms on real data and little guidance on how to select and parameterize the various algorithms. Here, we address this problem by analyzing different combinations of network normalization and propagation methods and by demonstrating schemes for the identification of optimal parameter settings on real proteome and transcriptome data. Our work highlights the risk of a 'topology bias' caused by the incorrect use of network normalization approaches. Capitalizing on the fact that network propagation is a regularization approach, we show that minimizing the bias-variance tradeoff can be utilized for selecting optimal parameters. The application to real multi-omics data demonstrated that optimal parameters could also be obtained by either maximizing the agreement between different omics layers (e.g. proteome and transcriptome) or by maximizing the consistency between biological replicates. Furthermore, we exemplified the utility and robustness of network propagation on multi-omics datasets for identifying ageing-associated genes in brain and liver tissues of rats and for elucidating molecular mechanisms underlying prostate cancer progression. Overall, this work compares different network propagation approaches and it presents strategies for how to use network propagation algorithms to optimally address a specific research question at hand.

Computational biology has gained traction as an independent scientific discipline over the last years in South America. However, there is still a growing need for bioscientists, from different backgrounds, with different levels, to acquire programming skills, which could reduce the time from data to insights and bridge communication between life scientists and computer scientists. Python is a programming language extensively used in bioinformatics and data science, which is particularly suitable for beginners. Here, we describe the conception, organization, and implementation of the Brazilian Python Workshop for Biological Data. This workshop has been organized by graduate and undergraduate students and supported, mostly in administrative matters, by experienced faculty members since 2017. The workshop was conceived for teaching bioscientists, mainly students in Brazil, on how to program in a biological context. The goal of this article was to share our experience with the 2020 edition of the workshop in its virtual format due to the Coronavirus Disease 2019 (COVID-19) pandemic and to compare and contrast this year's experience with the previous in-person editions. We described a hands-on and live coding workshop model for teaching introductory Python programming. We also highlighted the adaptations made from in-person to online format in 2020, the participants' assessment of learning progression, and general workshop management. Lastly, we provided a summary and reflections from our personal experiences from the workshops of the last 4 years. Our takeaways included the benefits of the learning from learners' feedback (LLF) that allowed us to improve the workshop in real time, in the short, and likely in the long term. We concluded that the Brazilian Python Workshop for Biological Data is a highly effective workshop model for teaching a programming language that allows bioscientists to go beyond an initial exploration of programming skills for data analysis in the medium to long term.

Although osteosarcoma (OS) is a rare cancer, it is the most common primary malignant bone tumor in children and adolescents. BRCAness is a phenotypical trait in tumors with a defect in homologous recombination repair, resembling tumors with inactivation of BRCA1/2, rendering these tumors sensitive to poly (ADP)-ribose polymerase inhibitors (PARPi). Recently, OS was shown to exhibit molecular features of BRCAness. Our goal was to develop a method complementing existing genomic methods to aid clinical decision making on administering PARPi in OS patients. OS samples with DNA-methylation data were divided to BRCAness-positive and negative groups based on the degree of their genomic instability (n = 41). Methylation probes were ranked according to decreasing variance difference between two groups. The top 2000 probes were selected for training and cross-validation of the random forest algorithm. Two-thirds of available OS RNA-Seq samples (n = 17) from the top and bottom of the sample list ranked according to genome instability score were subjected to differential expression and, subsequently, to gene set enrichment analysis (GSEA). The combined accuracy of trained random forest was 85% and the average area under the ROC curve (AUC) was 0.95. There were 449 upregulated and 1,079 downregulated genes in the BRCAness-positive group (fdr < 0.05). GSEA of upregulated genes detected enrichment of DNA replication and mismatch repair and homologous recombination signatures (FWER < 0.05). Validation of the BRCAness classifier with an independent OS set (n = 20) collected later in the course of study showed AUC of 0.87 with an accuracy of 90%. GSEA signatures computed for this test set were matching the ones observed in the training set enrichment analysis. In conclusion, we developed a new classifier based on DNA-methylation patterns that detects BRCAness in OS samples with high accuracy. GSEA identified genome instability signatures. Machine-learning and gene expression approaches add new epigenomic and transcriptomic aspects to already established genomic methods for evaluation of BRCAness in osteosarcoma and can be extended to cancers characterized by genome instability.

Isogenic cells cultured together show heterogeneity in their proliferation rate. To determine the differences between fast and slow-proliferating cells, we developed a method to sort cells by proliferation rate, and performed RNA-seq on slow and fast proliferating subpopulations of pluripotent mouse embryonic stem cells (mESCs) and mouse fibroblasts. We found that slowly proliferating mESCs have a more naïve pluripotent character. We identified an evolutionarily conserved proliferation-correlated transcriptomic signature that is common to all eukaryotes: fast cells have higher expression of genes for protein synthesis and protein degradation. This signature accurately predicted growth rate in yeast and cancer cells, and identified lineage-specific proliferation dynamics during development, using C. elegans scRNA-seq data. In contrast, sorting by mitochondria membrane potential revealed a highly cell-type specific mitochondria-state related transcriptome. mESCs with hyperpolarized mitochondria are fast proliferating, while the opposite is true for fibroblasts. The mitochondrial electron transport chain inhibitor antimycin affected slow and fast subpopulations differently. While a major transcriptional-signature associated with cell-to-cell heterogeneity in proliferation is conserved, the metabolic and energetic dependency of cell proliferation is cell-type specific.

Functional, usable, and maintainable open-source software is increasingly essential to scientific research, but there is a large variation in formal training for software development and maintainability. Here, we propose 10 "rules" centered on 2 best practice components: clean code and testing. These 2 areas are relatively straightforward and provide substantial utility relative to the learning investment. Adopting clean code practices helps to standardize and organize software code in order to enhance readability and reduce cognitive load for both the initial developer and subsequent contributors; this allows developers to concentrate on core functionality and reduce errors. Clean coding styles make software code more amenable to testing, including unit tests that work best with modular and consistent software code. Unit tests interrogate specific and isolated coding behavior to reduce coding errors and ensure intended functionality, especially as code increases in complexity; unit tests also implicitly provide example usages of code. Other forms of testing are geared to discover erroneous behavior arising from unexpected inputs or emerging from the interaction of complex codebases. Although conforming to coding styles and designing tests can add time to the software development project in the short term, these foundational tools can help to improve the correctness, quality, usability, and maintainability of open-source scientific software code. They also advance the principal point of scientific research: producing accurate results in a reproducible way. In addition to suggesting several tips for getting started with clean code and testing practices, we recommend numerous tools for the popular open-source scientific software languages Python, R, and Julia.

The growing number of next-generation sequencing (NGS) data presents a unique opportunity to study the combined impact of mitochondrial and nuclear-encoded genetic variation in complex disease. Mitochondrial DNA variants and in particular, heteroplasmic variants, are critical for determining human disease severity. While there are approaches for obtaining mitochondrial DNA variants from NGS data, these software do not account for the unique characteristics of mitochondrial genetics and can be inaccurate even for homoplasmic variants. We introduce MitoScape, a novel, big-data, software for extracting mitochondrial DNA sequences from NGS. MitoScape adopts a novel departure from other algorithms by using machine learning to model the unique characteristics of mitochondrial genetics. We also employ a novel approach of using rho-zero (mitochondrial DNA-depleted) data to model nuclear-encoded mitochondrial sequences. We showed that MitoScape produces accurate heteroplasmy estimates using gold-standard mitochondrial DNA data. We provide a comprehensive comparison of the most common tools for obtaining mtDNA variants from NGS and showed that MitoScape had superior performance to compared tools in every statistically category we compared, including false positives and false negatives. By applying MitoScape to common disease examples, we illustrate how MitoScape facilitates important heteroplasmy-disease association discoveries by expanding upon a reported association between hypertrophic cardiomyopathy and mitochondrial haplogroup T in men (adjusted p-value = 0.003). The improved accuracy of mitochondrial DNA variants produced by MitoScape will be instrumental in diagnosing disease in the context of personalized medicine and clinical diagnostics.

Cell assemblies are thought to be the substrate of memory in the brain. Theoretical studies have previously shown that assemblies can be formed in networks with multiple types of plasticity. But how exactly they are formed and how they encode information is yet to be fully understood. One possibility is that memories are stored in silent assemblies. Here we used a computational model to study the formation of silent assemblies in a network of spiking neurons with excitatory and inhibitory plasticity. We found that even though the formed assemblies were silent in terms of mean firing rate, they had an increased coefficient of variation of inter-spike intervals. We also found that this spiking irregularity could be read out with support of short-term plasticity, and that it could contribute to the longevity of memories.

Genome-scale metabolic models (GEMs) provide a powerful framework for simulating the entire set of biochemical reactions in a cell using a constraint-based modeling strategy called flux balance analysis (FBA). FBA relies on an assumed metabolic objective for generating metabolic fluxes using GEMs. But, the most appropriate metabolic objective is not always obvious for a given condition and is likely context-specific, which often complicate the estimation of metabolic flux alterations between conditions. Here, we propose a new method, called ΔFBA (deltaFBA), that integrates differential gene expression data to evaluate directly metabolic flux differences between two conditions. Notably, ΔFBA does not require specifying the cellular objective. Rather, ΔFBA seeks to maximize the consistency and minimize inconsistency between the predicted flux differences and differential gene expression. We showcased the performance of ΔFBA through several case studies involving the prediction of metabolic alterations caused by genetic and environmental perturbations in Escherichia coli and caused by Type-2 diabetes in human muscle. Importantly, in comparison to existing methods, ΔFBA gives a more accurate prediction of flux differences.

The first stage of the metastatic cascade often involves motile cells emerging from a primary tumor either as single cells or as clusters. These cells enter the circulation, transit to other parts of the body and finally are responsible for growth of secondary tumors in distant organs. The mode of dissemination is believed to depend on the EMT nature (epithelial, hybrid or mesenchymal) of the cells. Here, we calculate the cluster size distribution of these migrating cells, using a mechanistic computational model, in presence of different degree of EMT-ness of the cells; EMT is treated as given rise to changes in their active motile forces (μ) and cell-medium surface tension (Γ). We find that, for (μ > μmin, Γ > 1), when the cells are hybrid in nature, the mean cluster size, [Formula: see text], where μmin increases with increase in Γ. For Γ ≤ 0, [Formula: see text], the cells behave as completely mesenchymal. In presence of spectrum of hybrid states with different degree of EMT-ness (motility) in primary tumor, the cells which are relatively more mesenchymal (higher μ) in nature, form larger clusters, whereas the smaller clusters are relatively more epithelial (lower μ). Moreover, the heterogeneity in μ is comparatively higher for smaller clusters with respect to that for larger clusters. We also observe that more extended cell shapes promote the formation of smaller clusters. Overall, this study establishes a framework which connects the nature and size of migrating clusters disseminating from a primary tumor with the phenotypic composition of the tumor, and can lead to the better understanding of metastasis.

A common feature of morphogenesis is the formation of three-dimensional structures from the folding of two-dimensional epithelial sheets, aided by cell shape changes at the cellular-level. Changes in cell shape must be studied in the context of cell-polarised biomechanical processes within the epithelial sheet. In epithelia with highly curved surfaces, finding single-cell alignment along a biological axis can be difficult to automate in silico. We present 'Origami', a MATLAB-based image analysis pipeline to compute direction-variant cell shape features along the epithelial apico-basal axis. Our automated method accurately computed direction vectors denoting the apico-basal axis in regions with opposing curvature in synthetic epithelia and fluorescence images of zebrafish embryos. As proof of concept, we identified different cell shape signatures in the developing zebrafish inner ear, where the epithelium deforms in opposite orientations to form different structures. Origami is designed to be user-friendly and is generally applicable to fluorescence images of curved epithelia.

The spleen, the largest secondary lymphoid organ in humans, not only fulfils a broad range of immune functions, but also plays an important role in red blood cell's (RBC) life cycle. Although much progress has been made to elucidate the critical biological processes involved in the maturation of young RBCs (reticulocytes) as well as removal of senescent RBCs in the spleen, the underlying mechanisms driving these processes are still obscure. Herein, we perform a computational study to simulate the passage of RBCs through interendothelial slits (IES) in the spleen at different stages of their lifespan and investigate the role of the spleen in facilitating the maturation of reticulocytes and in clearing the senescent RBCs. Our simulations reveal that at the beginning of the RBC life cycle, intracellular non-deformable particles in reticulocytes can be biomechanically expelled from the cell upon passage through IES, an insightful explanation of why this peculiar "pitting" process is spleen-specific. Our results also show that immature RBCs shed surface area by releasing vesicles after crossing IES and progressively acquire the biconcave shape of mature RBCs. These findings likely explain why RBCs from splenectomized patients are significantly larger than those from nonsplenectomized subjects. Finally, we show that at the end of their life span, senescent RBCs are not only retained by IES due to reduced deformability but also become susceptible to mechanical lysis under shear stress. This finding supports the recent hypothesis that transformation into a hemolyzed ghost is a prerequisite for phagocytosis of senescent RBCs. Altogether, our computational investigation illustrates critical biological processes in the spleen that cannot be observed in vivo or in vitro and offer insights into the role of the spleen in the RBC physiology.

The study of microbial communities and their interactions has attracted the interest of the scientific community, because of their potential for applications in biotechnology, ecology and medicine. The complexity of interspecies interactions, which are key for the macroscopic behavior of microbial communities, cannot be studied easily experimentally. For this reason, the modeling of microbial communities has begun to leverage the knowledge of established constraint-based methods, which have long been used for studying and analyzing the microbial metabolism of individual species based on genome-scale metabolic reconstructions of microorganisms. A main problem of genome-scale metabolic reconstructions is that they usually contain metabolic gaps due to genome misannotations and unknown enzyme functions. This problem is traditionally solved by using gap-filling algorithms that add biochemical reactions from external databases to the metabolic reconstruction, in order to restore model growth. However, gap-filling algorithms could evolve by taking into account metabolic interactions among species that coexist in microbial communities. In this work, a gap-filling method that resolves metabolic gaps at the community level was developed. The efficacy of the algorithm was tested by analyzing its ability to resolve metabolic gaps on a synthetic community of auxotrophic Escherichia coli strains. Subsequently, the algorithm was applied to resolve metabolic gaps and predict metabolic interactions in a community of Bifidobacterium adolescentis and Faecalibacterium prausnitzii, two species present in the human gut microbiota, and in an experimentally studied community of Dehalobacter and Bacteroidales species of the ACT-3 community. The community gap-filling method can facilitate the improvement of metabolic models and the identification of metabolic interactions that are difficult to identify experimentally in microbial communities.

In biology, we are often confronted with information-rich, large-scale trajectory data, but exploring and communicating patterns in such data can be a cumbersome task. Ideally, the data should be wrapped with an interactive visualisation in one concise packet that makes it straightforward to create and test hypotheses collaboratively. To address these challenges, we have developed a tool, linus, which makes the process of exploring and sharing 3D trajectories as easy as browsing a website. We provide a python script that reads trajectory data, enriches them with additional features such as edge bundling or custom axes, and generates an interactive web-based visualisation that can be shared online. linus facilitates the collaborative discovery of patterns in complex trajectory data.

Sensory information from different modalities is processed in parallel, and then integrated in associative brain areas to improve object identification and the interpretation of sensory experiences. The Superior Colliculus (SC) is a midbrain structure that plays a critical role in integrating visual, auditory, and somatosensory input to assess saliency and promote action. Although the response properties of the individual SC neurons to visuoauditory stimuli have been characterized, little is known about the spatial and temporal dynamics of the integration at the population level. Here we recorded the response properties of SC neurons to spatially restricted visual and auditory stimuli using large-scale electrophysiology. We then created a general, population-level model that explains the spatial, temporal, and intensity requirements of stimuli needed for sensory integration. We found that the mouse SC contains topographically organized visual and auditory neurons that exhibit nonlinear multisensory integration. We show that nonlinear integration depends on properties of auditory but not visual stimuli. We also find that a heuristically derived nonlinear modulation function reveals conditions required for sensory integration that are consistent with previously proposed models of sensory integration such as spatial matching and the principle of inverse effectiveness.

Regulation of cytosolic calcium (Ca2+) dynamics is fundamental to microglial function. Temporal and spatial Ca2+ fluxes are induced from a complicated signal transduction pathway linked to brain ionic homeostasis. In this paper, we develop a novel biophysical model of Ca2+ and sodium (Na+) dynamics in human microglia and evaluate the contribution of purinergic receptors (P2XRs) to both intracellular Ca2+ and Na+ levels in response to agonist/ATP binding. This is the first comprehensive model that integrates P2XRs to predict intricate Ca2+ and Na+ transient responses in microglia. Specifically, a novel compact biophysical model is proposed for the capture of whole-cell patch-clamp currents associated with P2X4 and P2X7 receptors, which is composed of only four state variables. The entire model shows that intricate intracellular ion dynamics arise from the coupled interaction between P2X4 and P2X7 receptors, the Na+/Ca2+ exchanger (NCX), Ca2+ extrusion by the plasma membrane Ca2+ ATPase (PMCA), and Ca2+ and Na+ leak channels. Both P2XRs are modelled as two separate adenosine triphosphate (ATP) gated Ca2+ and Na+ conductance channels, where the stoichiometry is the removal of one Ca2+ for the hydrolysis of one ATP molecule. Two unique sets of model parameters were determined using an evolutionary algorithm to optimise fitting to experimental data for each of the receptors. This allows the proposed model to capture both human P2X7 and P2X4 data (hP2X7 and hP2X4). The model architecture enables a high degree of simplicity, accuracy and predictability of Ca2+ and Na+ dynamics thus providing quantitative insights into different behaviours of intracellular Na+ and Ca2+ which will guide future experimental research. Understanding the interactions between these receptors and other membrane-bound transporters provides a step forward in resolving the qualitative link between purinergic receptors and microglial physiology and their contribution to brain pathology.

Information flow within and between cells depends significantly on calcium (Ca2+) signaling dynamics. However, the biophysical mechanisms that govern emergent patterns of Ca2+ signaling dynamics at the organ level remain elusive. Recent experimental studies in developing Drosophila wing imaginal discs demonstrate the emergence of four distinct patterns of Ca2+ activity: Ca2+ spikes, intercellular Ca2+ transients, tissue-level Ca2+ waves, and a global "fluttering" state. Here, we used a combination of computational modeling and experimental approaches to identify two different populations of cells within tissues that are connected by gap junction proteins. We term these two subpopulations "initiator cells," defined by elevated levels of Phospholipase C (PLC) activity, and "standby cells," which exhibit baseline activity. We found that the type and strength of hormonal stimulation and extent of gap junctional communication jointly determine the predominate class of Ca2+ signaling activity. Further, single-cell Ca2+ spikes are stimulated by insulin, while intercellular Ca2+ waves depend on Gαq activity. Our computational model successfully reproduces how the dynamics of Ca2+ transients varies during organ growth. Phenotypic analysis of perturbations to Gαq and insulin signaling support an integrated model of cytoplasmic Ca2+ as a dynamic reporter of overall tissue growth. Further, we show that perturbations to Ca2+ signaling tune the final size of organs. This work provides a platform to further study how organ size regulation emerges from the crosstalk between biochemical growth signals and heterogeneous cell signaling states.