Cell-cell interactions are vital for numerous biological processes including development, differentiation, and response to inflammation. Currently, most methods for studying interactions on scRNA-seq level are based on curated databases of ligands and receptors. While those methods are useful, they are limited to our current biological knowledge. Recent advances in single cell protocols have allowed for physically interacting cells to be captured, and as such we have the potential to study interactions in a complemantary way without relying on prior knowledge. We introduce a new method based on Latent Dirichlet Allocation (LDA) for detecting genes that change as a result of interaction. We apply our method to synthetic datasets to demonstrate its ability to detect genes that change in an interacting population compared to a reference population. Next, we apply our approach to two datasets of physically interacting cells to identify the genes that change as a result of interaction, examples include adhesion and co-stimulatory molecules which confirm physical interaction between cells. For each dataset we produce a ranking of genes that are changing in subpopulations of the interacting cells. In addition to the genes discussed in the original publications, we highlight further candidates for interaction in the top 100 and 300 ranked genes. Lastly, we apply our method to a dataset generated by a standard droplet-based protocol not designed to capture interacting cells, and discuss its suitability for analysing interactions. We present a method that streamlines detection of interactions and does not require prior clustering and generation of synthetic reference profiles to detect changes in expression.

Population size has long been considered an important driver of cultural diversity and complexity. Results from population genetics, however, demonstrate that in populations with complex demographic structure or mode of inheritance, it is not the census population size, N, but the effective size of a population, Ne, that determines important evolutionary parameters. Here, we examine the concept of effective population size for traits that evolve culturally, through processes of innovation and social learning. We use mathematical and computational modeling approaches to investigate how cultural Ne and levels of diversity depend on (1) the way traits are learned, (2) population connectedness, and (3) social network structure. We show that one-to-many and frequency-dependent transmission can temporally or permanently lower effective population size compared to census numbers. We caution that migration and cultural exchange can have counter-intuitive effects on Ne. Network density in random networks leaves Ne unchanged, scale-free networks tend to decrease and small-world networks tend to increase Ne compared to census numbers. For one-to-many transmission and different network structures, larger effective sizes are closely associated with higher cultural diversity. For connectedness, however, even small amounts of migration and cultural exchange result in high diversity independently of Ne. Extending previous work, our results highlight the importance of carefully defining effective population size for cultural systems and show that inferring Ne requires detailed knowledge about underlying cultural and demographic processes.

Red blood cells can withstand the harsh mechanical conditions in the vasculature only because the bending rigidity of their plasma membrane is complemented by the shear elasticity of the underlying spectrin-actin network. During an infection by the malaria parasite Plasmodium falciparum, the parasite mines host actin from the junctional complexes and establishes a system of adhesive knobs, whose main structural component is the knob-associated histidine rich protein (KAHRP) secreted by the parasite. Here we aim at a mechanistic understanding of this dramatic transformation process. We have developed a particle-based computational model for the cytoskeleton of red blood cells and simulated it with Brownian dynamics to predict the mechanical changes resulting from actin mining and KAHRP-clustering. Our simulations include the three-dimensional conformations of the semi-flexible spectrin chains, the capping of the actin protofilaments and several established binding sites for KAHRP. For the healthy red blood cell, we find that incorporation of actin protofilaments leads to two regimes in the shear response. Actin mining decreases the shear modulus, but knob formation increases it. We show that dynamical changes in KAHRP binding affinities can explain the experimentally observed relocalization of KAHRP from ankyrin to actin complexes and demonstrate good qualitative agreement with experiments by measuring pair cross-correlations both in the computer simulations and in super-resolution imaging experiments.

Cellular barcoding methods offer the exciting possibility of 'infinite-pseudocolor' anatomical reconstruction-i.e., assigning each neuron its own random unique barcoded 'pseudocolor,' and then using these pseudocolors to trace the microanatomy of each neuron. Here we use simulations, based on densely-reconstructed electron microscopy microanatomy, with signal structure matched to real barcoding data, to quantify the feasibility of this procedure. We develop a new blind demixing approach to recover the barcodes that label each neuron, and validate this method on real data with known barcodes. We also develop a neural network which uses the recovered barcodes to reconstruct the neuronal morphology from the observed fluorescence imaging data, 'connecting the dots' between discontiguous barcode amplicon signals. We find that accurate recovery should be feasible, provided that the barcode signal density is sufficiently high. This study suggests the possibility of mapping the morphology and projection pattern of many individual neurons simultaneously, at high resolution and at large scale, via conventional light microscopy.

Robust and accurate behavioral tracking is essential for ethological studies. Common methods for tracking and extracting behavior rely on user adjusted heuristics that can significantly vary across different individuals, environments, and experimental conditions. As a result, they are difficult to implement in large-scale behavioral studies with complex, heterogenous environmental conditions. Recently developed deep-learning methods for object recognition such as Faster R-CNN have advantages in their speed, accuracy, and robustness. Here, we show that Faster R-CNN can be employed for identification and detection of Caenorhabditis elegans in a variety of life stages in complex environments. We applied the algorithm to track animal speeds during development, fecundity rates and spatial distribution in reproductive adults, and behavioral decline in aging populations. By doing so, we demonstrate the flexibility, speed, and scalability of Faster R-CNN across a variety of experimental conditions, illustrating its generalized use for future large-scale behavioral studies.

Many pathogenic missense mutations are found in protein positions that are neither well-conserved nor fall in any known functional domains. Consequently, we lack any mechanistic underpinning of dysfunction caused by such mutations. We explored the disruption of allosteric dynamic coupling between these positions and the known functional sites as a possible mechanism for pathogenesis. In this study, we present an analysis of 591 pathogenic missense variants in 144 human enzymes that suggests that allosteric dynamic coupling of mutated positions with known active sites is a plausible biophysical mechanism and evidence of their functional importance. We illustrate this mechanism in a case study of β-Glucocerebrosidase (GCase) in which a vast majority of 94 sites harboring Gaucher disease-associated missense variants are located some distance away from the active site. An analysis of the conformational dynamics of GCase suggests that mutations on these distal sites cause changes in the flexibility of active site residues despite their distance, indicating a dynamic communication network throughout the protein. The disruption of the long-distance dynamic coupling caused by missense mutations may provide a plausible general mechanistic explanation for biological dysfunction and disease.

Cortical actin networks are highly dynamic and play critical roles in shaping the mechanical properties of cells. The actin cytoskeleton undergoes significant reorganization in many different contexts, including during directed cell migration and over the course of the cell cycle, when cortical actin can transition between different configurations such as open patched meshworks, homogeneous distributions, and aligned bundles. Several types of myosin motor proteins, characterized by different kinetic parameters, have been involved in this reorganization of actin filaments. Given the limitations in studying the interactions of actin with myosin in vivo, we propose stochastic agent-based models and develop a set of data analysis measures to assess how myosin motor proteins mediate various actin organizations. In particular, we identify individual motor parameters, such as motor binding rate and step size, that generate actin networks with different levels of contractility and different patterns of myosin motor localization, which have previously been observed experimentally. In simulations where two motor populations with distinct kinetic parameters interact with the same actin network, we find that motors may act in a complementary way, by tuning the actin network organization, or in an antagonistic way, where one motor emerges as dominant. This modeling and data analysis framework also uncovers parameter regimes where spatial segregation between motor populations is achieved. By allowing for changes in kinetic rates during the actin-myosin dynamic simulations, our work suggests that certain actin-myosin organizations may require additional regulation beyond mediation by motor proteins in order to reconfigure the cytoskeleton network on experimentally-observed timescales.

We evaluate the efficiency of various heuristic strategies for allocating vaccines against COVID-19 and compare them to strategies found using optimal control theory. Our approach is based on a mathematical model which tracks the spread of disease among different age groups and across different geographical regions, and we introduce a method to combine age-specific contact data to geographical movement data. As a case study, we model the epidemic in the population of mainland Finland utilizing mobility data from a major telecom operator. Our approach allows to determine which geographical regions and age groups should be targeted first in order to minimize the number of deaths. In the scenarios that we test, we find that distributing vaccines demographically and in an age-descending order is not optimal for minimizing deaths and the burden of disease. Instead, more lives could be saved by using strategies which emphasize high-incidence regions and distribute vaccines in parallel to multiple age groups. The level of emphasis that high-incidence regions should be given depends on the overall transmission rate in the population. This observation highlights the importance of updating the vaccination strategy when the effective reproduction number changes due to the general contact patterns changing and new virus variants entering.

Broadly neutralizing antibodies (bNAbs) are promising agents to prevent HIV infection and achieve HIV remission without antiretroviral therapy (ART). As with ART, bNAb combinations are likely needed to cover HIV's extensive diversity. Not all bNAbs are identical in terms of their breadth, potency, and in vivo longevity (half-life). Given these differences, it is important to optimally select the composition, or dose ratio, of combination bNAb therapies for future clinical studies. We developed a model that synthesizes 1) pharmacokinetics, 2) potency against a wide HIV diversity, 3) interaction models for how drugs work together, and 4) correlates that translate in vitro potency to clinical protection. We found optimization requires drug-specific balances between potency, longevity, and interaction type. As an example, tradeoffs between longevity and potency are shown by comparing a combination therapy to a bi-specific antibody (a single protein merging both bNAbs) that takes the better potency but the worse longevity of the two components. Then, we illustrate a realistic dose ratio optimization of a triple combination of VRC07, 3BNC117, and 10-1074 bNAbs. We apply protection estimates derived from both a non-human primate (NHP) challenge study meta-analysis and the human antibody mediated prevention (AMP) trials. In both cases, we find a 2:1:1 dose emphasizing VRC07 is nearly optimal. Our approach can be immediately applied to optimize the next generation of combination antibody prevention and cure studies.

With the great advancements in experimental data, computational power and learning algorithms, artificial intelligence (AI) based drug design has begun to gain momentum recently. AI-based drug design has great promise to revolutionize pharmaceutical industries by significantly reducing the time and cost in drug discovery processes. However, a major issue remains for all AI-based learning model that is efficient molecular representations. Here we propose Dowker complex (DC) based molecular interaction representations and Riemann Zeta function based molecular featurization, for the first time. Molecular interactions between proteins and ligands (or others) are modeled as Dowker complexes. A multiscale representation is generated by using a filtration process, during which a series of DCs are generated at different scales. Combinatorial (Hodge) Laplacian matrices are constructed from these DCs, and the Riemann zeta functions from their spectral information can be used as molecular descriptors. To validate our models, we consider protein-ligand binding affinity prediction. Our DC-based machine learning (DCML) models, in particular, DC-based gradient boosting tree (DC-GBT), are tested on three most-commonly used datasets, i.e., including PDBbind-2007, PDBbind-2013 and PDBbind-2016, and extensively compared with other existing state-of-the-art models. It has been found that our DC-based descriptors can achieve the state-of-the-art results and have better performance than all machine learning models with traditional molecular descriptors. Our Dowker complex based machine learning models can be used in other tasks in AI-based drug design and molecular data analysis.

Epigenetic cell memory allows distinct gene expression patterns to persist in different cell types despite a common genotype. Although different patterns can be maintained by the concerted action of transcription factors (TFs), it was proposed that long-term persistence hinges on chromatin state. Here, we study how the dynamics of chromatin state affect memory, and focus on a biologically motivated circuit motif, among histones and DNA modifications, that mediates the action of TFs on gene expression. Memory arises from time-scale separation among three circuit's constituent processes: basal erasure, auto and cross-catalysis, and recruited erasure of modifications. When the two latter processes are sufficiently faster than the former, the circuit exhibits bistability and hysteresis, allowing active and repressed gene states to coexist and persist after TF stimulus removal. The duration of memory is stochastic with a mean value that increases as time-scale separation increases, but more so for the repressed state. This asymmetry stems from the cross-catalysis between repressive histone modifications and DNA methylation and is enhanced by the relatively slower decay rate of the latter. Nevertheless, TF-mediated positive autoregulation can rebalance this asymmetry and even confers robustness of active states to repressive stimuli. More generally, by wiring positively autoregulated chromatin modification circuits under time scale separation, long-term distinct gene expression patterns arise, which are also robust to failure in the regulatory links.

Residue-level coarse-grained (CG) models have become one of the most popular tools in biomolecular simulations in the trade-off between modeling accuracy and computational efficiency. To investigate large-scale biological phenomena in molecular dynamics (MD) simulations with CG models, unified treatments of proteins and nucleic acids, as well as efficient parallel computations, are indispensable. In the GENESIS MD software, we implement several residue-level CG models, covering structure-based and context-based potentials for both well-folded biomolecules and intrinsically disordered regions. An amino acid residue in protein is represented as a single CG particle centered at the Cα atom position, while a nucleotide in RNA or DNA is modeled with three beads. Then, a single CG particle represents around ten heavy atoms in both proteins and nucleic acids. The input data in CG MD simulations are treated as GROMACS-style input files generated from a newly developed toolbox, GENESIS-CG-tool. To optimize the performance in CG MD simulations, we utilize multiple neighbor lists, each of which is attached to a different nonbonded interaction potential in the cell-linked list method. We found that random number generations for Gaussian distributions in the Langevin thermostat are one of the bottlenecks in CG MD simulations. Therefore, we parallelize the computations with message-passing-interface (MPI) to improve the performance on PC clusters or supercomputers. We simulate Herpes simplex virus (HSV) type 2 B-capsid and chromatin models containing more than 1,000 nucleosomes in GENESIS as examples of large-scale biomolecular simulations with residue-level CG models. This framework extends accessible spatial and temporal scales by multi-scale simulations to study biologically relevant phenomena, such as genome-scale chromatin folding or phase-separated membrane-less condensations.

Microfluidic capacities for both recreating and monitoring cell cultures have opened the door to the use of Data Science and Machine Learning tools for understanding and simulating tumor evolution under controlled conditions. In this work, we show how these techniques could be applied to study Glioblastoma, the deadliest and most frequent primary brain tumor. In particular, we study Glioblastoma invasion using the recent concept of Physically-Guided Neural Networks with Internal Variables (PGNNIV), able to combine data obtained from microfluidic devices and some physical knowledge governing the tumor evolution. The physics is introduced in the network structure by means of a nonlinear advection-diffusion-reaction partial differential equation that models the Glioblastoma evolution. On the other hand, multilayer perceptrons combined with a nodal deconvolution technique are used for learning the go or grow metabolic behavior which characterises the Glioblastoma invasion. The PGNNIV is here trained using synthetic data obtained from in silico tests created under different oxygenation conditions, using a previously validated model. The unravelling capacity of PGNNIV enables discovering complex metabolic processes in a non-parametric way, thus giving explanatory capacity to the networks, and, as a consequence, surpassing the predictive power of any parametric approach and for any kind of stimulus. Besides, the possibility of working, for a particular tumor, with different boundary and initial conditions, permits the use of PGNNIV for defining virtual therapies and for drug design, thus making the first steps towards in silico personalised medicine.

We present a spatial model describing the growth of a photosynthetic microalgae biofilm. In this 2D-model we consider photosynthesis, cell carbon accumulation, extracellular matrix excretion, and mortality. The rate of each of these mechanisms is given by kinetic laws regulated by light, nitrate, oxygen and inorganic carbon. The model is based on mixture theory and the behaviour of each component is defined on one hand by mass conservation, which takes into account biological features of the system, and on the other hand by conservation of momentum, which expresses the physical properties of the components. The model simulates the biofilm structural dynamics following an initial colonization phase. It shows that a 75 μm thick active region drives the biofilm development. We then determine the optimal harvesting period and biofilm height which maximize productivity. Finally, different harvesting patterns are tested and their effect on biofilm structure are discussed. The optimal strategy differs whether the objective is to recover the total biofilm or just the algal biomass.

Arousal levels strongly affect task performance. Yet, what arousal level is optimal for a task depends on its difficulty. Easy task performance peaks at higher arousal levels, whereas performance on difficult tasks displays an inverted U-shape relationship with arousal, peaking at medium arousal levels, an observation first made by Yerkes and Dodson in 1908. It is commonly proposed that the noradrenergic locus coeruleus system regulates these effects on performance through a widespread release of noradrenaline resulting in changes of cortical gain. This account, however, does not explain why performance decays with high arousal levels only in difficult, but not in simple tasks. Here, we present a mechanistic model that revisits the Yerkes-Dodson effect from a sensory perspective: a deep convolutional neural network augmented with a global gain mechanism reproduced the same interaction between arousal state and task difficulty in its performance. Investigating this model revealed that global gain states differentially modulated sensory information encoding across the processing hierarchy, which explained their differential effects on performance on simple versus difficult tasks. These findings offer a novel hierarchical sensory processing account of how, and why, arousal state affects task performance.

The drivers behind regional differences of SARS-CoV-2 spread on finer spatio-temporal scales are yet to be fully understood. Here we develop a data-driven modelling approach based on an age-structured compartmental model that compares 116 Austrian regions to a suitably chosen control set of regions to explain variations in local transmission rates through a combination of meteorological factors, non-pharmaceutical interventions and mobility. We find that more than 60% of the observed regional variations can be explained by these factors. Decreasing temperature and humidity, increasing cloudiness, precipitation and the absence of mitigation measures for public events are the strongest drivers for increased virus transmission, leading in combination to a doubling of the transmission rates compared to regions with more favourable weather. We conjecture that regions with little mitigation measures for large events that experience shifts toward unfavourable weather conditions are particularly predisposed as nucleation points for the next seasonal SARS-CoV-2 waves.

For cortical motor activity, the relationships between different body part representations is unknown. Through reciprocal body part relationships, functionality of cortical motor areas with respect to whole body motor control can be characterized. In the current study, we investigate the relationship between body part representations within individual neuronal populations in motor cortices, following a 7 Tesla fMRI 18-body-part motor experiment in combination with our newly developed non-rigid population Response Field (pRF) model and graph theory. The non-rigid pRF metrics reveal somatotopic structures in all included motor cortices covering frontal, parietal, medial and insular cortices and that neuronal populations in primary sensorimotor cortex respond to fewer body parts than secondary motor cortices. Reciprocal body part relationships are estimated in terms of uniqueness, clique-formation, and influence. We report unique response profiles for the knee, a clique of body parts surrounding the ring finger, and a central role for the shoulder and wrist. These results reveal associations among body parts from the perspective of the central nervous system, while being in agreement with intuitive notions of body part usage.

Ratiometric time-lapse FRET analysis requires a robust and accurate processing pipeline to eliminate bias in intensity measurements on fluorescent images before further quantitative analysis can be conducted. This level of robustness can only be achieved by supplementing automated tools with built-in flexibility for manual ad-hoc adjustments. FRET-IBRA is a modular and fully parallelized configuration file-based tool written in Python. It simplifies the FRET processing pipeline to achieve accurate, registered, and unified ratio image stacks. The flexibility of this tool to handle discontinuous image frame sequences with tailored configuration parameters further streamlines the processing of outliers and time-varying effects in the original microscopy images. FRET-IBRA offers cluster-based channel background subtraction, photobleaching correction, and ratio image construction in an all-in-one solution without the need for multiple applications, image format conversions, and/or plug-ins. The package accepts a variety of input formats and outputs TIFF image stacks along with performance measures to detect both the quality and failure of the background subtraction algorithm on a per frame basis. Furthermore, FRET-IBRA outputs images with superior signal-to-noise ratio and accuracy in comparison to existing background subtraction solutions, whilst maintaining a fast runtime. We have used the FRET-IBRA package extensively to quantify the spatial distribution of calcium ions during pollen tube growth under mechanical constraints. Benchmarks against existing tools clearly demonstrate the need for FRET-IBRA in extracting reliable insights from FRET microscopy images of dynamic physiological processes at high spatial and temporal resolution. The source code for Linux and Mac operating systems is released under the BSD license and, along with installation instructions, test images, example configuration files, and a step-by-step tutorial, is freely available at github.com/gmunglani/fret-ibra.

Connecting protein sequence to function is becoming increasingly relevant since high-throughput sequencing studies accumulate large amounts of genomic data. In order to go beyond the existing database annotation, it is fundamental to understand the mechanisms underlying functional inheritance and divergence. If the homology relationship between proteins is known, can we determine whether the function diverged? In this work, we analyze different possibilities of protein sequence evolution after gene duplication and identify "inter-paralog inversions", i.e., sites where the relationship between the ancestry and the functional signal is decoupled. The amino acids in these sites are masked from being recognized by other prediction tools. Still, they play a role in functional divergence and could indicate a shift in protein function. We develop a method to specifically recognize inter-paralog amino acid inversions in a phylogeny and test it on real and simulated datasets. In a dataset built from the Epidermal Growth Factor Receptor (EGFR) sequences found in 88 fish species, we identify 19 amino acid sites that went through inversion after gene duplication, mostly located at the ligand-binding extracellular domain. Our work uncovers an outcome of protein duplications with direct implications in protein functional annotation and sequence evolution. The developed method is optimized to work with large protein datasets and can be readily included in a targeted protein analysis pipeline.

For de novo mutational signature analysis, the critical first step is to decide how many signatures should be expected in a cancer genomics study. An incorrect number could mislead downstream analyses. Here we present SUITOR (Selecting the nUmber of mutatIonal signaTures thrOugh cRoss-validation), an unsupervised cross-validation method that requires little assumptions and no numerical approximations to select the optimal number of signatures without overfitting the data. In vitro studies and in silico simulations demonstrated that SUITOR can correctly identify signatures, some of which were missed by other widely used methods. Applied to 2,540 whole-genome sequenced tumors across 22 cancer types, SUITOR selected signatures with the smallest prediction errors and almost all signatures of breast cancer selected by SUITOR were validated in an independent breast cancer study. SUITOR is a powerful tool to select the optimal number of mutational signatures, facilitating downstream analyses with etiological or therapeutic importance.

Update of bioRxiv. 2021 Sep 23;:

Variant allele frequencies (VAF) encode ongoing evolution and subclonal selection in growing tumours. However, existing methods that utilize VAF information for cancer evolutionary inference are compressive, slow, or incorrectly specify the underlying cancer evolutionary dynamics. Here, we provide a proof-of-principle synthetic supervised learning method, TumE, that integrates simulated models of cancer evolution with Bayesian neural networks, to infer ongoing selection in bulk-sequenced single tumour biopsies. Analyses in synthetic and patient tumours show that TumE significantly improves both accuracy and inference time per sample when detecting positive selection, deconvoluting selected subclonal populations, and estimating subclone frequency. Importantly, we show how transfer learning can leverage stored knowledge within TumE models for related evolutionary inference tasks-substantially reducing data and computational time for further model development and providing a library of recyclable deep learning models for the cancer evolution community. This extensible framework provides a foundation and future directions for harnessing progressive computational methods for the benefit of cancer genomics and, in turn, the cancer patient.

Cortical circuits process information by rich recurrent interactions between excitatory neurons and inhibitory interneurons. One of the prime functions of interneurons is to stabilize the circuit by feedback inhibition, but the level of specificity on which inhibitory feedback operates is not fully resolved. We hypothesized that inhibitory circuits could enable separate feedback control loops for different synaptic input streams, by means of specific feedback inhibition to different neuronal compartments. To investigate this hypothesis, we adopted an optimization approach. Leveraging recent advances in training spiking network models, we optimized the connectivity and short-term plasticity of interneuron circuits for compartment-specific feedback inhibition onto pyramidal neurons. Over the course of the optimization, the interneurons diversified into two classes that resembled parvalbumin (PV) and somatostatin (SST) expressing interneurons. Using simulations and mathematical analyses, we show that the resulting circuit can be understood as a neural decoder that inverts the nonlinear biophysical computations performed within the pyramidal cells. Our model provides a proof of concept for studying structure-function relations in cortical circuits by a combination of gradient-based optimization and biologically plausible phenomenological models.

Myocardial ischemia, injury and infarction (MI) are the three stages of acute coronary syndrome (ACS). In the past two decades, a great number of studies focused on myocardial ischemia and MI individually, and showed that the occurrence of reentrant arrhythmias is often associated with myocardial ischemia or MI. However, arrhythmogenic mechanisms in the tissue with various degrees of remodeling in the ischemic heart have not been fully understood. In this study, biophysical detailed single-cell models of ischemia 1a, 1b, and MI were developed to mimic the electrophysiological remodeling at different stages of ACS. 2D tissue models with different distributions of ischemia and MI areas were constructed to investigate the mechanisms of the initiation of reentrant waves during the progression of ischemia. Simulation results in 2D tissues showed that the vulnerable windows (VWs) in simultaneous presence of multiple ischemic conditions were associated with the dynamics of wave propagation in the tissues with each single pathological condition. In the tissue with multiple pathological conditions, reentrant waves were mainly induced by two different mechanisms: one is the heterogeneity along the excitation wavefront, especially the abrupt variation in conduction velocity (CV) across the border of ischemia 1b and MI, and the other is the decreased safe factor (SF) for conduction at the edge of the tissue in MI region which is attributed to the increased excitation threshold of MI region. Finally, the reentrant wave was observed in a 3D model with a scar reconstructed from MRI images of a MI patient. These comprehensive findings provide novel insights for understanding the arrhythmic risk during the progression of myocardial ischemia and highlight the importance of the multiple pathological stages in designing medical therapies for arrhythmias in ischemia.

Collective, coordinated cellular motions underpin key processes in all multicellular organisms, yet it has been difficult to simultaneously express the 'rules' behind these motions in clear, interpretable forms that effectively capture high-dimensional cell-cell interaction dynamics in a manner that is intuitive to the researcher. Here we apply deep attention networks to analyze several canonical living tissues systems and present the underlying collective migration rules for each tissue type using only cell migration trajectory data. We use these networks to learn the behaviors of key tissue types with distinct collective behaviors-epithelial, endothelial, and metastatic breast cancer cells-and show how the results complement traditional biophysical approaches. In particular, we present attention maps indicating the relative influence of neighboring cells to the learned turning decisions of a 'focal cell'-the primary cell of interest in a collective setting. Colloquially, we refer to this learned relative influence as 'attention', as it serves as a proxy for the physical parameters modifying the focal cell's future motion as a function of each neighbor cell. These attention networks reveal distinct patterns of influence and attention unique to each model tissue. Endothelial cells exhibit tightly focused attention on their immediate forward-most neighbors, while cells in more expansile epithelial tissues are more broadly influenced by neighbors in a relatively large forward sector. Attention maps of ensembles of more mesenchymal, metastatic cells reveal completely symmetric attention patterns, indicating the lack of any particular coordination or direction of interest. Moreover, we show how attention networks are capable of detecting and learning how these rules change based on biophysical context, such as location within the tissue and cellular crowding. That these results require only cellular trajectories and no modeling assumptions highlights the potential of attention networks for providing further biological insights into complex cellular systems.

K-Ras activating mutations are significantly associated with tumor progression and aggressive metastatic behavior in various human cancers including pancreatic cancer. So far, despite a large number of concerted efforts, targeting of mutant-type K-Ras has not been successful. In this regard, we aimed to target this oncogene by a combinational approach consisting of small peptide and small molecule inhibitors. Based on a comprehensive analysis of structural and physicochemical properties of predominantly K-Ras mutants, an anti-cancer peptide library and a small molecule library were screened to simultaneously target oncogenic mutations and functional domains of mutant-type K-Ras located in the P-loop, switch I, and switch II regions. The selected peptide and small molecule showed notable binding affinities to their corresponding binding sites, and hindered the growth of tumor cells carrying K-RasG12D and K-RasG12C mutations. Of note, the expression of K-Ras downstream genes (i.e., CTNNB1, CCND1) was diminished in the treated Kras-positive cells. In conclusion, our combinational platform signifies a new potential for blockade of oncogenic K-Ras and thereby prevention of tumor progression and metastasis. However, further validations are still required regarding the in vitro and in vivo efficacy and safety of this approach.

Studying isoform expression at the microscopic level has always been a challenging task. A classical example is kidney, where glomerular and tubulo-interstitial compartments carry out drastically different physiological functions and thus presumably their isoform expression also differs. We aim at developing an experimental and computational pipeline for identifying isoforms at microscopic structure-level. We microdissected glomerular and tubulo-interstitial compartments from healthy human kidney tissues from two cohorts. The two compartments were separately sequenced with the PacBio RS II platform. These transcripts were then validated using transcripts of the same samples by the traditional Illumina RNA-Seq protocol, distinct Illumina RNA-Seq short reads from European Renal cDNA Bank (ERCB) samples, and annotated GENCODE transcript list, thus identifying novel transcripts. We identified 14,739 and 14,259 annotated transcripts, and 17,268 and 13,118 potentially novel transcripts in the glomerular and tubulo-interstitial compartments, respectively. Of note, relying solely on either short or long reads would have resulted in many erroneous identifications. We identified distinct pathways involved in glomerular and tubulo-interstitial compartments at the isoform level, creating an important experimental and computational resource for the kidney research community.

B cells undergo rapid mutation and selection for antibody binding affinity when producing antibodies capable of neutralizing pathogens. This evolutionary process can be intermixed with migration between tissues, differentiation between cellular subsets, and switching between functional isotypes. B cell receptor (BCR) sequence data has the potential to elucidate important information about these processes. However, there is currently no robust, generalizable framework for making such inferences from BCR sequence data. To address this, we develop three parsimony-based summary statistics to characterize migration, differentiation, and isotype switching along B cell phylogenetic trees. We use simulations to demonstrate the effectiveness of this approach. We then use this framework to infer patterns of cellular differentiation and isotype switching from high throughput BCR sequence datasets obtained from patients in a study of HIV infection and a study of food allergy. These methods are implemented in the R package dowser, available at https://dowser.readthedocs.io.

In HIV-1-infected individuals, transmitted/founder (TF) virus contributes to establish new infection and expands during the acute phase of infection, while chronic control (CC) virus emerges during the chronic phase of infection. TF viruses are more resistant to interferon-alpha (IFN-α)-mediated antiviral effects than CC virus, however, its virological relevance in infected individuals remains unclear. Here we perform an experimental-mathematical investigation and reveal that IFN-α strongly inhibits cell-to-cell infection by CC virus but only weakly affects that by TF virus. Surprisingly, IFN-α enhances cell-free infection of HIV-1, particularly that of CC virus, in a virus-cell density-dependent manner. We further demonstrate that LY6E, an IFN-stimulated gene, can contribute to the density-dependent enhancement of cell-free HIV-1 infection. Altogether, our findings suggest that the major difference between TF and CC viruses can be explained by their resistance to IFN-α-mediated inhibition of cell-to-cell infection and their sensitivity to IFN-α-mediated enhancement of cell-free infection.

Eukaryotic genomes often encode multiple transporters for the same nutrient. For example, budding yeast has 17 hexose transporters (HXTs), all of which potentially transport glucose. Using mathematical modelling, we show that transporters that use either facilitated diffusion or symport can have a rate-affinity tradeoff, where an increase in the maximal rate of transport decreases the transporter's apparent affinity. These changes affect the import flux non-monotonically, and for a given concentration of extracellular nutrient there is one transporter, characterised by its affinity, that has a higher import flux than any other. Through encoding multiple transporters, cells can therefore mitigate the tradeoff by expressing those transporters with higher affinities in lower concentrations of nutrients. We verify our predictions using fluorescent tagging of seven HXT genes in budding yeast and follow their expression over time in batch culture. Using the known affinities of the corresponding transporters, we show that their regulation in glucose is broadly consistent with a rate-affinity tradeoff: as glucose falls, the levels of the different transporters peak in an order that mostly follows their affinity for glucose. More generally, evolution is constrained by tradeoffs. Our findings indicate that one such tradeoff often occurs in the cellular transport of nutrients.

Dendritic spines are highly dynamic neuronal compartments that control the synaptic transmission between neurons. Spines form ultrastructural units, coupling synaptic contact sites to the dendritic shaft and often harbor a spine apparatus organelle, composed of smooth endoplasmic reticulum, which is responsible for calcium sequestration and release into the spine head and neck. The spine apparatus has recently been linked to synaptic plasticity in adult human cortical neurons. While the morphological heterogeneity of spines and their intracellular organization has been extensively demonstrated in animal models, the influence of spine apparatus organelles on critical signaling pathways, such as calcium-mediated dynamics, is less well known in human dendritic spines. In this study we used serial transmission electron microscopy to anatomically reconstruct nine human cortical spines in detail as a basis for modeling and simulation of the calcium dynamics between spine and dendrite. The anatomical study of reconstructed human dendritic spines revealed that the size of the postsynaptic density correlates with spine head volume and that the spine apparatus volume is proportional to the spine volume. Using a newly developed simulation pipeline, we have linked these findings to spine-to-dendrite calcium communication. While the absence of a spine apparatus, or the presence of a purely passive spine apparatus did not enable any of the reconstructed spines to relay a calcium signal to the dendritic shaft, the calcium-induced calcium release from this intracellular organelle allowed for finely tuned "all-or-nothing" spine-to-dendrite calcium coupling; controlled by spine morphology, neck plasticity, and ryanodine receptors. Our results suggest that spine apparatus organelles are strategically positioned in the neck of human dendritic spines and demonstrate their potential relevance to the maintenance and regulation of spine-to-dendrite calcium communication.

Natural compounds constitute a rich resource of potential small molecule therapeutics. While experimental access to this resource is limited due to its vast diversity and difficulties in systematic purification, computational assessment of structural similarity with known therapeutic molecules offers a scalable approach. Here, we assessed functional similarity between natural compounds and approved drugs by combining multiple chemical similarity metrics and physicochemical properties using a machine-learning approach. We computed pairwise similarities between 1410 drugs for training classification models and used the drugs shared protein targets as class labels. The best performing models were random forest which gave an average area under the ROC of 0.9, Matthews correlation coefficient of 0.35, and F1 score of 0.33, suggesting that it captured the structure-activity relation well. The models were then used to predict protein targets of circa 11k natural compounds by comparing them with the drugs. This revealed therapeutic potential of several natural compounds, including those with support from previously published sources as well as those hitherto unexplored. We experimentally validated one of the predicted pair's activities, viz., Cox-1 inhibition by 5-methoxysalicylic acid, a molecule commonly found in tea, herbs and spices. In contrast, another natural compound, 4-isopropylbenzoic acid, with the highest similarity score when considering most weighted similarity metric but not picked by our models, did not inhibit Cox-1. Our results demonstrate the utility of a machine-learning approach combining multiple chemical features for uncovering protein binding potential of natural compounds.

Tumors often harbor orders of magnitude more mutations than healthy tissues. The increased number of mutations may be due to an elevated mutation rate or frequent cell death and correspondingly rapid cell turnover, or a combination of the two. It is difficult to disentangle these two mechanisms based on widely available bulk sequencing data, where sequences from individual cells are intermixed and, thus, the cell lineage tree of the tumor cannot be resolved. Here we present a method that can simultaneously estimate the cell turnover rate and the rate of mutations from bulk sequencing data. Our method works by simulating tumor growth and finding the parameters with which the observed data can be reproduced with maximum likelihood. Applying this method to a real tumor sample, we find that both the mutation rate and the frequency of death may be high.

The transformation of synaptic input into action potential output is a fundamental single-cell computation resulting from the complex interaction of distinct cellular morphology and the unique expression profile of ion channels that define the cellular phenotype. Experimental studies aimed at uncovering the mechanisms of the transfer function have led to important insights, yet are limited in scope by technical feasibility, making biophysical simulations an attractive complementary approach to push the boundaries in our understanding of cellular computation. Here we take a data-driven approach by utilizing high-resolution morphological reconstructions and patch-clamp electrophysiology data together with a multi-objective optimization algorithm to build two populations of biophysically detailed models of murine hippocampal CA3 pyramidal neurons based on the two principal cell types that comprise this region. We evaluated the performance of these models and find that our approach quantitatively matches the cell type-specific firing phenotypes and recapitulate the intrinsic population-level variability in the data. Moreover, we confirm that the conductance values found by the optimization algorithm are consistent with differentially expressed ion channel genes in single-cell transcriptomic data for the two cell types. We then use these models to investigate the cell type-specific biophysical properties involved in the generation of complex-spiking output driven by synaptic input through an information-theoretic treatment of their respective transfer functions. Our simulations identify a host of cell type-specific biophysical mechanisms that define the morpho-functional phenotype to shape the cellular transfer function and place these findings in the context of a role for bursting in CA3 recurrent network synchronization dynamics.

We investigated the effects of updating age-specific social contact matrices to match evolving demography on vaccine impact estimates. We used a dynamic transmission model of tuberculosis in India as a case study. We modelled four incremental methods to update contact matrices over time, where each method incorporated its predecessor: fixed contact matrix (M0), preserved contact reciprocity (M1), preserved contact assortativity (M2), and preserved average contacts per individual (M3). We updated the contact matrices of a deterministic compartmental model of tuberculosis transmission, calibrated to epidemiologic data between 2000 and 2019 derived from India. We additionally calibrated the M0, M2, and M3 models to the 2050 TB incidence rate projected by the calibrated M1 model. We stratified age into three groups, children (<15y), adults (≥15y, <65y), and the elderly (≥65y), using World Population Prospects demographic data, between which we applied POLYMOD-derived social contact matrices. We simulated an M72-AS01E-like tuberculosis vaccine delivered from 2027 and estimated the per cent TB incidence rate reduction (IRR) in 2050 under each update method. We found that vaccine impact estimates in all age groups remained relatively stable between the M0-M3 models, irrespective of vaccine-targeting by age group. The maximum difference in impact, observed following adult-targeted vaccination, was 7% in the elderly, in whom we observed IRRs of 19% (uncertainty range 13-32), 20% (UR 13-31), 22% (UR 14-37), and 26% (UR 18-38) following M0, M1, M2 and M3 updates, respectively. We found that model-based TB vaccine impact estimates were relatively insensitive to demography-matched contact matrix updates in an India-like demographic and epidemiologic scenario. Current model-based TB vaccine impact estimates may be reasonably robust to the lack of contact matrix updates, but further research is needed to confirm and generalise this finding.

The coactivator KIX of CBP uses two binding surfaces to recognize multiple activators and exhibits allostery in ternary complex formation. Activator•coactivator interactions are central to transcriptional regulation, yet the microscopic origins of allostery in dynamic proteins like KIX are largely unknown. Here, we investigate the molecular recognition and allosteric manifestations involved in two KIX ternary systems c-Myb•KIX•MLL and pKID•KIX•MLL. Exploring the hypothesis that binary complex formation prepays an entropic cost for positive cooperativity, we utilize molecular dynamics simulations, side chain methyl order parameters, and differential scanning fluorimetry (DSF) to explore conformational entropy changes in KIX. The protein's configurational micro-states from structural clustering highlight the utility of protein plasticity in molecular recognition and allostery. We find that apo KIX occupies a wide distribution of lowly-populated configurational states. Each binding partner has its own suite of KIX states that it selects, building a model of molecular recognition fingerprints. Allostery is maximized with MLL pre-binding, which corresponds to the observation of a significant reduction in KIX micro-states observed when MLL binds. With all binding partners, the changes in KIX conformational entropy arise predominantly from changes in the most flexible loop. Likewise, we find that a small molecule and mutations allosterically inhibit/enhance activator binding by tuning loop dynamics, suggesting that loop-targeting chemical probes could be developed to alter KIX•activator interactions. Experimentally capturing KIX stabilization is challenging, particularly because of the disordered nature of particular activators. However, DSF melting curves allow for inference of relative entropic changes that occur across complexes, which we compare to our computed entropy changes using simulation methyl order parameters.

Machine learning-based classification approaches are widely used to predict host phenotypes from microbiome data. Classifiers are typically employed by considering operational taxonomic units or relative abundance profiles as input features. Such types of data are intrinsically sparse, which opens the opportunity to make predictions from the presence/absence rather than the relative abundance of microbial taxa. This also poses the question whether it is the presence rather than the abundance of particular taxa to be relevant for discrimination purposes, an aspect that has been so far overlooked in the literature. In this paper, we aim at filling this gap by performing a meta-analysis on 4,128 publicly available metagenomes associated with multiple case-control studies. At species-level taxonomic resolution, we show that it is the presence rather than the relative abundance of specific microbial taxa to be important when building classification models. Such findings are robust to the choice of the classifier and confirmed by statistical tests applied to identifying differentially abundant/present taxa. Results are further confirmed at coarser taxonomic resolutions and validated on 4,026 additional 16S rRNA samples coming from 30 public case-control studies.

The design of vaccines against highly mutable pathogens, such as HIV and influenza, requires a detailed understanding of how the adaptive immune system responds to encountering multiple variant antigens (Ags). Here, we describe a multiscale model of B cell receptor (BCR) affinity maturation that employs actual BCR nucleotide sequences and treats BCR/Ag interactions in atomistic detail. We apply the model to simulate the maturation of a broadly neutralizing Ab (bnAb) against HIV. Starting from a germline precursor sequence of the VRC01 anti-HIV Ab, we simulate BCR evolution in response to different vaccination protocols and different Ags, which were previously designed by us. The simulation results provide qualitative guidelines for future vaccine design and reveal unique insights into bnAb evolution against the CD4 binding site of HIV. Our model makes possible direct comparisons of simulated BCR populations with results of deep sequencing data, which will be explored in future applications.

Erratum in PLoS Comput Biol. 2022 Nov 15;18(11):e1010703.

Biochemical interactions in systems and synthetic biology are often modeled with chemical reaction networks (CRNs). CRNs provide a principled modeling environment capable of expressing a huge range of biochemical processes. In this paper, we present a software toolbox, written in Python, that compiles high-level design specifications represented using a modular library of biochemical parts, mechanisms, and contexts to CRN implementations. This compilation process offers four advantages. First, the building of the actual CRN representation is automatic and outputs Systems Biology Markup Language (SBML) models compatible with numerous simulators. Second, a library of modular biochemical components allows for different architectures and implementations of biochemical circuits to be represented succinctly with design choices propagated throughout the underlying CRN automatically. This prevents the often occurring mismatch between high-level designs and model dynamics. Third, high-level design specification can be embedded into diverse biomolecular environments, such as cell-free extracts and in vivo milieus. Finally, our software toolbox has a parameter database, which allows users to rapidly prototype large models using very few parameters which can be customized later. By using BioCRNpyler, users ranging from expert modelers to novice script-writers can easily build, manage, and explore sophisticated biochemical models using diverse biochemical implementations, environments, and modeling assumptions.

The development of high-throughput genomic technologies associated with recent genetic perturbation techniques such as short hairpin RNA (shRNA), gene trapping, or gene editing (CRISPR/Cas9) has made it possible to obtain large perturbation data sets. These data sets are invaluable sources of information regarding the function of genes, and they offer unique opportunities to reverse engineer gene regulatory networks in specific cell types. Modular response analysis (MRA) is a well-accepted mathematical modeling method that is precisely aimed at such network inference tasks, but its use has been limited to rather small biological systems so far. In this study, we show that MRA can be employed on large systems with almost 1,000 network components. In particular, we show that MRA performance surpasses general-purpose mutual information-based algorithms. Part of these competitive results was obtained by the application of a novel heuristic that pruned MRA-inferred interactions a posteriori. We also exploited a block structure in MRA linear algebra to parallelize large system resolutions.

Recent advances in experimental and computational protein structure determination have provided access to high-quality structures for most human proteins and mutants thereof. However, linking changes in structure in protein mutants to functional impact remains an active area of method development. If successful, such methods can ultimately assist physicians in taking appropriate treatment decisions. This work presents three artificial neural network (ANN)-based predictive models that classify four key functional parameters of KCNQ1 variants as normal or dysfunctional using PSSM-based evolutionary and/or biophysical descriptors. Recent advances in predicting protein structure and variant properties with artificial intelligence (AI) rely heavily on the availability of evolutionary features and thus fail to directly assess the biophysical underpinnings of a change in structure and/or function. The central goal of this work was to develop an ANN model based on structure and physiochemical properties of KCNQ1 potassium channels that performs comparably or better than algorithms using only on PSSM-based evolutionary features. These biophysical features highlight the structure-function relationships that govern protein stability, function, and regulation. The input sensitivity algorithm incorporates the roles of hydrophobicity, polarizability, and functional densities on key functional parameters of the KCNQ1 channel. Inclusion of the biophysical features outperforms exclusive use of PSSM-based evolutionary features in predicting activation voltage dependence and deactivation time. As AI is increasingly applied to problems in biology, biophysical understanding will be critical with respect to 'explainable AI', i.e., understanding the relation of sequence, structure, and function of proteins. Our model is available at www.kcnq1predict.org.

The dynamical and physiological basis of alpha band activity and 1/fβ noise in the EEG are the subject of continued speculation. Here we conjecture, on the basis of empirical data analysis, that both of these features may be economically accounted for through a single process if the resting EEG is conceived of being the sum of multiple stochastically perturbed alpha band damped linear oscillators with a distribution of dampings (relaxation rates). The modulation of alpha-band and 1/fβ noise activity by changes in damping is explored in eyes closed (EC) and eyes open (EO) resting state EEG. We aim to estimate the distribution of dampings by solving an inverse problem applied to EEG power spectra. The characteristics of the damping distribution are examined across subjects, sensors and recording condition (EC/EO). We find that there are robust changes in the damping distribution between EC and EO recording conditions across participants. The estimated damping distributions are found to be predominantly bimodal, with the number and position of the modes related to the sharpness of the alpha resonance and the scaling (β) of the power spectrum (1/fβ). The results suggest that there exists an intimate relationship between resting state alpha activity and 1/fβ noise with changes in both governed by changes to the damping of the underlying alpha oscillatory processes. In particular, alpha-blocking is observed to be the result of the most weakly damped distribution mode becoming more heavily damped. The results suggest a novel way of characterizing resting EEG power spectra and provides new insight into the central role that damped alpha-band activity may play in characterising the spatio-temporal features of resting state EEG.

Pertussis has resurfaced in the UK, with incidence levels not seen since the 1980s. While the fundamental causes of this resurgence remain the subject of much conjecture, the study of historical patterns of pathogen diffusion can be illuminating. Here, we examined time series of pertussis incidence in the boroughs of Greater London from 1982 to 2013 to document the spatial epidemiology of this bacterial infection and to identify the potential drivers of its percolation. The incidence of pertussis over this period is characterized by 3 distinct stages: a period exhibiting declining trends with 4-year inter-epidemic cycles from 1982 to 1994, followed by a deep trough until 2006 and the subsequent resurgence. We observed systematic temporal trends in the age distribution of cases and the fade-out profile of pertussis coincident with increasing national vaccine coverage from 1982 to 1990. To quantify the hierarchy of epidemic phases across the boroughs of London, we used the Hilbert transform. We report a consistent pattern of spatial organization from 1982 to the early 1990s, with some boroughs consistently leading epidemic waves and others routinely lagging. To determine the potential drivers of these geographic patterns, a comprehensive parallel database of borough-specific features was compiled, comprising of demographic, movement and socio-economic factors that were used in statistical analyses to predict epidemic phase relationships among boroughs. Specifically, we used a combination of a feed-forward neural network (FFNN), and SHapley Additive exPlanations (SHAP) values to quantify the contribution of each covariate to model predictions. Our analyses identified a number of predictors of a borough's historical epidemic phase, specifically the age composition of households, the number of agricultural and skilled manual workers, latitude, the population of public transport commuters and high-occupancy households. Univariate regression analysis of the 2012 epidemic identified the ratio of cumulative unvaccinated children to the total population and population of Pakistan-born population to have moderate positive and negative association, respectively, with the timing of epidemic. In addition to providing a comprehensive overview of contemporary pertussis transmission in a large metropolitan population, this study has identified the characteristics that determine the spatial spread of this bacterium across the boroughs of London.

Genome wide association studies (GWASs) for complex traits have implicated thousands of genetic loci. Most GWAS-nominated variants lie in noncoding regions, complicating the systematic translation of these findings into functional understanding. Here, we leverage convolutional neural networks to assist in this challenge. Our computational framework, peaBrain, models the transcriptional machinery of a tissue as a two-stage process: first, predicting the mean tissue specific abundance of all genes and second, incorporating the transcriptomic consequences of genotype variation to predict individual abundance on a subject-by-subject basis. We demonstrate that peaBrain accounts for the majority (>50%) of variance observed in mean transcript abundance across most tissues and outperforms regularized linear models in predicting the consequences of individual genotype variation. We highlight the validity of the peaBrain model by calculating non-coding impact scores that correlate with nucleotide evolutionary constraint that are also predictive of disease-associated variation and allele-specific transcription factor binding. We further show how these tissue-specific peaBrain scores can be leveraged to pinpoint functional tissues underlying complex traits, outperforming methods that depend on colocalization of eQTL and GWAS signals. We subsequently: (a) derive continuous dense embeddings of genes for downstream applications; (b) highlight the utility of the model in predicting transcriptomic impact of small molecules and shRNA (on par with in vitro experimental replication of external test sets); (c) explore how peaBrain can be used to model difficult-to-study processes (such as neural induction); and (d) identify putatively functional eQTLs that are missed by high-throughput experimental approaches.

Segmenting three-dimensional (3D) microscopy images is essential for understanding phenomena like morphogenesis, cell division, cellular growth, and genetic expression patterns. Recently, deep learning (DL) pipelines have been developed, which claim to provide high accuracy segmentation of cellular images and are increasingly considered as the state of the art for image segmentation problems. However, it remains difficult to define their relative performances as the concurrent diversity and lack of uniform evaluation strategies makes it difficult to know how their results compare. In this paper, we first made an inventory of the available DL methods for 3D cell segmentation. We next implemented and quantitatively compared a number of representative DL pipelines, alongside a highly efficient non-DL method named MARS. The DL methods were trained on a common dataset of 3D cellular confocal microscopy images. Their segmentation accuracies were also tested in the presence of different image artifacts. A specific method for segmentation quality evaluation was adopted, which isolates segmentation errors due to under- or oversegmentation. This is complemented with a 3D visualization strategy for interactive exploration of segmentation quality. Our analysis shows that the DL pipelines have different levels of accuracy. Two of them, which are end-to-end 3D and were originally designed for cell boundary detection, show high performance and offer clear advantages in terms of adaptability to new data.

The types of epileptiform activity occurring in the sclerotic hippocampus with highest incidence are interictal-like events (II) and periodic ictal spiking (PIS). These activities are classified according to their event rates, but it is still unclear if these rate differences are consequences of underlying physiological mechanisms. Identifying new and more specific information related to these two activities may bring insights to a better understanding about the epileptogenic process and new diagnosis. We applied Poincaré map analysis and Recurrence Quantification Analysis (RQA) onto 35 in vitro electrophysiological signals recorded from slices of 12 hippocampal tissues surgically resected from patients with pharmacoresistant temporal lobe epilepsy. These analyzes showed that the II activity is related to chaotic dynamics, whereas the PIS activity is related to deterministic periodic dynamics. Additionally, it indicates that their different rates are consequence of different endogenous dynamics. Finally, by using two computational models we were able to simulate the transition between II and PIS activities. The RQA was applied to different periods of these simulations to compare the recurrences between artificial and real signals, showing that different ranges of regularity-chaoticity can be directly associated with the generation of PIS and II activities.

Single-cell mass cytometry, also known as cytometry by time of flight (CyTOF) is a powerful high-throughput technology that allows analysis of up to 50 protein markers per cell for the quantification and classification of single cells. Traditional manual gating utilized to identify new cell populations has been inadequate, inefficient, unreliable, and difficult to use, and no algorithms to identify both calibration and new cell populations has been well established. A deep learning with graphic cluster (DGCyTOF) visualization is developed as a new integrated embedding visualization approach in identifying canonical and new cell types. The DGCyTOF combines deep-learning classification and hierarchical stable-clustering methods to sequentially build a tri-layer construct for known cell types and the identification of new cell types. First, deep classification learning is constructed to distinguish calibration cell populations from all cells by softmax classification assignment under a probability threshold, and graph embedding clustering is then used to identify new cell populations sequentially. In the middle of two-layer, cell labels are automatically adjusted between new and unknown cell populations via a feedback loop using an iteration calibration system to reduce the rate of error in the identification of cell types, and a 3-dimensional (3D) visualization platform is finally developed to display the cell clusters with all cell-population types annotated. Utilizing two benchmark CyTOF databases comprising up to 43 million cells, we compared accuracy and speed in the identification of cell types among DGCyTOF, DeepCyTOF, and other technologies including dimension reduction with clustering, including Principal Component Analysis (PCA), Factor Analysis (FA), Independent Component Analysis (ICA), Isometric Feature Mapping (Isomap), t-distributed Stochastic Neighbor Embedding (t-SNE), and Uniform Manifold Approximation and Projection (UMAP) with k-means clustering and Gaussian mixture clustering. We observed the DGCyTOF represents a robust complete learning system with high accuracy, speed and visualization by eight measurement criteria. The DGCyTOF displayed F-scores of 0.9921 for CyTOF1 and 0.9992 for CyTOF2 datasets, whereas those scores were only 0.507 and 0.529 for the t-SNE+k-means; 0.565 and 0.59, for UMAP+ k-means. Comparison of DGCyTOF with t-SNE and UMAP visualization in accuracy demonstrated its approximately 35% superiority in predicting cell types. In addition, observation of cell-population distribution was more intuitive in the 3D visualization in DGCyTOF than t-SNE and UMAP visualization. The DGCyTOF model can automatically assign known labels to single cells with high accuracy using deep-learning classification assembling with traditional graph-clustering and dimension-reduction strategies. Guided by a calibration system, the model seeks optimal accuracy balance among calibration cell populations and unknown cell types, yielding a complete and robust learning system that is highly accurate in the identification of cell populations compared to results using other methods in the analysis of single-cell CyTOF data. Application of the DGCyTOF method to identify cell populations could be extended to the analysis of single-cell RNASeq data and other omics data.

We find that epidemic resurgence, defined as an upswing in the effective reproduction number (R) of the contagion from subcritical to supercritical values, is fundamentally difficult to detect in real time. Inherent latencies in pathogen transmission, coupled with smaller and intrinsically noisier case incidence across periods of subcritical spread, mean that resurgence cannot be reliably detected without significant delays of the order of the generation time of the disease, even when case reporting is perfect. In contrast, epidemic suppression (where R falls from supercritical to subcritical values) may be ascertained 5-10 times faster due to the naturally larger incidence at which control actions are generally applied. We prove that these innate limits on detecting resurgence only worsen when spatial or demographic heterogeneities are incorporated. Consequently, we argue that resurgence is more effectively handled proactively, potentially at the expense of false alarms. Timely responses to recrudescent infections or emerging variants of concern are more likely to be possible when policy is informed by a greater quality and diversity of surveillance data than by further optimisation of the statistical models used to process routine outbreak data.

The 3-dimensional fold of an RNA molecule is largely determined by patterns of intramolecular hydrogen bonds between bases. Predicting the base pairing network from the sequence, also referred to as RNA secondary structure prediction or RNA folding, is a nondeterministic polynomial-time (NP)-complete computational problem. The structure of the molecule is strongly predictive of its functions and biochemical properties, and therefore the ability to accurately predict the structure is a crucial tool for biochemists. Many methods have been proposed to efficiently sample possible secondary structure patterns. Classic approaches employ dynamic programming, and recent studies have explored approaches inspired by evolutionary and machine learning algorithms. This work demonstrates leveraging quantum computing hardware to predict the secondary structure of RNA. A Hamiltonian written in the form of a Binary Quadratic Model (BQM) is derived to drive the system toward maximizing the number of consecutive base pairs while jointly maximizing the average length of the stems. A Quantum Annealer (QA) is compared to a Replica Exchange Monte Carlo (REMC) algorithm programmed with the same objective function, with the QA being shown to be highly competitive at rapidly identifying low energy solutions. The method proposed in this study was compared to three algorithms from literature and, despite its simplicity, was found to be competitive on a test set containing known structures with pseudoknots.

The pathogenesis of many inflammatory diseases is a coordinated process involving metabolic dysfunctions and immune response-usually modulated by the production of cytokines and associated inflammatory molecules. In this work, we seek to understand how genes involved in pathogenesis which are often not associated with the immune system in an obvious way communicate with the immune system. We have embedded a network of human protein-protein interactions (PPI) from the STRING database with 14,707 human genes using feature learning that captures high confidence edges. We have found that our predicted Association Scores derived from the features extracted from STRING's high confidence edges are useful for predicting novel connections between genes, thus enabling the construction of a full map of predicted associations for all possible pairs between 14,707 human genes. In particular, we analyzed the pattern of associations for 126 cytokines and found that the six patterns of cytokine interaction with human genes are consistent with their functional classifications. To define the disease-specific roles of cytokines we have collected gene sets for 11,944 diseases from DisGeNET. We used these gene sets to predict disease-specific gene associations with cytokines by calculating the normalized average Association Scores between disease-associated gene sets and the 126 cytokines; this creates a unique profile of inflammatory genes (both known and predicted) for each disease. We validated our predicted cytokine associations by comparing them to known associations for 171 diseases. The predicted cytokine profiles correlate (p-value<0.0003) with the known ones in 95 diseases. We further characterized the profiles of each disease by calculating an "Inflammation Score" that summarizes different modes of immune responses. Finally, by analyzing subnetworks formed between disease-specific pathogenesis genes, hormones, receptors, and cytokines, we identified the key genes responsible for interactions between pathogenesis and inflammatory responses. These genes and the corresponding cytokines used by different immune disorders suggest unique targets for drug discovery.

Scientific research is shedding light on the interaction of the gut microbiome with the human host and on its role in human health. Existing machine learning methods have shown great potential in discriminating healthy from diseased microbiome states. Most of them leverage shotgun metagenomic sequencing to extract gut microbial species-relative abundances or strain-level markers. Each of these gut microbial profiling modalities showed diagnostic potential when tested separately; however, no existing approach combines them in a single predictive framework. Here, we propose the Multimodal Variational Information Bottleneck (MVIB), a novel deep learning model capable of learning a joint representation of multiple heterogeneous data modalities. MVIB achieves competitive classification performance while being faster than existing methods. Additionally, MVIB offers interpretable results. Our model adopts an information theoretic interpretation of deep neural networks and computes a joint stochastic encoding of different input data modalities. We use MVIB to predict whether human hosts are affected by a certain disease by jointly analysing gut microbial species-relative abundances and strain-level markers. MVIB is evaluated on human gut metagenomic samples from 11 publicly available disease cohorts covering 6 different diseases. We achieve high performance (0.80 < ROC AUC < 0.95) on 5 cohorts and at least medium performance on the remaining ones. We adopt a saliency technique to interpret the output of MVIB and identify the most relevant microbial species and strain-level markers to the model's predictions. We also perform cross-study generalisation experiments, where we train and test MVIB on different cohorts of the same disease, and overall we achieve comparable results to the baseline approach, i.e. the Random Forest. Further, we evaluate our model by adding metabolomic data derived from mass spectrometry as a third input modality. Our method is scalable with respect to input data modalities and has an average training time of < 1.4 seconds. The source code and the datasets used in this work are publicly available.

Accurate measurements of metabolic fluxes in living cells are central to metabolism research and metabolic engineering. The gold standard method is model-based metabolic flux analysis (MFA), where fluxes are estimated indirectly from mass isotopomer data with the use of a mathematical model of the metabolic network. A critical step in MFA is model selection: choosing what compartments, metabolites, and reactions to include in the metabolic network model. Model selection is often done informally during the modelling process, based on the same data that is used for model fitting (estimation data). This can lead to either overly complex models (overfitting) or too simple ones (underfitting), in both cases resulting in poor flux estimates. Here, we propose a method for model selection based on independent validation data. We demonstrate in simulation studies that this method consistently chooses the correct model in a way that is independent on errors in measurement uncertainty. This independence is beneficial, since estimating the true magnitude of these errors can be difficult. In contrast, commonly used model selection methods based on the χ2-test choose different model structures depending on the believed measurement uncertainty; this can lead to errors in flux estimates, especially when the magnitude of the error is substantially off. We present a new approach for quantification of prediction uncertainty of mass isotopomer distributions in other labelling experiments, to check for problems with too much or too little novelty in the validation data. Finally, in an isotope tracing study on human mammary epithelial cells, the validation-based model selection method identified pyruvate carboxylase as a key model component. Our results argue that validation-based model selection should be an integral part of MFA model development.

Comparing SARS-CoV-2 infection-induced gene expression signatures to drug treatment-induced gene expression signatures is a promising bioinformatic tool to repurpose existing drugs against SARS-CoV-2. The general hypothesis of signature-based drug repurposing is that drugs with inverse similarity to a disease signature can reverse disease phenotype and thus be effective against it. However, in the case of viral infection diseases, like SARS-CoV-2, infected cells also activate adaptive, antiviral pathways, so that the relationship between effective drug and disease signature can be more ambiguous. To address this question, we analysed gene expression data from in vitro SARS-CoV-2 infected cell lines, and gene expression signatures of drugs showing anti-SARS-CoV-2 activity. Our extensive functional genomic analysis showed that both infection and treatment with in vitro effective drugs leads to activation of antiviral pathways like NFkB and JAK-STAT. Based on the similarity-and not inverse similarity-between drug and infection-induced gene expression signatures, we were able to predict the in vitro antiviral activity of drugs. We also identified SREBF1/2, key regulators of lipid metabolising enzymes, as the most activated transcription factors by several in vitro effective antiviral drugs. Using a fluorescently labeled cholesterol sensor, we showed that these drugs decrease the cholesterol levels of plasma-membrane. Supplementing drug-treated cells with cholesterol reversed the in vitro antiviral effect, suggesting the depleting plasma-membrane cholesterol plays a key role in virus inhibitory mechanism. Our results can help to more effectively repurpose approved drugs against SARS-CoV-2, and also highlights key mechanisms behind their antiviral effect.

Understanding of the pairing statistics in solutions populated by a large number of distinct solute species with mutual interactions is a challenging topic, relevant in modeling the complexity of real biological systems. Here we describe, both experimentally and theoretically, the formation of duplexes in a solution of random-sequence DNA (rsDNA) oligomers of length L = 8, 12, 20 nucleotides. rsDNA solutions are formed by 4L distinct molecular species, leading to a variety of pairing motifs that depend on sequence complementarity and range from strongly bound, fully paired defectless helices to weakly interacting mismatched duplexes. Experiments and theory coherently combine revealing a hybridization statistics characterized by a prevalence of partially defected duplexes, with a distribution of type and number of pairing errors that depends on temperature. We find that despite the enormous multitude of inter-strand interactions, defectless duplexes are formed, involving a fraction up to 15% of the rsDNA chains at the lowest temperatures. Experiments and theory are limited here to equilibrium conditions.

Protein-protein interactions (PPIs) are key drivers of cell function and evolution. While it is widely assumed that most permanent PPIs are important for cellular function, it remains unclear whether transient PPIs are equally important. Here, we estimate and compare dispensable content among transient PPIs and permanent PPIs in human. Starting with a human reference interactome mapped by experiments, we construct a human structural interactome by building three-dimensional structural models for PPIs, and then distinguish transient PPIs from permanent PPIs using several structural and biophysical properties. We map common mutations from healthy individuals and disease-causing mutations onto the structural interactome, and perform structure-based calculations of the probabilities for common mutations (assumed to be neutral) and disease mutations (assumed to be mildly deleterious) to disrupt transient PPIs and permanent PPIs. Using Bayes' theorem we estimate that a similarly small fraction (<~20%) of both transient and permanent PPIs are completely dispensable, i.e., effectively neutral upon disruption. Hence, transient and permanent interactions are subject to similarly strong selective constraints in the human interactome.

The spread of COVID-19 caused by the SARS-CoV-2 virus has become a worldwide problem with devastating consequences. Here, we implement a comprehensive contact tracing and network analysis to find an optimized quarantine protocol to dismantle the chain of transmission of coronavirus with minimal disruptions to society. We track billions of anonymized GPS human mobility datapoints to monitor the evolution of the contact network of disease transmission before and after mass quarantines. As a consequence of the lockdowns, people's mobility decreases by 53%, which results in a drastic disintegration of the transmission network by 90%. However, this disintegration did not halt the spreading of the disease. Our analysis indicates that superspreading k-core structures persist in the transmission network to prolong the pandemic. Once the k-cores are identified, an optimized strategy to break the chain of transmission is to quarantine a minimal number of 'weak links' with high betweenness centrality connecting the large k-cores.

The paper presents a statistical-mechanics model for the kinetic selection of viral RNA molecules by packaging signals during the nucleation stage of the assembly of small RNA viruses. The effects of the RNA secondary structure and folding geometry of the packaging signals on the assembly activation energy barrier are encoded by a pair of characteristics: the wrapping number and the maximum ladder distance. Kinetic selection is found to be optimal when assembly takes place under conditions of supersaturation and also when the concentration ratio of capsid protein and viral RNA concentrations equals the stoichiometric ratio of assembled viral particles. As a function of the height of the activation energy barrier, there is a form of order-disorder transition such that for sufficiently low activation energy barriers, kinetic selectivity is erased by entropic effects associated with the number of assembly pathways.

Advances in single-cell RNA sequencing (scRNA-seq) have led to successes in discovering novel cell types and understanding cellular heterogeneity among complex cell populations through cluster analysis. However, cluster analysis is not able to reveal continuous spectrum of states and underlying gene expression programs (GEPs) shared across cell types. We introduce scAAnet, an autoencoder for single-cell non-linear archetypal analysis, to identify GEPs and infer the relative activity of each GEP across cells. We use a count distribution-based loss term to account for the sparsity and overdispersion of the raw count data and add an archetypal constraint to the loss function of scAAnet. We first show that scAAnet outperforms existing methods for archetypal analysis across different metrics through simulations. We then demonstrate the ability of scAAnet to extract biologically meaningful GEPs using publicly available scRNA-seq datasets including a pancreatic islet dataset, a lung idiopathic pulmonary fibrosis dataset and a prefrontal cortex dataset.

Application of epicardial patches constructed from human-induced pluripotent stem cell- derived cardiomyocytes (hiPSC-CMs) has been proposed as a long-term therapy to treat scarred hearts post myocardial infarction (MI). Understanding electrical interaction between engineered heart tissue patches (EHT) and host myocardium represents a key step toward a successful patch engraftment. EHT retain different electrical properties with respect to the host heart tissue due to the hiPSC-CMs immature phenotype, which may lead to increased arrhythmia risk. We developed a modelling framework to examine the influence of patch design on electrical activation at the engraftment site. We performed an in silico investigation of different patch design approaches to restore pre-MI activation properties and evaluated the associated arrhythmic risk. We developed an in silico cardiac electrophysiology model of a transmural cross section of host myocardium. The model featured an infarct region, an epicardial patch spanning the infarct region and a bath region. The patch is modelled as a layer of hiPSC-CM, combined with a layer of conductive polymer (CP). Tissue and patch geometrical dimensions and conductivities were incorporated through 10 modifiable model parameters. We validated our model against 4 independent experimental studies and showed that it can qualitatively reproduce their findings. We performed a global sensitivity analysis (GSA) to isolate the most important parameters, showing that the stimulus propagation is mainly governed by the scar depth, radius and conductivity when the scar is not transmural, and by the EHT patch conductivity when the scar is transmural. We assessed the relevance of small animal studies to humans by comparing simulations of rat, rabbit and human myocardium. We found that stimulus propagation paths and GSA sensitivity indices are consistent across species. We explored which EHT design variables have the potential to restore physiological propagation. Simulations predict that increasing EHT conductivity from 0.28 to 1-1.1 S/m recovered physiological activation in rat, rabbit and human. Finally, we assessed arrhythmia risk related to increasing EHT conductivity and tested increasing the EHT Na+ channel density as an alternative strategy to match healthy activation. Our results revealed a greater arrhythmia risk linked to increased EHT conductivity compared to increased Na+ channel density. We demonstrated that our modeling framework could capture the interaction between host and EHT patches observed in in vitro experiments. We showed that large (patch and tissue dimensions) and small (cardiac myocyte electrophysiology) scale differences between small animals and humans do not alter EHT patch effect on infarcted tissue. Our model revealed that only when the scar is transmural do EHT properties impact activation times and isolated the EHT conductivity as the main parameter influencing propagation. We predicted that restoring physiological activation by tuning EHT conductivity is possible but may promote arrhythmic behavior. Finally, our model suggests that acting on hiPSC-CMs low action potential upstroke velocity and lack of IK1 may restore pre-MI activation while not promoting arrhythmia.

As the most widespread viral infection transmitted by the Aedes mosquitoes, dengue has been estimated to cause 51 million febrile disease cases globally each year. Although sustained vector control remains key to reducing the burden of dengue, current understanding of the key factors that explain the observed variation in the short- and long-term vector control effectiveness across different transmission settings remains limited. We used a detailed individual-based model to simulate dengue transmission with and without sustained vector control over a 30-year time frame, under different transmission scenarios. Vector control effectiveness was derived for different time windows within the 30-year intervention period. We then used the extreme gradient boosting algorithm to predict the effectiveness of vector control given the simulation parameters, and the resulting machine learning model was interpreted using Shapley Additive Explanations. According to our simulation outputs, dengue transmission would be nearly eliminated during the early stage of sustained and intensive vector control, but over time incidence would gradually bounce back to the pre-intervention level unless the intervention is implemented at a very high level of intensity. The time point at which intervention ceases to be effective is strongly influenced not only by the intensity of vector control, but also by the pre-intervention transmission intensity and the individual-level heterogeneity in biting risk. Moreover, the impact of many transmission model parameters on the intervention effectiveness is shown to be modified by the intensity of vector control, as well as to vary over time. Our study has identified some of the critical drivers for the difference in the time-varying effectiveness of sustained vector control across different dengue endemic settings, and the insights obtained will be useful to inform future model-based studies that seek to predict the impact of dengue vector control in their local contexts.

The large taxonomic variability of microbial community composition is a consequence of the combination of environmental variability, mediated through ecological interactions, and stochasticity. Most of the analysis aiming to infer the biological factors determining this difference in community structure start by quantifying how much communities are similar in their composition, trough beta-diversity metrics. The central role that these metrics play in microbial ecology does not parallel with a quantitative understanding of their relationships and statistical properties. In particular, we lack a framework that reproduces the empirical statistical properties of beta-diversity metrics. Here we take a macroecological approach and introduce a model to reproduce the statistical properties of community similarity. The model is based on the statistical properties of individual communities and on a single tunable parameter, the correlation of species' carrying capacities across communities, which sets the difference of two communities. The model reproduces quantitatively the empirical values of several commonly-used beta-diversity metrics, as well as the relationships between them. In particular, this modeling framework naturally reproduces the negative correlation between overlap and dissimilarity, which has been observed in both empirical and experimental communities and previously related to the existence of universal features of community dynamics. In this framework, such correlation naturally emerges due to the effect of random sampling.