In the natural world, stimulus-outcome associations are often ambiguous, and most associations are highly complex and situation-dependent. Learning to disambiguate these complex associations to identify which specific outcomes will occur in which situations is critical for survival. Pavlovian occasion setters are stimuli that determine whether other stimuli will result in a specific outcome. Occasion setting is a well-established phenomenon, but very little investigation has been conducted on how occasion setters are disambiguated when they themselves are ambiguous (i.e., when they do not consistently signal whether another stimulus will be reinforced). In two preregistered studies, we investigated the role of higher-order Pavlovian occasion setting in humans. We developed and tested the first computational model predicting direct associative learning, traditional occasion setting (i.e., 1st-order occasion setting), and 2nd-order occasion setting. This model operationalizes stimulus ambiguity as a mechanism to engage in higher-order Pavlovian learning. Both behavioral and computational modeling results suggest that 2nd-order occasion setting was learned, as evidenced by lack and presence of transfer of occasion setting properties when expected and the superior fit of our 2nd-order occasion setting model compared to the 1st-order occasion setting or direct associations models. These results provide a controlled investigation into highly complex associative learning and may ultimately lead to improvements in the treatment of Pavlovian-based mental health disorders (e.g., anxiety disorders, substance use).

Unraveling the network of interactions in ecological communities is a daunting task. Common methods to infer interspecific interactions from cross-sectional data are based on co-occurrence measures. For instance, interactions in the human microbiome are often inferred from correlations between the abundances of bacterial phylogenetic groups across subjects. We tested whether such correlation-based methods are indeed reliable for inferring interaction networks. For this purpose, we simulated bacterial communities by means of the generalized Lotka-Volterra model, with variation in model parameters representing variability among hosts. Our results show that correlations can be indicative for presence of bacterial interactions, but only when measurement noise is low relative to the variation in interaction strengths between hosts. Indication of interaction was affected by type of interaction network, process noise and sampling under non-equilibrium conditions. The sign of a correlation mostly coincided with the nature of the strongest pairwise interaction, but this is not necessarily the case. For instance, under rare conditions of identical interaction strength, we found that competitive and exploitative interactions can result in positive as well as negative correlations. Thus, cross-sectional abundance data carry limited information on specific interaction types. Correlations in abundance may hint at interactions but require independent validation.

Competitive endogenous RNA (ceRNA) represents a novel mechanism of gene regulation that controls several biological and pathological processes. Recently, an increasing number of in silico methods have been developed to accelerate the identification of such regulatory events. However, there is still a need for a tool supporting the hypothesis that ceRNA regulatory events only occur at specific miRNA expression levels. To this end, we present an R package, ceRNAR, which allows identification and analysis of ceRNA-miRNA triplets via integration of miRNA and RNA expression data. The ceRNAR package integrates three main steps: (i) identification of ceRNA pairs based on a rank-based correlation between pairs that considers the impact of miRNA and a running sum correlation statistic, (ii) sample clustering based on gene-gene correlation by circular binary segmentation, and (iii) peak merging to identify the most relevant sample patterns. In addition, ceRNAR also provides downstream analyses of identified ceRNA-miRNA triplets, including network analysis, functional annotation, survival analysis, external validation, and integration of different tools. The performance of our proposed approach was validated through simulation studies of different scenarios. Compared with several published tools, ceRNAR was able to identify true ceRNA triplets with high sensitivity, low false-positive rates, and acceptable running time. In real data applications, the ceRNAs common to two lung cancer datasets were identified in both datasets. The bridging miRNA for one of these, the ceRNA for MAP4K3, was identified by ceRNAR as hsa-let-7c-5p. Since similar cancer subtypes do share some biological patterns, these results demonstrated that our proposed algorithm was able to identify potential ceRNA targets in real patients. In summary, ceRNAR offers a novel algorithm and a comprehensive pipeline to identify and analyze ceRNA regulation. The package is implemented in R and is available on GitHub (https://github.com/ywhsiao/ceRNAR).

Activation of gene expression in response to environmental cues results in substantial phenotypic heterogeneity between cells that can impact a wide range of outcomes including differentiation, viral activation, and drug resistance. An important source of gene expression noise is transcriptional bursting, or the process by which transcripts are produced during infrequent bursts of promoter activity. Chromatin accessibility impacts transcriptional bursting by regulating the assembly of transcription factor and polymerase complexes on promoters, suggesting that the effect of an activating signal on transcriptional noise will depend on the initial chromatin state at the promoter. To explore this possibility, we simulated transcriptional activation using a transcriptional cycling model with three promoter states that represent chromatin remodeling, polymerase binding and pause release. We initiated this model over a large parameter range representing target genes with different chromatin environments, and found that, upon increasing the polymerase pause release rate to activate transcription, changes in gene expression noise varied significantly across initial promoter states. This model captured phenotypic differences in activation of latent HIV viruses integrated at different chromatin locations and mediated by the transcription factor NF-κB. Activating transcription in the model via increasing one or more of the transcript production rates, as occurs following NF-κB activation, reproduced experimentally measured transcript distributions for four different latent HIV viruses, as well as the bimodal pattern of HIV protein expression that leads to a subset of reactivated virus. Importantly, the parameter 'activation path' differentially affected gene expression noise, and ultimately viral activation, in line with experimental observations. This work demonstrates how upstream signaling pathways can be connected to biological processes that underlie transcriptional bursting, resulting in target gene-specific noise profiles following stimulation of a single upstream pathway.

Simulating complex biological models consisting of multiple ordinary differential equations can aid in the prediction of the pharmacological/biological responses; however, they are often hampered by the availability of reliable kinetic parameters. In the present study, we aimed to discover the properties of behaviors without determining an optimal combination of kinetic parameter values (parameter set). The key idea was to collect as many parameter sets as possible. Given that many systems are biologically stable and resilient (BSR), we focused on the dynamics around the steady state and formulated objective functions for BSR by partial linear approximation of the focused region. Using the objective functions and modified global cluster Newton method, we developed an algorithm for a thorough exploration of the allowable parameter space for biological systems (TEAPS). We first applied TEAPS to the NF-κB signaling model. This system shows a damped oscillation after stimulation and seems to fit the BSR constraint. By applying TEAPS, we found several directions in parameter space which stringently determines the BSR property. In such directions, the experimentally fitted parameter values were included in the range of the obtained parameter sets. The arachidonic acid metabolic pathway model was used as a model related to pharmacological responses. The pharmacological effects of nonsteroidal anti-inflammatory drugs were simulated using the parameter sets obtained by TEAPS. The structural properties of the system were partly extracted by analyzing the distribution of the obtained parameter sets. In addition, the simulations showed inter-drug differences in prostacyclin to thromboxane A2 ratio such that aspirin treatment tends to increase the ratio, while rofecoxib treatment tends to decrease it. These trends are comparable to the clinical observations. These results on real biological models suggest that the parameter sets satisfying the BSR condition can help in finding biologically plausible parameter sets and understanding the properties of biological systems.

Intrinsically disordered proteins (IDPs) are highly dynamic systems that play an important role in cell signaling processes and their misfunction often causes human disease. Proper understanding of IDP function not only requires the realistic characterization of their three-dimensional conformational ensembles at atomic-level resolution but also of the time scales of interconversion between their conformational substates. Large sets of experimental data are often used in combination with molecular modeling to restrain or bias models to improve agreement with experiment. It is shown here for the N-terminal transactivation domain of p53 (p53TAD) and Pup, which are two IDPs that fold upon binding to their targets, how the latest advancements in molecular dynamics (MD) simulations methodology produces native conformational ensembles by combining replica exchange with series of microsecond MD simulations. They closely reproduce experimental data at the global conformational ensemble level, in terms of the distribution properties of the radius of gyration tensor, and at the local level, in terms of NMR properties including 15N spin relaxation, without the need for reweighting. Further inspection revealed that 10-20% of the individual MD trajectories display the formation of secondary structures not observed in the experimental NMR data. The IDP ensembles were analyzed by graph theory to identify dominant inter-residue contact clusters and characteristic amino-acid contact propensities. These findings indicate that modern MD force fields with residue-specific backbone potentials can produce highly realistic IDP ensembles sampling a hierarchy of nano- and picosecond time scales providing new insights into their biological function.

Prior work on abiogenesis, the emergence of life from non-life, suggests that it requires chemical reaction networks that contain self-amplifying motifs, namely, autocatalytic cores. However, little is known about how the presence of multiple autocatalytic cores might allow for the gradual accretion of complexity on the path to life. To explore this problem, we develop the concept of a seed-dependent autocatalytic system (SDAS), which is a subnetwork that can autocatalytically self-maintain given a flux of food, but cannot be initiated by food alone. Rather, initiation of SDASs requires the transient introduction of chemical "seeds." We show that, depending on the topological relationship of SDASs in a chemical reaction network, a food-driven system can accrete complexity in a historically contingent manner, governed by rare seeding events. We develop new algorithms for detecting and analyzing SDASs in chemical reaction databases and describe parallels between multi-SDAS networks and biological ecosystems. Applying our algorithms to both an abiotic reaction network and a biochemical one, each driven by a set of simple food chemicals, we detect SDASs that are organized as trophic tiers, of which the higher tier can be seeded by relatively simple chemicals if the lower tier is already activated. This indicates that sequential activation of trophically organized SDASs by seed chemicals that are not much more complex than what already exist could be a mechanism of gradual complexification from relatively simple abiotic reactions to more complex life-like systems. Interestingly, in both reaction networks, higher-tier SDASs include chemicals that might alter emergent features of chemical systems and could serve as early targets of selection. Our analysis provides computational tools for analyzing very large chemical/biochemical reaction networks and suggests new approaches to studying abiogenesis in the lab.

Hemiparesis, defined as unilateral muscle weakness, often occurs in people post-stroke or people with cerebral palsy, however it is difficult to understand how this hemiparesis affects movement patterns as it often presents alongside a variety of other neuromuscular impairments. Predictive musculoskeletal modeling presents an opportunity to investigate how impairments affect gait performance assuming a particular cost function. Here, we use predictive simulation to quantify the spatiotemporal asymmetries and changes to metabolic cost that emerge when muscle strength is unilaterally reduced and how reducing spatiotemporal symmetry affects metabolic cost. We modified a 2-D musculoskeletal model by uniformly reducing the peak isometric muscle force unilaterally. We then solved optimal control simulations of walking across a range of speeds by minimizing the sum of the cubed muscle excitations. Lastly, we ran additional optimizations to test if reducing spatiotemporal asymmetry would result in an increase in metabolic cost. Our results showed that the magnitude and direction of effort-optimal spatiotemporal asymmetries depends on both the gait speed and level of weakness. Also, the optimal speed was 1.25 m/s for the symmetrical and 20% weakness models but slower (1.00 m/s) for the 40% and 60% weakness models, suggesting that hemiparesis can account for a portion of the slower gait speed seen in people with hemiparesis. Modifying the cost function to minimize spatiotemporal asymmetry resulted in small increases (~4%) in metabolic cost. Overall, our results indicate that spatiotemporal asymmetry may be optimal for people with hemiparesis. Additionally, the effect of speed and the level of weakness on spatiotemporal asymmetry may help explain the well-known heterogenous distribution of spatiotemporal asymmetries observed in the clinic. Future work could extend our results by testing the effects of other neuromuscular impairments on optimal gait strategies, and therefore build a more comprehensive understanding of the gait patterns observed in clinical populations.

Sustainable research on computational models of neuronal networks requires published models to be understandable, reproducible, and extendable. Missing details or ambiguities about mathematical concepts and assumptions, algorithmic implementations, or parameterizations hinder progress. Such flaws are unfortunately frequent and one reason is a lack of readily applicable standards and tools for model description. Our work aims to advance complete and concise descriptions of network connectivity but also to guide the implementation of connection routines in simulation software and neuromorphic hardware systems. We first review models made available by the computational neuroscience community in the repositories ModelDB and Open Source Brain, and investigate the corresponding connectivity structures and their descriptions in both manuscript and code. The review comprises the connectivity of networks with diverse levels of neuroanatomical detail and exposes how connectivity is abstracted in existing description languages and simulator interfaces. We find that a substantial proportion of the published descriptions of connectivity is ambiguous. Based on this review, we derive a set of connectivity concepts for deterministically and probabilistically connected networks and also address networks embedded in metric space. Beside these mathematical and textual guidelines, we propose a unified graphical notation for network diagrams to facilitate an intuitive understanding of network properties. Examples of representative network models demonstrate the practical use of the ideas. We hope that the proposed standardizations will contribute to unambiguous descriptions and reproducible implementations of neuronal network connectivity in computational neuroscience.

Improving biological plausibility and functional capacity are two important goals for brain models that connect low-level neural details to high-level behavioral phenomena. We develop a method called "oracle-supervised Neural Engineering Framework" (osNEF) to train biologically-detailed spiking neural networks that realize a variety of cognitively-relevant dynamical systems. Specifically, we train networks to perform computations that are commonly found in cognitive systems (communication, multiplication, harmonic oscillation, and gated working memory) using four distinct neuron models (leaky-integrate-and-fire neurons, Izhikevich neurons, 4-dimensional nonlinear point neurons, and 4-compartment, 6-ion-channel layer-V pyramidal cell reconstructions) connected with various synaptic models (current-based synapses, conductance-based synapses, and voltage-gated synapses). We show that osNEF networks exhibit the target dynamics by accounting for nonlinearities present within the neuron models: performance is comparable across all four systems and all four neuron models, with variance proportional to task and neuron model complexity. We also apply osNEF to build a model of working memory that performs a delayed response task using a combination of pyramidal cells and inhibitory interneurons connected with NMDA and GABA synapses. The baseline performance and forgetting rate of the model are consistent with animal data from delayed match-to-sample tasks (DMTST): we observe a baseline performance of 95% and exponential forgetting with time constant τ = 8.5s, while a recent meta-analysis of DMTST performance across species observed baseline performances of 58 - 99% and exponential forgetting with time constants of τ = 2.4 - 71s. These results demonstrate that osNEF can train functional brain models using biologically-detailed components and open new avenues for investigating the relationship between biophysical mechanisms and functional capabilities.

The development of increasingly efficient and cost-effective high throughput DNA sequencing techniques has enhanced the possibility of studying complex microbial systems. Recently, researchers have shown great interest in studying the microorganisms that characterise different ecological niches. Differential abundance analysis aims to find the differences in the abundance of each taxa between two classes of subjects or samples, assigning a significance value to each comparison. Several bioinformatic methods have been specifically developed, taking into account the challenges of microbiome data, such as sparsity, the different sequencing depth constraint between samples and compositionality. Differential abundance analysis has led to important conclusions in different fields, from health to the environment. However, the lack of a known biological truth makes it difficult to validate the results obtained. In this work we exploit metaSPARSim, a microbial sequencing count data simulator, to simulate data with differential abundance features between experimental groups. We perform a complete comparison of recently developed and established methods on a common benchmark with great effort to the reliability of both the simulated scenarios and the evaluation metrics. The performance overview includes the investigation of numerous scenarios, studying the effect on methods' results on the main covariates such as sample size, percentage of differentially abundant features, sequencing depth, feature variability, normalisation approach and ecological niches. Mainly, we find that methods show a good control of the type I error and, generally, also of the false discovery rate at high sample size, while recall seem to depend on the dataset and sample size.

Models of consumer effects on a shared resource environment have helped clarify how the interplay of consumer traits and resource supply impact stable coexistence. Recent models generalize this picture to include the exchange of resources alongside resource competition. These models exemplify the fact that although consumers shape the resource environment, the outcome of consumer interactions is context-dependent: such models can have either stable or unstable equilibria, depending on the resource supply. However, these recent models focus on a simplified version of microbial metabolism where the depletion of resources always leads to consumer growth. Here, we model an arbitrarily large system of consumers governed by Liebig's law, where species require and deplete multiple resources, but each consumer's growth rate is only limited by a single one of these resources. Resources that are taken up but not incorporated into new biomass are leaked back into the environment, possibly transformed by intracellular reactions, thereby tying the mismatch between depletion and growth to cross-feeding. For this set of dynamics, we show that feasible equilibria can be either stable or unstable, again depending on the resource environment. We identify special consumption and production networks which protect the community from instability when resources are scarce. Using simulations, we demonstrate that the qualitative stability patterns derived analytically apply to a broader class of network structures and resource inflow profiles, including cases where multiple species coexist on only one externally supplied resource. Our stability criteria bear some resemblance to classic stability results for pairwise interactions, but also demonstrate how environmental context can shape coexistence patterns when resource limitation and exchange are modeled directly.

Measles is one the best-documented and most-mechanistically-studied non-linear infectious disease dynamical systems. However, systematic investigation into the comparative performance of traditional mechanistic models and machine learning approaches in forecasting the transmission dynamics of this pathogen are still rare. Here, we compare one of the most widely used semi-mechanistic models for measles (TSIR) with a commonly used machine learning approach (LASSO), comparing performance and limits in predicting short to long term outbreak trajectories and seasonality for both regular and less regular measles outbreaks in England and Wales (E&W) and the United States. First, our results indicate that the proposed LASSO model can efficiently use data from multiple major cities and achieve similar short-to-medium term forecasting performance to semi-mechanistic models for E&W epidemics. Second, interestingly, the LASSO model also captures annual to biennial bifurcation of measles epidemics in E&W caused by susceptible response to the late 1940s baby boom. LASSO may also outperform TSIR for predicting less-regular dynamics such as those observed in major cities in US between 1932-45. Although both approaches capture short-term forecasts, accuracy suffers for both methods as we attempt longer-term predictions in highly irregular, post-vaccination outbreaks in E&W. Finally, we illustrate that the LASSO model can both qualitatively and quantitatively reconstruct mechanistic assumptions, notably susceptible dynamics, in the TSIR model. Our results characterize the limits of predictability of infectious disease dynamics for strongly immunizing pathogens with both mechanistic and machine learning models, and identify connections between these two approaches.

Constraint-based modeling is a powerful framework for studying cellular metabolism, with applications ranging from predicting growth rates and optimizing production of high value metabolites to identifying enzymes in pathogens that may be targeted for therapeutic interventions. Results from modeling experiments can be affected at least in part by the quality of the metabolic models used. Reconstructing a metabolic network manually can produce a high-quality metabolic model but is a time-consuming task. At the same time, current methods for automating the process typically transfer metabolic function based on sequence similarity, a process known to produce many false positives. We created Architect, a pipeline for automatic metabolic model reconstruction from protein sequences. First, it performs enzyme annotation through an ensemble approach, whereby a likelihood score is computed for an EC prediction based on predictions from existing tools; for this step, our method shows both increased precision and recall compared to individual tools. Next, Architect uses these annotations to construct a high-quality metabolic network which is then gap-filled based on likelihood scores from the ensemble approach. The resulting metabolic model is output in SBML format, suitable for constraints-based analyses. Through comparisons of enzyme annotations and curated metabolic models, we demonstrate improved performance of Architect over other state-of-the-art tools, notably with higher precision and recall on the eukaryote C. elegans and when compared to UniProt annotations in two bacterial species. Code for Architect is available at https://github.com/ParkinsonLab/Architect. For ease-of-use, Architect can be readily set up and utilized using its Docker image, maintained on Docker Hub.

Genetic risk for complex traits is strongly enriched in non-coding genomic regions involved in gene regulation, especially enhancers. However, we lack adequate tools to connect the characteristics of these disruptions to genetic risk. Here, we propose RWAS (Regulome Wide Association Study), a new application of the MAGMA software package to identify the characteristics of enhancers that contribute to genetic risk for disease. RWAS involves three steps: (i) assign genotyped SNPs to cell type- or tissue-specific regulatory features (e.g., enhancers); (ii) test associations of each regulatory feature with a trait of interest for which genome-wide association study (GWAS) summary statistics are available; (iii) perform enhancer-set enrichment analyses to identify quantitative or categorical features of regulatory elements that are associated with the trait. These steps are implemented as a novel application of MAGMA, a tool originally developed for gene-based GWAS analyses. Applying RWAS to interrogate genetic risk for schizophrenia, we discovered a class of risk-associated AT-rich enhancers that are active in the developing brain and harbor binding sites for multiple transcription factors with neurodevelopmental functions. RWAS utilizes open-source software, and we provide a comprehensive collection of annotations for tissue-specific enhancer locations and features, including their evolutionary conservation, AT content, and co-localization with binding sites for hundreds of TFs. RWAS will enable researchers to characterize properties of regulatory elements associated with any trait of interest for which GWAS summary statistics are available.

Drug-induced toxicity damages the health and is one of the key factors causing drug withdrawal from the market. It is of great significance to identify drug-induced target-organ toxicity, especially the detailed pathological findings, which are crucial for toxicity assessment, in the early stage of drug development process. A large variety of studies have devoted to identify drug toxicity. However, most of them are limited to single organ or only binary toxicity. Here we proposed a novel multi-label learning model named Att-RethinkNet, for predicting drug-induced pathological findings targeted on liver and kidney based on toxicogenomics data. The Att-RethinkNet is equipped with a memory structure and can effectively use the label association information. Besides, attention mechanism is embedded to focus on the important features and obtain better feature presentation. Our Att-RethinkNet is applicable in multiple organs and takes account the compound type, dose, and administration time, so it is more comprehensive and generalized. And more importantly, it predicts multiple pathological findings at the same time, instead of predicting each pathology separately as the previous model did. To demonstrate the effectiveness of the proposed model, we compared the proposed method with a series of state-of-the-arts methods. Our model shows competitive performance and can predict potential hepatotoxicity and nephrotoxicity in a more accurate and reliable way. The implementation of the proposed method is available at https://github.com/RanSuLab/Drug-Toxicity-Prediction-MultiLabel.

ATP7B is a human copper-transporting P1B-type ATPase that is involved in copper homeostasis and resistance to platinum drugs in cancer cells. ATP7B consists of a copper-transporting core and a regulatory N-terminal tail that contains six metal-binding domains (MBD1-6) connected by linker regions. The MBDs can bind copper, which changes the dynamics of the regulatory domain and activates the protein, but the underlying mechanism remains unknown. To identify possible copper-specific structural dynamics involved in transport regulation, we constructed a model of ATP7B spanning the N-terminal tail and core catalytic domains and performed molecular dynamics (MD) simulations with (holo) and without (apo) copper ions bound to the MBDs. In the holo protein, MBD2, MBD3 and MBD5 showed enhanced mobilities, which resulted in a more extended N-terminal regulatory region. The observed separation of MBD2 and MBD3 from the core protein supports a mechanism where copper binding activates the ATP7B protein by reducing interactions among MBD1-3 and between MBD1-3 and the core protein. We also observed an increased interaction between MBD5 and the core protein that brought the copper-binding site of MBD5 closer to the high-affinity internal copper-binding site in the core protein. The simulation results assign specific, mechanistic roles to the metal-binding domains involved in ATP7B regulation that are testable in experimental settings.

The first year of the COVID-19 pandemic put considerable strain on healthcare systems worldwide. In order to predict the effect of the local epidemic on hospital capacity in England, we used a variety of data streams to inform the construction and parameterisation of a hospital progression model, EpiBeds, which was coupled to a model of the generalised epidemic. In this model, individuals progress through different pathways (e.g. may recover, die, or progress to intensive care and recover or die) and data from a partially complete patient-pathway line-list was used to provide initial estimates of the mean duration that individuals spend in the different hospital compartments. We then fitted EpiBeds using complete data on hospital occupancy and hospital deaths, enabling estimation of the proportion of individuals that follow the different clinical pathways, the reproduction number of the generalised epidemic, and to make short-term predictions of hospital bed demand. The construction of EpiBeds makes it straightforward to adapt to different patient pathways and settings beyond England. As part of the UK response to the pandemic, EpiBeds provided weekly forecasts to the NHS for hospital bed occupancy and admissions in England, Wales, Scotland, and Northern Ireland at national and regional scales.

The self-organising global dynamics underlying brain states emerge from complex recursive nonlinear interactions between interconnected brain regions. Until now, most efforts of capturing the causal mechanistic generating principles have supposed underlying stationarity, being unable to describe the non-stationarity of brain dynamics, i.e. time-dependent changes. Here, we present a novel framework able to characterise brain states with high specificity, precisely by modelling the time-dependent dynamics. Through describing a topological structure associated to the brain state at each moment in time (its attractor or 'information structure'), we are able to classify different brain states by using the statistics across time of these structures hitherto hidden in the neuroimaging dynamics. Proving the strong potential of this framework, we were able to classify resting-state BOLD fMRI signals from two classes of post-comatose patients (minimally conscious state and unresponsive wakefulness syndrome) compared with healthy controls with very high precision.

When the mixture solution of cyanobacterial proteins, KaiA, KaiB, and KaiC, is incubated with ATP in vitro, the phosphorylation level of KaiC shows stable oscillations with the temperature-compensated circadian period. Elucidating this temperature compensation is essential for understanding the KaiABC circadian clock, but its mechanism has remained a mystery. We analyzed the KaiABC temperature compensation by developing a theoretical model describing the feedback relations among reactions and structural transitions in the KaiC molecule. The model showed that the reduced structural cooperativity should weaken the negative feedback coupling among reactions and structural transitions, which enlarges the oscillation amplitude and period, explaining the observed significant period extension upon single amino-acid residue substitution. We propose that an increase in thermal fluctuations similarly attenuates the reaction-structure feedback, explaining the temperature compensation in the KaiABC clock. The model explained the experimentally observed responses of the oscillation phase to the temperature shift or the ADP-concentration change and suggested that the ATPase reactions in the CI domain of KaiC affect the period depending on how the reaction rates are modulated. The KaiABC clock provides a unique opportunity to analyze how the reaction-structure coupling regulates the system-level synchronized oscillations of molecules.

Comprehensive molecular characterization of cancer subtypes is essential for predicting clinical outcomes and searching for personalized treatments. We present bnClustOmics, a statistical model and computational tool for multi-omics unsupervised clustering, which serves a dual purpose: Clustering patient samples based on a Bayesian network mixture model and learning the networks of omics variables representing these clusters. The discovered networks encode interactions among all omics variables and provide a molecular characterization of each patient subgroup. We conducted simulation studies that demonstrated the advantages of our approach compared to other clustering methods in the case where the generative model is a mixture of Bayesian networks. We applied bnClustOmics to a hepatocellular carcinoma (HCC) dataset comprising genome (mutation and copy number), transcriptome, proteome, and phosphoproteome data. We identified three main HCC subtypes together with molecular characteristics, some of which are associated with survival even when adjusting for the clinical stage. Cluster-specific networks shed light on the links between genotypes and molecular phenotypes of samples within their respective clusters and suggest targets for personalized treatments.

BACKGROUND: Based on the principles of equity and effectiveness, the World Health Organization and COVAX formulate vaccine allocation as a mathematical optimization problem. This study aims to solve the optimization problem using agent-based simulations. METHODS: We built open-sourced agent-based models to simulate virus transition among a demographically representative sample of 198 million people in 148 countries using advanced computational services. All countries continuing their current vaccine progress is defined as the baseline scenario. Comparison scenarios include achieving minimum vaccination rates and allocating vaccines based on pandemic levels. FINDINGS: The simulations are fitted using the pandemic data from 148 countries from January 2020 to June 2021. Under the baseline scenario, the world will add 24.36 million cases and 468,945 deaths during the projection period of three months. Inoculating at least 10%, 20%, and 26% of populations in all countries requires 1.12, 3.31, and 5.00 million additional vaccine doses every day, respectively. Achieving these benchmarks reduces new cases by 0.56, 2.74, and 3.32 million, respectively. If allocated by the current global distribution, 5.00 million additional vaccine doses will only avert 1.45 million new cases. If those 5.00 million vaccines are allocated based on projected cases in each country, the averted cases will increase more than six-fold to 9.20 million. Similar differences between allocation methods are observed in averted deaths. CONCLUSION: The global distribution of COVID-19 vaccines can be optimized to achieve better outcomes in terms of both equity and effectiveness. Alternative vaccine allocation methods may avert several times more cases and deaths than the current global distribution. With reasonable requirements on additional vaccines, COVAX could adopt alternative allocation strategies that reduce cross-country inequity and save more lives.

Convolutional neural networks trained on object recognition derive inspiration from the neural architecture of the visual system in mammals, and have been used as models of the feedforward computation performed in the primate ventral stream. In contrast to the deep hierarchical organization of primates, the visual system of the mouse has a shallower arrangement. Since mice and primates are both capable of visually guided behavior, this raises questions about the role of architecture in neural computation. In this work, we introduce a novel framework for building a biologically constrained convolutional neural network model of the mouse visual cortex. The architecture and structural parameters of the network are derived from experimental measurements, specifically the 100-micrometer resolution interareal connectome, the estimates of numbers of neurons in each area and cortical layer, and the statistics of connections between cortical layers. This network is constructed to support detailed task-optimized models of mouse visual cortex, with neural populations that can be compared to specific corresponding populations in the mouse brain. Using a well-studied image classification task as our working example, we demonstrate the computational capability of this mouse-sized network. Given its relatively small size, MouseNet achieves roughly 2/3rds the performance level on ImageNet as VGG16. In combination with the large scale Allen Brain Observatory Visual Coding dataset, we use representational similarity analysis to quantify the extent to which MouseNet recapitulates the neural representation in mouse visual cortex. Importantly, we provide evidence that optimizing for task performance does not improve similarity to the corresponding biological system beyond a certain point. We demonstrate that the distributions of some physiological quantities are closer to the observed distributions in the mouse brain after task training. We encourage the use of the MouseNet architecture by making the code freely available.

Robustness in developing and homeostatic tissues is supported by various types of spatiotemporal cell-to-cell interactions. Although live imaging and cell tracking are powerful in providing direct evidence of cell coordination rules, extracting and comparing these rules across many tissues with potentially different length and timescales of coordination requires a versatile framework of analysis. Here we demonstrate that graph neural network (GNN) models are suited for this purpose, by showing how they can be applied to predict cell fate in tissues and utilized to infer the cell interactions governing the multicellular dynamics. Analyzing the live mammalian epidermis data, where spatiotemporal graphs constructed from cell tracks and cell contacts are given as inputs, GNN discovers distinct neighbor cell fate coordination rules that depend on the region of the body. This approach demonstrates how the GNN framework is powerful in inferring general cell interaction rules from live data without prior knowledge of the signaling involved.

The widespread, and in many countries unprecedented, use of non-pharmaceutical interventions (NPIs) during the COVID-19 pandemic has highlighted the need for mathematical models which can estimate the impact of these measures while accounting for the highly heterogeneous risk profile of COVID-19. Models accounting either for age structure or the household structure necessary to explicitly model many NPIs are commonly used in infectious disease modelling, but models incorporating both levels of structure present substantial computational and mathematical challenges due to their high dimensionality. Here we present a modelling framework for the spread of an epidemic that includes explicit representation of age structure and household structure. Our model is formulated in terms of tractable systems of ordinary differential equations for which we provide an open-source Python implementation. Such tractability leads to significant benefits for model calibration, exhaustive evaluation of possible parameter values, and interpretability of results. We demonstrate the flexibility of our model through four policy case studies, where we quantify the likely benefits of the following measures which were either considered or implemented in the UK during the current COVID-19 pandemic: control of within- and between-household mixing through NPIs; formation of support bubbles during lockdown periods; out-of-household isolation (OOHI); and temporary relaxation of NPIs during holiday periods. Our ordinary differential equation formulation and associated analysis demonstrate that multiple dimensions of risk stratification and social structure can be incorporated into infectious disease models without sacrificing mathematical tractability. This model and its software implementation expand the range of tools available to infectious disease policy analysts.

Aging involves a transition from youthful vigor to geriatric infirmity and death. Individuals who remain vigorous longer tend to live longer, and within isogenic populations of C. elegans the timing of age-associated vigorous movement cessation (VMC) is highly correlated with lifespan. Yet, many mutations and interventions in aging alter the proportion of lifespan spent moving vigorously, appearing to "uncouple" youthful vigor from lifespan. To clarify the relationship between vigorous movement cessation, death, and the physical declines that determine their timing, we developed a new version of the imaging platform called "The Lifespan Machine". This technology allows us to compare behavioral aging and lifespan at an unprecedented scale. We find that behavioral aging involves a time-dependent increase in the risk of VMC, reminiscent of the risk of death. Furthermore, we find that VMC times are inversely correlated with remaining lifespan across a wide range of genotypes and environmental conditions. Measuring and modelling a variety of lifespan-altering interventions including a new RNA-polymerase II auxin-inducible degron system, we find that vigorous movement and lifespan are best described as emerging from the interplay between at least two distinct physical declines whose rates co-vary between individuals. In this way, we highlight a crucial limitation of predictors of lifespan like VMC-in organisms experiencing multiple, distinct, age-associated physical declines, correlations between mid-life biomarkers and late-life outcomes can arise from the contextual influence of confounding factors rather than a reporting by the biomarker of a robustly predictive biological age.

Recent state-of-the-art multiplex imaging techniques have expanded the depth of information that can be captured within a single tissue sample by allowing for panels with dozens of markers. Despite this increase in capacity, space on the panel is still limited due to technical artifacts, tissue loss, and long imaging acquisition time. As such, selecting which markers to include on a panel is important, since removing important markers will result in a loss of biologically relevant information, but identifying redundant markers will provide a room for other markers. To address this, we propose computational approaches to determine the amount of shared information between markers and select an optimally reduced panel that captures maximum amount of information with the fewest markers. Here we examine several panel selection approaches and evaluate them based on their ability to reconstruct the full panel images and information within breast cancer tissue microarray datasets using cyclic immunofluorescence as a proof of concept. We show that all methods perform adequately and can re-capture cell types using only 18 of 25 markers (72% of the original panel size). The correlation-based selection methods achieved the best single-cell marker mean intensity predictions with a Spearman correlation of 0.90 with the reduced panel. Using the proposed methods shown here, it is possible for researchers to design more efficient multiplex imaging panels that maximize the amount of information retained with the limited number of markers with respect to certain evaluation metrics and architecture biases.

Despite increased understanding about psoriasis pathophysiology, currently there is a lack of predictive computational models. We developed a personalisable ordinary differential equations model of human epidermis and psoriasis that incorporates immune cells and cytokine stimuli to regulate the transition between two stable steady states of clinically healthy (non-lesional) and disease (lesional psoriasis, plaque) skin. In line with experimental data, an immune stimulus initiated transition from healthy skin to psoriasis and apoptosis of immune and epidermal cells induced by UVB phototherapy returned the epidermis back to the healthy state. Notably, our model was able to distinguish disease flares. The flexibility of our model permitted the development of a patient-specific "UVB sensitivity" parameter that reflected subject-specific sensitivity to apoptosis and enabled simulation of individual patients' clinical response trajectory. In a prospective clinical study of 94 patients, serial individual UVB doses and clinical response (Psoriasis Area Severity Index) values collected over the first three weeks of UVB therapy informed estimation of the "UVB sensitivity" parameter and the prediction of individual patient outcome at the end of phototherapy. An important advance of our model is its potential for direct clinical application through early assessment of response to UVB therapy, and for individualised optimisation of phototherapy regimes to improve clinical outcome. Additionally by incorporating the complex interaction of immune cells and epidermal keratinocytes, our model provides a basis to study and predict outcomes to biologic therapies in psoriasis.

Identification of specific species in metagenomic samples is critical for several key applications, yet many tools available require large computational power and are often prone to false positive identifications. Here we describe High-AccuracY and Scalable Taxonomic Assignment of MetagenomiC data (HAYSTAC), which can estimate the probability that a specific taxon is present in a metagenome. HAYSTAC provides a user-friendly tool to construct databases, based on publicly available genomes, that are used for competitive read mapping. It then uses a novel Bayesian framework to infer the abundance and statistical support for each species identification and provide per-read species classification. Unlike other methods, HAYSTAC is specifically designed to efficiently handle both ancient and modern DNA data, as well as incomplete reference databases, making it possible to run highly accurate hypothesis-driven analyses (i.e., assessing the presence of a specific species) on variably sized reference databases while dramatically improving processing speeds. We tested the performance and accuracy of HAYSTAC using simulated Illumina libraries, both with and without ancient DNA damage, and compared the results to other currently available methods (i.e., Kraken2/Bracken, KrakenUniq, MALT/HOPS, and Sigma). HAYSTAC identified fewer false positives than both Kraken2/Bracken, KrakenUniq and MALT in all simulations, and fewer than Sigma in simulations of ancient data. It uses less memory than Kraken2/Bracken, KrakenUniq as well as MALT both during database construction and sample analysis. Lastly, we used HAYSTAC to search for specific pathogens in two published ancient metagenomic datasets, demonstrating how it can be applied to empirical datasets. HAYSTAC is available from https://github.com/antonisdim/HAYSTAC.

A ubiquitous pattern in ecological systems is that more abundant species tend to be more generalist; that is, they interact with more species or can occur in wider range of habitats. However, there is no consensus on whether generalism drives abundance (a selection process) or abundance drives generalism (a drift process). As it is difficult to conduct direct experiments to solve this chicken-and-egg dilemma, previous studies have used a causal discovery method based on formal logic and have found that abundance drives generalism. Here, we refine this method by correcting its bias regarding skewed distributions, and employ two other independent causal discovery methods based on nonparametric regression and on information theory, respectively. Contrary to previous work, all three independent methods strongly indicate that generalism drives abundance when applied to datasets on plant-hummingbird communities and reef fishes. Furthermore, we find that selection processes are more important than drift processes in structuring multispecies systems when the environment is variable. Our results showcase the power of the computational causal discovery approach to aid ecological research.

Long-term and culture-specific experience of music shapes rhythm perception, leading to enculturated expectations that make certain rhythms easier to track and more conducive to synchronized movement. However, the influence of enculturated bias on the moment-to-moment dynamics of rhythm tracking is not well understood. Recent modeling work has formulated entrainment to rhythms as a formal inference problem, where phase is continuously estimated based on precise event times and their correspondence to timing expectations: PIPPET (Phase Inference from Point Process Event Timing). Here we propose that the problem of optimally tracking a rhythm also requires an ongoing process of inferring which pattern of event timing expectations is most suitable to predict a stimulus rhythm. We formalize this insight as an extension of PIPPET called pPIPPET (PIPPET with pattern inference). The variational solution to this problem introduces terms representing the likelihood that a stimulus is based on a particular member of a set of event timing patterns, which we initialize according to culturally-learned prior expectations of a listener. We evaluate pPIPPET in three experiments. First, we demonstrate that pPIPPET can qualitatively reproduce enculturated bias observed in human tapping data for simple two-interval rhythms. Second, we simulate categorization of a continuous three-interval rhythm space by Western-trained musicians through derivation of a comprehensive set of priors for pPIPPET from metrical patterns in a sample of Western rhythms. Third, we simulate iterated reproduction of three-interval rhythms, and show that models configured with notated rhythms from different cultures exhibit both universal and enculturated biases as observed experimentally in listeners from those cultures. These results suggest the influence of enculturated timing expectations on human perceptual and motor entrainment can be understood as approximating optimal inference about the rhythmic stimulus, with respect to prototypical patterns in an empirical sample of rhythms that represent the music-cultural environment of the listener.

In many biological populations, such as human groups, individuals face a complex strategic setting, where they need to make strategic decisions over a diverse set of issues and their behavior in one strategic context can affect their decisions in another. This raises the question of how the interaction between different strategic contexts affects individuals' strategic choices and social norms? To address this question, I introduce a framework where individuals play two games with different structures and decide upon their strategy in a second game based on their knowledge of their opponent's strategy in the first game. I consider both multistage games, where the same opponents play the two games consecutively, and reputation-based model, where individuals play their two games with different opponents but receive information about their opponent's strategy. By considering a case where the first game is a social dilemma, I show that when the second game is a coordination or anti-coordination game, the Nash equilibrium of the coupled game can be decomposed into two classes, a defective equilibrium which is composed of two simple equilibrium of the two games, and a cooperative equilibrium, in which coupling between the two games emerge and sustain cooperation in the social dilemma. For the existence of the cooperative equilibrium, the cost of cooperation should be smaller than a value determined by the structure of the second game. Investigation of the evolutionary dynamics shows that a cooperative fixed point exists when the second game belongs to coordination or anti-coordination class in a mixed population. However, the basin of attraction of the cooperative fixed point is much smaller for the coordination class, and this fixed point disappears in a structured population. When the second game belongs to the anti-coordination class, the system possesses a spontaneous symmetry-breaking phase transition above which the symmetry between cooperation and defection breaks. A set of cooperation supporting moral norms emerges according to which cooperation stands out as a valuable trait. Notably, the moral system also brings a more efficient allocation of resources in the second game. This observation suggests a moral system has two different roles: Promotion of cooperation, which is against individuals' self-interest but beneficial for the population, and promotion of organization and order, which is at both the population's and the individual's self-interest. Interestingly, the latter acts like a Trojan horse: Once established out of individuals' self-interest, it brings the former with itself. Importantly, the fact that the evolution of moral norms depends only on the cost of cooperation and is independent of the benefit of cooperation implies that moral norms can be harmful and incur a pure collective cost, yet they are just as effective in promoting order and organization. Finally, the model predicts that recognition noise can have a surprisingly positive effect on the evolution of moral norms and facilitates cooperation in the Snow Drift game in structured populations.

Selection protocols such as SELEX, where molecules are selected over multiple rounds for their ability to bind to a target of interest, are popular methods for obtaining binders for diagnostic and therapeutic purposes. We show that Restricted Boltzmann Machines (RBMs), an unsupervised two-layer neural network architecture, can successfully be trained on sequence ensembles from single rounds of SELEX experiments for thrombin aptamers. RBMs assign scores to sequences that can be directly related to their fitnesses estimated through experimental enrichment ratios. Hence, RBMs trained from sequence data at a given round can be used to predict the effects of selection at later rounds. Moreover, the parameters of the trained RBMs are interpretable and identify functional features contributing most to sequence fitness. To exploit the generative capabilities of RBMs, we introduce two different training protocols: one taking into account sequence counts, capable of identifying the few best binders, and another based on unique sequences only, generating more diverse binders. We then use RBMs model to generate novel aptamers with putative disruptive mutations or good binding properties, and validate the generated sequences with gel shift assay experiments. Finally, we compare the RBM's performance with different supervised learning approaches that include random forests and several deep neural network architectures.

Our understanding of the firing behaviour of motoneuron (MN) pools during human voluntary muscle contractions is currently limited to electrophysiological findings from animal experiments extrapolated to humans, mathematical models of MN pools not validated for human data, and experimental results obtained from decomposition of electromyographical (EMG) signals. These approaches are limited in accuracy or provide information on only small partitions of the MN population. Here, we propose a method based on the combination of high-density EMG (HDEMG) data and realistic modelling for predicting the behaviour of entire pools of motoneurons in humans. The method builds on a physiologically realistic model of a MN pool which predicts, from the experimental spike trains of a smaller number of individual MNs identified from decomposed HDEMG signals, the unknown recruitment and firing activity of the remaining unidentified MNs in the complete MN pool. The MN pool model is described as a cohort of single-compartment leaky fire-and-integrate (LIF) models of MNs scaled by a physiologically realistic distribution of MN electrophysiological properties and driven by a spinal synaptic input, both derived from decomposed HDEMG data. The MN spike trains and effective neural drive to muscle, predicted with this method, have been successfully validated experimentally. A representative application of the method in MN-driven neuromuscular modelling is also presented. The proposed approach provides a validated tool for neuroscientists, experimentalists, and modelers to infer the firing activity of MNs that cannot be observed experimentally, investigate the neuromechanics of human MN pools, support future experimental investigations, and advance neuromuscular modelling for investigating the neural strategies controlling human voluntary contractions.

A few years ago, it was proposed to use the simultaneous quantification of unspliced and spliced messenger RNA (mRNA) to add a temporal dimension to high-throughput snapshots of single cell RNA sequencing data. This concept can yield additional insight into the transcriptional dynamics of the biological systems under study. However, current methods for inferring cell state velocities from such data (known as RNA velocities) are afflicted by several theoretical and computational problems, hindering realistic and reliable velocity estimation. We discuss these issues and propose new solutions for addressing some of the current challenges in consistency of data processing, velocity inference and visualisation. We translate our computational conclusion in two velocity analysis tools: one detailed method κ-velo and one heuristic method eco-velo, each of which uses a different set of assumptions about the data.

Imaging neural activity in a behaving animal presents unique challenges in part because motion from an animal's movement creates artifacts in fluorescence intensity time-series that are difficult to distinguish from neural signals of interest. One approach to mitigating these artifacts is to image two channels simultaneously: one that captures an activity-dependent fluorophore, such as GCaMP, and another that captures an activity-independent fluorophore such as RFP. Because the activity-independent channel contains the same motion artifacts as the activity-dependent channel, but no neural signals, the two together can be used to identify and remove the artifacts. However, existing approaches for this correction, such as taking the ratio of the two channels, do not account for channel-independent noise in the measured fluorescence. Here, we present Two-channel Motion Artifact Correction (TMAC), a method which seeks to remove artifacts by specifying a generative model of the two channel fluorescence that incorporates motion artifact, neural activity, and noise. We use Bayesian inference to infer latent neural activity under this model, thus reducing the motion artifact present in the measured fluorescence traces. We further present a novel method for evaluating ground-truth performance of motion correction algorithms by comparing the decodability of behavior from two types of neural recordings; a recording that had both an activity-dependent fluorophore and an activity-independent fluorophore (GCaMP and RFP) and a recording where both fluorophores were activity-independent (GFP and RFP). A successful motion correction method should decode behavior from the first type of recording, but not the second. We use this metric to systematically compare five models for removing motion artifacts from fluorescent time traces. We decode locomotion from a GCaMP expressing animal 20x more accurately on average than from control when using TMAC inferred activity and outperforms all other methods of motion correction tested, the best of which were ~8x more accurate than control.

With the aid of laboratory typing techniques, infectious disease surveillance networks have the opportunity to obtain powerful information on the emergence, circulation, and evolution of multiple genotypes, serotypes or other subtypes of pathogens, informing understanding of transmission dynamics and strategies for prevention and control. The volume of typing performed on clinical isolates is typically limited by its ability to inform clinical care, cost and logistical constraints, especially in comparison with the capacity to monitor clinical reports of disease occurrence, which remains the most widespread form of public health surveillance. Viewing clinical disease reports as arising from a latent mixture of pathogen subtypes, laboratory typing of a subset of clinical cases can provide inference on the proportion of clinical cases attributable to each subtype (i.e., the mixture components). Optimizing protocols for the selection of isolates for typing by weighting specific subpopulations, locations, time periods, or case characteristics (e.g., disease severity), may improve inference of the frequency and distribution of pathogen subtypes within and between populations. Here, we apply the Disease Surveillance Informatics Optimization and Simulation (DIOS) framework to simulate and optimize hand foot and mouth disease (HFMD) surveillance in a high-burden region of western China. We identify laboratory surveillance designs that significantly outperform the existing network: the optimal network reduced mean absolute error in estimated serotype-specific incidence rates by 14.1%; similarly, the optimal network for monitoring severe cases reduced mean absolute error in serotype-specific incidence rates by 13.3%. In both cases, the optimal network designs achieved improved inference without increasing subtyping effort. We demonstrate how the DIOS framework can be used to optimize surveillance networks by augmenting clinical diagnostic data with limited laboratory typing resources, while adapting to specific, local surveillance objectives and constraints.

Insect flight is a strongly nonlinear and actuated dynamical system. As such, strategies for understanding its control have typically relied on either model-based methods or linearizations thereof. Here we develop a framework that combines model predictive control on an established flight dynamics model and deep neural networks (DNN) to create an efficient method for solving the inverse problem of flight control. We turn to natural systems for inspiration since they inherently demonstrate network pruning with the consequence of yielding more efficient networks for a specific set of tasks. This bio-inspired approach allows us to leverage network pruning to optimally sparsify a DNN architecture in order to perform flight tasks with as few neural connections as possible, however, there are limits to sparsification. Specifically, as the number of connections falls below a critical threshold, flight performance drops considerably. We develop sparsification paradigms and explore their limits for control tasks. Monte Carlo simulations also quantify the statistical distribution of network weights during pruning given initial random weights of the DNNs. We demonstrate that on average, the network can be pruned to retain a small amount of original network weights and still perform comparably to its fully-connected counterpart. The relative number of remaining weights, however, is highly dependent on the initial architecture and size of the network. Overall, this work shows that sparsely connected DNNs are capable of predicting the forces required to follow flight trajectories. Additionally, sparsification has sharp performance limits.

Genetic variations affect behavior and cause disease but understanding how these variants drive complex traits is still an open question. A common approach is to link the genetic variants to intermediate molecular phenotypes such as the transcriptome using RNA-sequencing (RNA-seq). Paradoxically, these variants between the samples are usually ignored at the beginning of RNA-seq analyses of many model organisms. This can skew the transcriptome estimates that are used later for downstream analyses, such as expression quantitative trait locus (eQTL) detection. Here, we assessed the impact of reference-based analysis on the transcriptome and eQTLs in a widely-used mouse genetic population: the BXD panel of recombinant inbred lines. We highlight existing reference bias in the transcriptome data analysis and propose practical solutions which combine available genetic variants, genotypes, and genome reference sequence. The use of custom BXD line references improved downstream analysis compared to classical genome reference. These insights would likely benefit genetic studies with a transcriptomic component and demonstrate that genome references need to be reassessed and improved.

After years of experience, humans become experts at perceiving letters. Is this visual capacity attained by learning specialized letter features, or by reusing general visual features previously learned in service of object categorization? To explore this question, we first measured the perceptual similarity of letters in two behavioral tasks, visual search and letter categorization. Then, we trained deep convolutional neural networks on either 26-way letter categorization or 1000-way object categorization, as a way to operationalize possible specialized letter features and general object-based features, respectively. We found that the general object-based features more robustly correlated with the perceptual similarity of letters. We then operationalized additional forms of experience-dependent letter specialization by altering object-trained networks with varied forms of letter training; however, none of these forms of letter specialization improved the match to human behavior. Thus, our findings reveal that it is not necessary to appeal to specialized letter representations to account for perceptual similarity of letters. Instead, we argue that it is more likely that the perception of letters depends on domain-general visual features.

Proximal genetic variants are frequently correlated, implying that the corresponding effect sizes detected by genome-wide association studies (GWAS) are also not independent. Methods already exist to account for this when aggregating effects from a single GWAS across genes or pathways. Here we present a rigorous yet fast method for detecting genes with coherent association signals for two traits, facilitating cross-GWAS analyses. To this end, we devised a new significance test for the covariance of datapoints not drawn independently but with a known inter-sample covariance structure. We show that the distribution of its test statistic is a linear combination of χ2 distributions with positive and negative coefficients. The corresponding cumulative distribution function can be efficiently calculated with Davies' algorithm at high precision. We apply this general framework to test for dependence between SNP-wise effect sizes of two GWAS at the gene level. We extend this test to detect also gene-wise causal links. We demonstrate the utility of our method by uncovering potential shared genetic links between the severity of COVID-19 and (1) being prescribed class M05B medication (drugs affecting bone structure and mineralization), (2) rheumatoid arthritis, (3) vitamin D (25OHD), and (4) serum calcium concentrations. Our method detects a potential role played by chemokine receptor genes linked to TH1 versus TH2 immune response, a gene related to integrin beta-1 cell surface expression, and other genes potentially impacting the severity of COVID-19. Our approach will be useful for similar analyses involving datapoints with known auto-correlation structures.

Transcription factors (TFs) are crucial for regulating cell differentiation during the development of the immune system. However, the key TFs for orchestrating the specification of distinct immune cells are not fully understood. Here, we integrated the transcriptomic and epigenomic measurements in 73 mouse and 61 human primary cell types, respectively, that span the immune cell differentiation pathways. We constructed the cell-type-specific transcriptional regulatory network and assessed the global importance of TFs based on the Taiji framework, which is a method we have previously developed that can infer the global impact of TFs using integrated transcriptomic and epigenetic data. Integrative analysis across cell types revealed putative driver TFs in cell lineage-specific differentiation in both mouse and human systems. We have also identified TF combinations that play important roles in specific developmental stages. Furthermore, we validated the functions of predicted novel TFs in murine CD8+ T cell differentiation and showed the importance of Elf1 and Prdm9 in the effector versus memory T cell fate specification and Kdm2b and Tet3 in promoting differentiation of CD8+ tissue resident memory (Trm) cells, validating the approach. Thus, we have developed a bioinformatic approach that provides a global picture of the regulatory mechanisms that govern cellular differentiation in the immune system and aids the discovery of novel mechanisms in cell fate decisions.

The rise of SARS-CoV-2 variants and the history of outbreaks caused by zoonotic coronaviruses point to the need for next-generation vaccines that confer protection against variant strains. Here, we combined analyses of diverse sequences and structures of coronavirus spikes with data from deep mutational scanning to design SARS-CoV-2 variant antigens containing the most significant mutations that may emerge. We trained a neural network to predict RBD expression and ACE2 binding from sequence, which allowed us to determine that these antigens are stable and bind to ACE2. Thus, they represent viable variants. We then used a computational model of affinity maturation (AM) to study the antibody response to immunization with different combinations of the designed antigens. The results suggest that immunization with a cocktail of the antigens is likely to promote evolution of higher titers of antibodies that target SARS-CoV-2 variants than immunization or infection with the wildtype virus alone. Finally, our analysis of 12 coronaviruses from different genera identified the S2' cleavage site and fusion peptide as potential pan-coronavirus vaccine targets.

Fluorescence Recovery After Photobleaching (FRAP) has been extensively used to understand molecular dynamics in cells. This technique when applied to soluble, globular molecules driven by diffusion is easily interpreted and well understood. However, the classical methods of analysis cannot be applied to anisotropic structures subjected to directed transport, such as cytoskeletal filaments or elongated organelles transported along microtubule tracks. A new mathematical approach is needed to analyze FRAP data in this context and determine what information can be obtain from such experiments. To address these questions, we analyze fluorescence intensity profile curves after photobleaching of fluorescently labelled intermediate filaments anterogradely transported along microtubules. We apply the analysis to intermediate filament data to determine information about the filament motion. Our analysis consists of deriving equations for fluorescence intensity profiles and developing a mathematical model for the motion of filaments and simulating the model. Two closed forms for profile curves were derived, one for filaments of constant length and one for filaments with constant velocity, and three types of simulation were carried out. In the first type of simulation, the filaments have random velocities which are constant for the duration of the simulation. In the second type, filaments have random velocities which instantaneously change at random times. In the third type, filaments have random velocities and exhibit pausing between velocity changes. Our analysis shows: the most important distribution governing the shape of the intensity profile curves obtained from filaments is the distribution of the filament velocity. Furthermore, filament length which is constant during the experiment, had little impact on intensity profile curves. Finally, gamma distributions for the filament velocity with pauses give the best fit to asymmetric fluorescence intensity profiles of intermediate filaments observed in FRAP experiments performed in polarized migrating astrocytes. Our analysis also shows that the majority of filaments are stationary. Overall, our data give new insight into the regulation of intermediate filament dynamics during cell migration.

Inferring causes of the good and bad events that we experience is part of the process of building models of our own capabilities and of the world around us. Making such inferences can be difficult because of complex reciprocal relationships between attributions of the causes of particular events, and beliefs about the capabilities and skills that influence our role in bringing them about. Abnormal causal attributions have long been studied in connection with psychiatric disorders, notably depression and paranoia; however, the mechanisms behind attributional inferences and the way they can go awry are not fully understood. We administered a novel, challenging, game of skill to a substantial population of healthy online participants, and collected trial-by-trial time series of both their beliefs about skill and attributions about the causes of the success and failure of real experienced outcomes. We found reciprocal relationships that provide empirical confirmation of the attribution-self representation cycle theory. This highlights the dynamic nature of the processes involved in attribution, and validates a framework for developing and testing computational accounts of attribution-belief interactions.

Previous studies for cancer biomarker discovery based on pre-diagnostic blood DNA methylation (DNAm) profiles, either ignore the explicit modeling of the Time To Diagnosis (TTD), or provide inconsistent results. This lack of consistency is likely due to the limitations of standard EWAS approaches, that model the effect of DNAm at CpG sites on TTD independently. In this work, we aim to identify blood DNAm profiles associated with TTD, with the aim to improve the reliability of the results, as well as their biological meaningfulness. We argue that a global approach to estimate CpG sites effect profile should capture the complex (potentially non-linear) relationships interplaying between sites. To prove our concept, we develop a new Deep Learning-based approach assessing the relevance of individual CpG Islands (i.e., assigning a weight to each site) in determining TTD while modeling their combined effect in a survival analysis scenario. The algorithm combines a tailored sampling procedure with DNAm sites agglomeration, deep non-linear survival modeling and SHapley Additive exPlanations (SHAP) values estimation to aid robustness of the derived effects profile. The proposed approach deals with the common complexities arising from epidemiological studies, such as small sample size, noise, and low signal-to-noise ratio of blood-derived DNAm. We apply our approach to a prospective case-control study on breast cancer nested in the EPIC Italy cohort and we perform weighted gene-set enrichment analyses to demonstrate the biological meaningfulness of the obtained results. We compared the results of Deep Survival EWAS with those of a traditional EWAS approach, demonstrating that our method performs better than the standard approach in identifying biologically relevant pathways.

The metagenome embedded in urban sewage is an attractive new data source to understand urban ecology and assess human health status at scales beyond a single host. Analyzing the viral fraction of wastewater in the ongoing COVID-19 pandemic has shown the potential of wastewater as aggregated samples for early detection, prevalence monitoring, and variant identification of human diseases in large populations. However, using census-based population size instead of real-time population estimates can mislead the interpretation of data acquired from sewage, hindering assessment of representativeness, inference of prevalence, or comparisons of taxa across sites. Here, we show that taxon abundance and sub-species diversisty in gut-associated microbiomes are new feature space to utilize for human population estimation. Using a population-scale human gut microbiome sample of over 1,100 people, we found that taxon-abundance distributions of gut-associated multi-person microbiomes exhibited generalizable relationships with respect to human population size. Here and throughout this paper, the human population size is essentially the sample size from the wastewater sample. We present a new algorithm, MicrobiomeCensus, for estimating human population size from sewage samples. MicrobiomeCensus harnesses the inter-individual variability in human gut microbiomes and performs maximum likelihood estimation based on simultaneous deviation of multiple taxa's relative abundances from their population means. MicrobiomeCensus outperformed generic algorithms in data-driven simulation benchmarks and detected population size differences in field data. New theorems are provided to justify our approach. This research provides a mathematical framework for inferring population sizes in real time from sewage samples, paving the way for more accurate ecological and public health studies utilizing the sewage metagenome.

Update of medRxiv. 2020 Sep 23;:

Many neurons in the mammalian central nervous system have complex dendritic arborisations and active dendritic conductances that enable these cells to perform sophisticated computations. How dendritically targeted inhibition affects local dendritic excitability is not fully understood. Here we use computational models of branched dendrites to investigate where GABAergic synapses should be placed to minimise dendritic excitability over time. To do so, we formulate a metric we term the "Inhibitory Level" (IL), which quantifies the effectiveness of synaptic inhibition for reducing the depolarising effect of nearby excitatory input. GABAergic synaptic inhibition is dependent on the reversal potential for GABAA receptors (EGABA), which is primarily set by the transmembrane chloride ion (Cl-) concentration gradient. We, therefore, investigated how variable EGABA and dynamic chloride affects dendritic inhibition. We found that the inhibitory effectiveness of dendritic GABAergic synapses combines at an encircled branch junction. The extent of this inhibitory accumulation is dependent on the number of branches and location of synapses but is independent of EGABA. This inhibitory accumulation occurs even for very distally placed inhibitory synapses when they are hyperpolarising-but not when they are shunting. When accounting for Cl- fluxes and dynamics in Cl- concentration, we observed that Cl- loading is detrimental to inhibitory effectiveness. This enabled us to determine the most inhibitory distribution of GABAergic synapses which is close to-but not at-a shared branch junction. This distribution balances a trade-off between a stronger combined inhibitory influence when synapses closely encircle a branch junction with the deleterious effects of increased Cl- by loading that occurs when inhibitory synapses are co-located.

Detecting structures within the continuous diving behavior of marine animals is challenging, and no universal framework is available. We captured such diverse structures using chaos theory. By applying time-delay embedding to exceptionally long dive records (83 d) from the narwhal, we reconstructed the state-space portrait. Using measures of chaos, we detected a diurnal pattern and its seasonal modulation, classified data, and found how sea-ice appearance shifts time budgets. There is more near-surface rest but deeper dives at solar noon, and more intense diving during twilight and at night but to shallower depths (likely following squid); sea-ice appearance reduces rest. The introduced geometrical approach is simple to implement and potentially helpful for mapping and labeling long-term behavioral data, identifying differences between individual animals and species, and detecting perturbations.

Abnormal cell proliferation and epithelial-mesenchymal transition (EMT) are the essential events that induce cancer initiation and progression. A fundamental goal in cancer research is to develop an efficient method to detect mutational genes capable of driving cancer. Although several computational methods have been proposed to identify these key mutations, many of them focus on the association between genetic mutations and functional changes in relevant biological processes, but not their real causality. Causal effect inference provides a way to estimate the real induce effect of a certain mutation on vital biological processes of cancer initiation and progression, through addressing the confounder bias due to neutral mutations and unobserved latent variables. In this study, integrating genomic and transcriptomic data, we construct a novel causal inference model based on a deep variational autoencoder to identify key oncogenic somatic mutations. Applied to 10 cancer types, our method quantifies the causal effect of genetic mutations on cell proliferation and EMT by reducing both observed and unobserved confounding biases. The experimental results indicate that genes with higher mutation frequency do not necessarily mean they are more potent in inducing cancer and promoting cancer development. Moreover, our study fills a gap in the use of machine learning for causal inference to identify oncogenic mutations.

Since next-generation sequencing (NGS) has become widely available, large gene panels containing up to several hundred genes can be sequenced cost-efficiently. However, the interpretation of the often large numbers of sequence variants detected when using NGS is laborious, prone to errors and is often difficult to compare across laboratories. To overcome this challenge, the American College of Medical Genetics and Genomics and the Association for Molecular Pathology (ACMG/AMP) have introduced standards and guidelines for the interpretation of sequencing variants. Additionally, disease-specific refinements have been developed that include accurate thresholds for many criteria, enabling highly automated processing. This is of particular interest for common but heterogeneous disorders such as hearing impairment. With more than 200 genes associated with hearing disorders, the manual inspection of possible causative variants is particularly difficult and time-consuming. To this end, we developed the open-source bioinformatics tool GenOtoScope, which automates the analysis of all ACMG/AMP criteria that can be assessed without further individual patient information or human curator investigation, including the refined loss of function criterion ("PVS1"). Two types of interfaces are provided: (i) a command line application to classify sequence variants in batches for a set of patients and (ii) a user-friendly website to classify single variants. We compared the performance of our tool with two other variant classification tools using two hearing loss data sets, which were manually annotated either by the ClinGen Hearing Loss Gene Curation Expert Panel or the diagnostics unit of our human genetics department. GenOtoScope achieved the best average accuracy and precision for both data sets. Compared to the second-best tool, GenOtoScope improved the accuracy metric by 25.75% and 4.57% and precision metric by 52.11% and 12.13% on the two data sets, respectively. The web interface is accessible via: http://genotoscope.mh-hannover.de:5000 and the command line interface via: https://github.com/damianosmel/GenOtoScope.

The sequencing of antibody repertoires of B-cells at increasing coverage and depth has led to the identification of vast numbers of immunoglobulin heavy and light chains. However, the size and complexity of these Adaptive Immune Receptor Repertoire sequencing (AIRR-seq) datasets makes it difficult to perform exploratory analyses. To aid in data exploration, we have developed AIRRscape, an R Shiny-based interactive web browser application that enables B-cell receptor (BCR) and antibody feature discovery through comparisons among multiple repertoires. Using AIRR-seq data as input, AIRRscape starts by aggregating and sorting repertoires into interactive and explorable bins of germline V-gene, germline J-gene, and CDR3 length, providing a high-level view of the entire repertoire. Interesting subsets of repertoires can be quickly identified and selected, and then network topologies of CDR3 motifs can be generated for further exploration. Here we demonstrate AIRRscape using patient BCR repertoires and sequences of published monoclonal antibodies to investigate patterns of humoral immunity to three viral pathogens: SARS-CoV-2, HIV-1, and DENV (dengue virus). AIRRscape reveals convergent antibody sequences among datasets for all three pathogens, although HIV-1 antibody datasets display limited convergence and idiosyncratic responses. We have made AIRRscape available as a web-based Shiny application, along with code on GitHub to encourage its open development and use by immuno-informaticians, virologists, immunologists, vaccine developers, and other scientists that are interested in exploring and comparing multiple immune receptor repertoires.

In-vivo studies of pulmonary vascular disease and pulmonary hypertension (PH) have provided key insight into the progression of right ventricular (RV) dysfunction. Additional in-silico experiments using multiscale computational models have provided further details into biventricular mechanics and hemodynamic function in the presence of PH, yet few have assessed whether model parameters are practically identifiable prior to data collection. Moreover, none have used modeling to devise synergistic experimental designs. To address this knowledge gap, we conduct a practical identifiability analysis of a multiscale cardiovascular model across four simulated experimental designs. We determine a set of parameters using a combination of Morris screening and local sensitivity analysis, and test for practical identifiability using profile likelihood-based confidence intervals. We employ Markov chain Monte Carlo (MCMC) techniques to quantify parameter and model forecast uncertainty in the presence of noise corrupted data. Our results show that model calibration to only RV pressure suffers from practical identifiability issues and suffers from large forecast uncertainty in output space. In contrast, parameter and model forecast uncertainty is substantially reduced once additional left ventricular (LV) pressure and volume data is included. A comparison between single point systolic and diastolic LV data and continuous, time-dependent LV pressure-volume data reveals that at least some quantitative data from both ventricles should be included for future experimental studies.

With the recent approval by the FDA of the first disease-modifying drug for Alzheimer's Disease (AD), personalized medicine will be increasingly important for appropriate management and counseling of patients with AD and those at risk. The growing availability of clinical biomarker data and data-driven computational modeling techniques provide an opportunity for new approaches to individualized AD therapeutic planning. In this paper, we develop a new mathematical model, based on AD cognitive, cerebrospinal fluid (CSF) and MRI biomarkers, to provide a personalized optimal treatment plan for individuals. This model is parameterized by biomarker data from the AD Neuroimaging Initiative (ADNI) cohort, a large multi-institutional database monitoring the natural history of subjects with AD and mild cognitive impairment (MCI). Optimal control theory is used to incorporate time-varying treatment controls and side-effects into the model, based on recent clinical trial data, to provide a personalized treatment regimen with anti-amyloid-beta therapy. In-silico treatment studies were conducted on the approved treatment, aducanumab, as well as on another promising anti-amyloid-beta therapy under evaluation, donanemab. Clinical trial simulations were conducted over both short-term (78 weeks) and long-term (10 years) periods with low-dose (6 mg/kg) and high-dose (10 mg/kg) regimens for aducanumab, and a single-dose regimen (1400 mg) for donanemab. Results confirm those of actual clinical trials showing a large and sustained effect of both aducanumab and donanemab on amyloid beta clearance. The effect on slowing cognitive decline was modest for both treatments, but greater for donanemab. This optimal treatment computational modeling framework can be applied to other single and combination treatments for both prediction and optimization, as well as incorporate new clinical trial data as it becomes available.

The Polycomb system via the methylation of the lysine 27 of histone H3 (H3K27) plays central roles in the silencing of many lineage-specific genes during development. Recent experimental evidence suggested that the recruitment of histone modifying enzymes like the Polycomb repressive complex 2 (PRC2) at specific sites and their spreading capacities from these sites are key to the establishment and maintenance of a proper epigenomic landscape around Polycomb-target genes. Here, to test whether such mechanisms, as a minimal set of qualitative rules, are quantitatively compatible with data, we developed a mathematical model that can predict the locus-specific distributions of H3K27 modifications based on previous biochemical knowledge. Within the biological context of mouse embryonic stem cells, our model showed quantitative agreement with experimental profiles of H3K27 acetylation and methylation around Polycomb-target genes in wild-type and mutants. In particular, we demonstrated the key role of the reader-writer module of PRC2 and of the competition between the binding of activating and repressing enzymes in shaping the H3K27 landscape around transcriptional start sites. The predicted dynamics of establishment and maintenance of the repressive trimethylated H3K27 state suggest a slow accumulation, in perfect agreement with experiments. Our approach represents a first step towards a quantitative description of PcG regulation in various cellular contexts and provides a generic framework to better characterize epigenetic regulation in normal or disease situations.

The human brain generates a rich repertoire of spatio-temporal activity patterns, which support a wide variety of motor and cognitive functions. These patterns of activity change with age in a multi-factorial manner. One of these factors is the variations in the brain's connectomics that occurs along the lifespan. However, the precise relationship between high-order functional interactions and connnectomics, as well as their variations with age are largely unknown, in part due to the absence of mechanistic models that can efficiently map brain connnectomics to functional connectivity in aging. To investigate this issue, we have built a neurobiologically-realistic whole-brain computational model using both anatomical and functional MRI data from 161 participants ranging from 10 to 80 years old. We show that the differences in high-order functional interactions between age groups can be largely explained by variations in the connectome. Based on this finding, we propose a simple neurodegeneration model that is representative of normal physiological aging. As such, when applied to connectomes of young participant it reproduces the age-variations that occur in the high-order structure of the functional data. Overall, these results begin to disentangle the mechanisms by which structural changes in the connectome lead to functional differences in the ageing brain. Our model can also serve as a starting point for modeling more complex forms of pathological ageing or cognitive deficits.

The mapping from genotype to phenotype to fitness typically involves multiple nonlinearities that can transform the effects of mutations. For example, mutations may contribute additively to a phenotype, but their effects on fitness may combine non-additively because selection favors a low or intermediate value of that phenotype. This can cause incongruence between the topographical properties of a fitness landscape and its underlying genotype-phenotype landscape. Yet, genotype-phenotype landscapes are often used as a proxy for fitness landscapes to study the dynamics and predictability of evolution. Here, we use theoretical models and empirical data on transcription factor-DNA interactions to systematically study the incongruence of genotype-phenotype and fitness landscapes when selection favors a low or intermediate phenotypic value. Using the theoretical models, we prove a number of fundamental results. For example, selection for low or intermediate phenotypic values does not change simple sign epistasis into reciprocal sign epistasis, implying that genotype-phenotype landscapes with only simple sign epistasis motifs will always give rise to single-peaked fitness landscapes under such selection. More broadly, we show that such selection tends to create fitness landscapes that are more rugged than the underlying genotype-phenotype landscape, but this increased ruggedness typically does not frustrate adaptive evolution because the local adaptive peaks in the fitness landscape tend to be nearly as tall as the global peak. Many of these results carry forward to the empirical genotype-phenotype landscapes, which may help to explain why low- and intermediate-affinity transcription factor-DNA interactions are so prevalent in eukaryotic gene regulation.

Forecasts based on epidemiological modelling have played an important role in shaping public policy throughout the COVID-19 pandemic. This modelling combines knowledge about infectious disease dynamics with the subjective opinion of the researcher who develops and refines the model and often also adjusts model outputs. Developing a forecast model is difficult, resource- and time-consuming. It is therefore worth asking what modelling is able to add beyond the subjective opinion of the researcher alone. To investigate this, we analysed different real-time forecasts of cases of and deaths from COVID-19 in Germany and Poland over a 1-4 week horizon submitted to the German and Polish Forecast Hub. We compared crowd forecasts elicited from researchers and volunteers, against a) forecasts from two semi-mechanistic models based on common epidemiological assumptions and b) the ensemble of all other models submitted to the Forecast Hub. We found crowd forecasts, despite being overconfident, to outperform all other methods across all forecast horizons when forecasting cases (weighted interval score relative to the Hub ensemble 2 weeks ahead: 0.89). Forecasts based on computational models performed comparably better when predicting deaths (rel. WIS 1.26), suggesting that epidemiological modelling and human judgement can complement each other in important ways.

In Escherichia coli, translocation of RNA polymerase (RNAP) during transcription introduces supercoiling to DNA, which influences the initiation and elongation behaviors of RNAP. To quantify the role of supercoiling in transcription regulation, we developed a spatially resolved supercoiling model of transcription. The integrated model describes how RNAP activity feeds back with the local DNA supercoiling and how this mechanochemical feedback controls transcription, subject to topoisomerase activities and stochastic topological domain formation. This model establishes that transcription-induced supercoiling mediates the cooperation of co-transcribing RNAP molecules in highly expressed genes, and this cooperation is achieved under moderate supercoiling diffusion and high topoisomerase unbinding rates. It predicts that a topological domain could serve as a transcription regulator, generating substantial transcriptional noise. It also shows the relative orientation of two closely arranged genes plays an important role in regulating their transcription. The model provides a quantitative platform for investigating how genome organization impacts transcription.

Counting the number of species, items, or genes that are shared between two groups, sets, or communities is a simple calculation when sampling is complete. However, when only partial samples are available, quantifying the overlap between two communities becomes an estimation problem. Furthermore, to calculate normalized measures of β-diversity, such as the Jaccard and Sorenson-Dice indices, one must also estimate the total sizes of the communities being compared. Previous efforts to address these problems have assumed knowledge of total community sizes and then used Bayesian methods to produce unbiased estimates with quantified uncertainty. Here, we address communities of unknown size and show that this produces systematically better estimates-both in terms of central estimates and quantification of uncertainty in those estimates. We further show how to use species, item, or gene count data to refine estimates of community size in a Bayesian joint model of community size and overlap.

Understanding the segregation of cells is crucial to answer questions about tissue formation in embryos or tumor progression. Steinberg proposed that separation of cells can be compared to the separation of two liquids. Such a separation is well described by the Cahn-Hilliard (CH) equations and the segregation indices exhibit an algebraic decay with exponent 1/3 with respect to time. Similar exponents are also observed in cell-based models. However, the scaling behavior in these numerical models is usually only examined in the asymptotic regime and these models have not been directly applied to actual cell segregation data. In contrast, experimental data also reveals other scaling exponents and even slow logarithmic scaling laws. These discrepancies are commonly attributed to the effects of collective motion or velocity-dependent interactions. By calibrating a 2D cellular automaton (CA) model which efficiently implements a dynamic variant of the differential adhesion hypothesis to 2D experimental data from Méhes et al., we reproduce the biological cell segregation experiments with just adhesive forces. The segregation in the cellular automaton model follows a logarithmic scaling initially, which is in contrast to the proposed algebraic scaling with exponent 1/3. However, within the less than two orders of magnitudes in time which are observable in the experiments, a logarithmic scaling may appear as a pseudo-algebraic scaling. In particular, we demonstrate that the cellular automaton model can exhibit a range of exponents ≤1/3 for such a pseudo-algebraic scaling. Moreover, the time span of the experiment falls into the transitory regime of the cellular automaton rather than the asymptotic one. We additionally develop a method for the calibration of the 2D Cahn-Hilliard model and find a match with experimental data within the transitory regime of the Cahn-Hilliard model with exponent 1/4. On the one hand this demonstrates that the transitory behavior is relevant for the experiment rather than the asymptotic one. On the other hand this corroborates the ambiguity of the scaling behavior, when segregation processes can be only observed on short time spans.

Mathematical models of vector-borne infections, including malaria, often assume age-independent mortality rates of vectors, despite evidence that many insects senesce. In this study we present survival data on insecticide-resistant Anopheles gambiae s.l. from experiments in Côte d'Ivoire. We fit a constant mortality function and two age-dependent functions (logistic and Gompertz) to the data from mosquitoes exposed (treated) and not exposed (control) to insecticide-treated nets (ITNs), to establish biologically realistic survival functions. This enables us to explore the effects of insecticide exposure on mosquito mortality rates, and the extent to which insecticide resistance might impact the effectiveness of ITNs. We investigate this by calculating the expected number of infectious bites a mosquito will take in its lifetime, and by extension the vectorial capacity. Our results show that the predicted vectorial capacity is substantially lower in mosquitoes exposed to ITNs, despite the mosquitoes in the experiment being highly insecticide-resistant. The more realistic age-dependent functions provide a better fit to the experimental data compared to a constant mortality function and, hence, influence the predicted impact of ITNs on malaria transmission potential. In models with age-independent mortality, there is a great reduction for the vectorial capacity under exposure compared to no exposure. However, the two age-dependent functions predicted an even larger reduction due to exposure, highlighting the impact of incorporating age in the mortality rates. These results further show that multiple exposures to ITNs had a considerable effect on the vectorial capacity. Overall, the study highlights the importance of including age dependency in mathematical models of vector-borne disease transmission and in fully understanding the impact of interventions.

We introduce a novel, biologically plausible local learning rule that provably increases the robustness of neural dynamics to noise in nonlinear recurrent neural networks with homogeneous nonlinearities. Our learning rule achieves higher noise robustness without sacrificing performance on the task and without requiring any knowledge of the particular task. The plasticity dynamics-an integrable dynamical system operating on the weights of the network-maintains a multiplicity of conserved quantities, most notably the network's entire temporal map of input to output trajectories. The outcome of our learning rule is a synaptic balancing between the incoming and outgoing synapses of every neuron. This synaptic balancing rule is consistent with many known aspects of experimentally observed heterosynaptic plasticity, and moreover makes new experimentally testable predictions relating plasticity at the incoming and outgoing synapses of individual neurons. Overall, this work provides a novel, practical local learning rule that exactly preserves overall network function and, in doing so, provides new conceptual bridges between the disparate worlds of the neurobiology of heterosynaptic plasticity, the engineering of regularized noise-robust networks, and the mathematics of integrable Lax dynamical systems.

The mutual interactions between the synthetic gene circuits and the host growth could cause unexpected outcomes in the dynamical behaviors of the circuits. However, how the steady states and the stabilities of the gene circuits are affected by host cell growth is not fully understood. Here, we developed a mathematical model for nonlinear growth feedback based on published experimental data. The model analysis predicts that growth feedback could significantly change the qualitative states of the system. Bistability could emerge in a circuit without positive feedback, and high-order multistability (three or more steady states) arises in the self-activation and toggle switch circuits. Our results provide insight into the potential effects of ultrasensitive growth feedback on the emergence of qualitative states in synthetic circuits and the corresponding underlying mechanism.

Unraveling the complex demographic histories of natural populations is a central problem in population genetics. Understanding past demographic events is of general anthropological interest, but is also an important step in establishing accurate null models when identifying adaptive or disease-associated genetic variation. An important class of tools for inferring past population size changes from genomic sequence data are Coalescent Hidden Markov Models (CHMMs). These models make efficient use of the linkage information in population genomic datasets by using the local genealogies relating sampled individuals as latent states that evolve along the chromosome in an HMM framework. Extending these models to large sample sizes is challenging, since the number of possible latent states increases rapidly. Here, we present our method CHIMP (CHMM History-Inference Maximum-Likelihood Procedure), a novel CHMM method for inferring the size history of a population. It can be applied to large samples (hundreds of haplotypes) and only requires unphased genomes as input. The two implementations of CHIMP that we present here use either the height of the genealogical tree (TMRCA) or the total branch length, respectively, as the latent variable at each position in the genome. The requisite transition and emission probabilities are obtained by numerically solving certain systems of differential equations derived from the ancestral process with recombination. The parameters of the population size history are subsequently inferred using an Expectation-Maximization algorithm. In addition, we implement a composite likelihood scheme to allow the method to scale to large sample sizes. We demonstrate the efficiency and accuracy of our method in a variety of benchmark tests using simulated data and present comparisons to other state-of-the-art methods. Specifically, our implementation using TMRCA as the latent variable shows comparable performance and provides accurate estimates of effective population sizes in intermediate and ancient times. Our method is agnostic to the phasing of the data, which makes it a promising alternative in scenarios where high quality data is not available, and has potential applications for pseudo-haploid data.

Despite the immense progress recently witnessed in protein structure prediction, the modeling accuracy for proteins that lack sequence and/or structure homologs remains to be improved. We developed an open-source program, DeepFold, which integrates spatial restraints predicted by multi-task deep residual neural-networks along with a knowledge-based energy function to guide its gradient-descent folding simulations. The results on large-scale benchmark tests showed that DeepFold creates full-length models with accuracy significantly beyond classical folding approaches and other leading deep learning methods. Of particular interest is the modeling performance on the most difficult targets with very few homologous sequences, where DeepFold achieved an average TM-score that was 40.3% higher than trRosetta and 44.9% higher than DMPfold. Furthermore, the folding simulations for DeepFold were 262 times faster than traditional fragment assembly simulations. These results demonstrate the power of accurately predicted deep learning potentials to improve both the accuracy and speed of ab initio protein structure prediction.

Extracting information on the selective and demographic past of populations that is contained in samples of genome sequences requires a description of the distribution of the underlying genealogies. Using the Laplace transform, this distribution can be generated with a simple recursive procedure, regardless of model complexity. Assuming an infinite-sites mutation model, the probability of observing specific configurations of linked variants within small haplotype blocks can be recovered from the Laplace transform of the joint distribution of branch lengths. However, the repeated differentiation required to compute these probabilities has proven to be a serious computational bottleneck in earlier implementations. Here, I show that the state space diagram can be turned into a computational graph, allowing efficient evaluation of the Laplace transform by means of a graph traversal algorithm. This general algorithm can, for example, be applied to tabulate the likelihoods of mutational configurations in non-recombining blocks. This work provides a crucial speed up for existing composite likelihood approaches that rely on the joint distribution of branch lengths to fit isolation with migration models and estimate the parameters of selective sweeps. The associated software is available as an open-source Python library, agemo.

Adult hippocampal neurogenesis and hippocampus-dependent cognition in mice have been found to be adversely affected by hypomagnetic field exposure. The effect concurred with a reduction of reactive oxygen species in the absence of the geomagnetic field. A recent theoretical study suggests a mechanistic interpretation of this phenomenon in the framework of the Radical Pair Mechanism. According to this model, a flavin-superoxide radical pair, born in the singlet spin configuration, undergoes magnetic field-dependent spin dynamics such that the pair's recombination is enhanced as the applied magnetic field is reduced. This model has two ostensible weaknesses: a) the assumption of a singlet initial state is irreconcilable with known reaction pathways generating such radical pairs, and b) the model neglects the swift spin relaxation of free superoxide, which abolishes any magnetic sensitivity in geomagnetic/hypomagnetic fields. We here suggest that a model based on a radical triad and the assumption of a secondary radical scavenging reaction can, in principle, explain the phenomenon without unnatural assumptions, thus providing a coherent explanation of hypomagnetic field effects in biology.

Functional forms of biophysically-realistic neuron models are constrained by neurobiological and anatomical considerations, such as cell morphologies and the presence of known ion channels. Despite these constraints, neuron models still contain unknown static parameters which must be inferred from experiment. This inference task is most readily cast into the framework of state-space models, which systematically takes into account partial observability and measurement noise. Inferring only dynamical state variables such as membrane voltages is a well-studied problem, and has been approached with a wide range of techniques beginning with the well-known Kalman filter. Inferring both states and fixed parameters, on the other hand, is less straightforward. Here, we develop a method for joint parameter and state inference that combines traditional state space modeling with chaotic synchronization and optimal control. Our methods are tailored particularly to situations with considerable measurement noise, sparse observability, very nonlinear or chaotic dynamics, and highly uninformed priors. We illustrate our approach both in a canonical chaotic model and in a phenomenological neuron model, showing that many unknown parameters can be uncovered reliably and accurately from short and noisy observed time traces. Our method holds promise for estimation in larger-scale systems, given ongoing improvements in calcium reporters and genetically-encoded voltage indicators.

Thrombin is an enzyme produced during blood coagulation that is crucial to the formation of a stable clot. Thrombin cleaves soluble fibrinogen into fibrin, which polymerizes and forms an insoluble, stabilizing gel around the growing clot. A small fraction of circulating fibrinogen is the variant γA/γ', which has been associated with high-affinity thrombin binding and implicated as a risk factor for myocardial infarctions, deep vein thrombosis, and coronary artery disease. Thrombin is also known to be strongly sequestered by polymerized fibrin for extended periods of time in a way that is partially regulated by γA/γ'. However, the role of γA/γ'-thrombin interactions during fibrin polymerization is not fully understood. Here, we present a mathematical model of fibrin polymerization that considered the interactions between thrombin, fibrinogen, and fibrin, including those with γA/γ'. In our model, bivalent thrombin-fibrin binding greatly increased thrombin residency times and allowed for thrombin-trapping during fibrin polymerization. Results from the model showed that early in fibrin polymerization, γ' binding to thrombin served to localize the thrombin to the fibrin(ogen), which effectively enhanced the enzymatic conversion of fibrinogen to fibrin. When all the fibrin was fully generated, however, the fibrin-thrombin binding persisted but the effect of fibrin on thrombin switched quickly to serve as a sink, essentially removing all free thrombin from the system. This dual role for γ'-thrombin binding during polymerization led to a paradoxical decrease in trapped thrombin as the amount of γ' was increased. The model highlighted biochemical and biophysical roles for fibrin-thrombin interactions during polymerization and agreed well with experimental observations.

The ubiquitous use of computational work for data generation, processing, and modeling increased the importance of digital documentation in improving research quality and impact. Computational notebooks are files that contain descriptive text, as well as code and its outputs, in a single, dynamic, and visually appealing file that is easier to understand by nonspecialists. Traditionally used by data scientists when producing reports and informing decision-making, the use of this tool in research publication is not common, despite its potential to increase research impact and quality. For a single study, the content of such documentation partially overlaps with that of classical lab notebooks and that of the scientific manuscript reporting the study. Therefore, to minimize the amount of work required to manage all the files related to these contents and optimize their production, we present a starter kit to facilitate the implementation of computational notebooks in the research process, including publication. The kit contains the template of a computational notebook integrated into a research project that employs R, Python, or Julia. Using examples of ecological studies, we show how computational notebooks also foster the implementation of principles of Open Science, such as reproducibility and traceability. The kit is designed for beginners, but at the end we present practices that can be gradually implemented to develop a fully digital research workflow. Our hope is that such minimalist yet effective starter kit will encourage researchers to adopt this practice in their workflow, regardless of their computational background.

Honeybees protect their colony against vertebrates by mass stinging and they coordinate their actions during this crucial event thanks to an alarm pheromone carried directly on the stinger, which is therefore released upon stinging. The pheromone then recruits nearby bees so that more and more bees participate in the defence. However, a quantitative understanding of how an individual bee adapts its stinging response during the course of an attack is still a challenge: Typically, only the group behaviour is effectively measurable in experiment; Further, linking the observed group behaviour with individual responses requires a probabilistic model enumerating a combinatorial number of possible group contexts during the defence; Finally, extracting the individual characteristics from group observations requires novel methods for parameter inference. We first experimentally observed the behaviour of groups of bees confronted with a fake predator inside an arena and quantified their defensive reaction by counting the number of stingers embedded in the dummy at the end of a trial. We propose a biologically plausible model of this phenomenon, which transparently links the choice of each individual bee to sting or not, to its group context at the time of the decision. Then, we propose an efficient method for inferring the parameters of the model from the experimental data. Finally, we use this methodology to investigate the effect of group size on stinging initiation and alarm pheromone recruitment. Our findings shed light on how the social context influences stinging behaviour, by quantifying how the alarm pheromone concentration level affects the decision of each bee to sting or not in a given group size. We show that recruitment is curbed as group size grows, thus suggesting that the presence of nestmates is integrated as a negative cue by individual bees. Moreover, the unique integration of exact and statistical methods provides a quantitative characterisation of uncertainty associated to each of the inferred parameters.

The rodent vibrissal (whisker) system has been studied for decades as a model of active touch sensing. There are no sensors along the length of a whisker; all sensing occurs at the whisker base. Therefore, a large open question in many neuroscience studies is how an animal could estimate the three-dimensional (3D) location at which a whisker makes contact with an object. In the present work we simulated the shape of a real rat whisker to demonstrate the existence of several unique mappings from triplets of mechanical signals at the whisker base to the three-dimensional whisker-object contact point. We then used high speed video to record whisker deflections as an awake rat whisked against a peg, and used the mechanics resulting from those deflections to extract the contact points along the peg surface. These results demonstrate that measurement of specific mechanical triplets at the base of a biological whisker can enable 3D contact point determination during natural whisking behavior. The approach is viable even though the biological whisker has non-ideal, non-planar curvature, and even given the rat's real-world choices of whisking parameters. Visual intuition for the quality of the approach is provided in a video that shows the contour of the peg gradually emerging during active whisking behavior.

Magnesium (Mg2+) ions are frequently reported to regulate osteogenic activities of mesenchymal stem cells (MSCs). In this study, we propose a numerical model to study the regulatory importance of Mg2+ ions on MSCs osteoblastic differentiation in the presence of an inflammatory response. A fuzzy logic controller was formulated to receive the concentrations of Mg2+ ions and the inflammatory cytokines of TNF-α, IL-10, IL-1β, and IL-8 as cellular inputs and predict the cells' early and late differentiation rates. Five sets of empirical data obtained from published cell culture experiments were used to calibrate the model. The model successfully reproduced the empirical data regarding the concentration- and phase-dependent effect of Mg2+ ions on the differentiation process. In agreement with the experiments, the model showed the stimulatory role of Mg2+ ions on the early differentiation phase, once administered at low concentration, and their inhibitory role on the late differentiation phase. The numerical approach used in this study suggested 6-8 mM as the most effective concentration of Mg2+ ions in promoting the early differentiation process. Also, the proposed model sheds light on the fundamental differences in the behavioral properties of cells cultured in different experiments, e.g. differentiation rate and the sensitivity of the cultured cells to stimulatory signals such as Mg2+ ions. Thus, it can be used to interpret and compare different empirical findings. Moreover, the model successfully reproduced the nonlinearities in the concentration-dependent role of the inflammatory cytokines in early and late differentiation rates. Overall, the proposed model can be employed in studying the osteogenic properties of Mg-based implants in the presence of an inflammatory response.

In order to accurately parse the visual scene into distinct surfaces, it is essential to determine whether a local luminance edge is caused by a boundary between two surfaces or a shadow cast across a single surface. Previous studies have demonstrated that local chromatic cues may help to distinguish edges caused by shadows from those caused by surface boundaries, but the information potentially available in local achromatic cues like contrast, texture, and penumbral blur remains poorly understood. In this study, we develop and analyze a large database of hand-labeled achromatic shadow edges to better understand what image properties distinguish them from occlusion edges. We find that both the highest contrast as well as the lowest contrast edges are more likely to be occlusions than shadows, extending previous observations based on a more limited image set. We also find that contrast cues alone can reliably distinguish the two edge categories with nearly 70% accuracy at 40x40 resolution. Logistic regression on a Gabor Filter bank (GFB) modeling a population of V1 simple cells separates the categories with nearly 80% accuracy, and furthermore exhibits tuning to penumbral blur. A Filter-Rectify Filter (FRF) style neural network extending the GFB model performed at better than 80% accuracy, and exhibited blur tuning and greater sensitivity to texture differences. We compare human performance on our edge classification task to that of the FRF and GFB models, finding the best human observers attaining the same performance as the machine classifiers. Several analyses demonstrate both classifiers exhibit significant positive correlation with human behavior, although we find a slightly better agreement on an image-by-image basis between human performance and the FRF model than the GFB model, suggesting an important role for texture.

Accurately predicting contact between our bodies and environmental objects is paramount to our evolutionary survival. It has been hypothesized that multisensory neurons responding both to touch on the body, and to auditory or visual stimuli occurring near them-thus delineating our peripersonal space (PPS)-may be a critical player in this computation. However, we lack a normative account (i.e., a model specifying how we ought to compute) linking impact prediction and PPS encoding. Here, we leverage Bayesian Decision Theory to develop such a model and show that it recapitulates many of the characteristics of PPS. Namely, a normative model of impact prediction (i) delineates a graded boundary between near and far space, (ii) demonstrates an enlargement of PPS as the speed of incoming stimuli increases, (iii) shows stronger contact prediction for looming than receding stimuli-but critically is still present for receding stimuli when observation uncertainty is non-zero-, (iv) scales with the value we attribute to environmental objects, and finally (v) can account for the differing sizes of PPS for different body parts. Together, these modeling results support the conjecture that PPS reflects the computation of impact prediction, and make a number of testable predictions for future empirical studies.

BACKGROUND: Finding microbiome associations with possibly censored survival times is an important problem, especially as specific taxa could serve as biomarkers for disease prognosis or as targets for therapeutic interventions. The two existing methods for survival outcomes, MiRKAT-S and OMiSA, are restricted to testing associations at the community level and do not provide results at the individual taxon level. An ad hoc approach testing each taxon with a survival outcome using the Cox proportional hazard model may not perform well in the microbiome setting with sparse count data and small sample sizes. METHODS: We have previously developed the linear decomposition model (LDM) for testing continuous or discrete outcomes that unifies community-level and taxon-level tests into one framework. Here we extend the LDM to test survival outcomes. We propose to use the Martingale residuals or the deviance residuals obtained from the Cox model as continuous covariates in the LDM. We further construct tests that combine the results of analyzing each set of residuals separately. Finally, we extend PERMANOVA, the most commonly used distance-based method for testing community-level hypotheses, to handle survival outcomes in a similar manner. RESULTS: Using simulated data, we showed that the LDM-based tests preserved the false discovery rate for testing individual taxa and had good sensitivity. The LDM-based community-level tests and PERMANOVA-based tests had comparable or better power than MiRKAT-S and OMiSA. An analysis of data on the association of the gut microbiome and the time to acute graft-versus-host disease revealed several dozen associated taxa that would not have been achievable by any community-level test, as well as improved community-level tests by the LDM and PERMANOVA over those obtained using MiRKAT-S and OMiSA. CONCLUSIONS: Unlike existing methods, our new methods are capable of discovering individual taxa that are associated with survival times, which could be of important use in clinical settings.

A growing body of evidence highlights the intricate linkage of exteroceptive perception to the rhythmic activity of the visceral body. In parallel, interoceptive inference theories of affective perception and self-consciousness are on the rise in cognitive science. However, thus far no formal theory has emerged to integrate these twin domains; instead, most extant work is conceptual in nature. Here, we introduce a formal model of cardiac active inference, which explains how ascending cardiac signals entrain exteroceptive sensory perception and uncertainty. Through simulated psychophysics, we reproduce the defensive startle reflex and commonly reported effects linking the cardiac cycle to affective behaviour. We further show that simulated 'interoceptive lesions' blunt affective expectations, induce psychosomatic hallucinations, and exacerbate biases in perceptual uncertainty. Through synthetic heart-rate variability analyses, we illustrate how the balance of arousal-priors and visceral prediction errors produces idiosyncratic patterns of physiological reactivity. Our model thus offers a roadmap for computationally phenotyping disordered brain-body interaction.

The over-expression of the Bcl-2 protein is a common feature of many solid cancers and hematological malignancies, and it is typically associated with poor prognosis and resistance to chemotherapy. Bcl-2-specific inhibitors, such as venetoclax, have recently been approved for the treatment of chronic lymphocytic leukemia and small lymphocytic lymphoma, and they are showing promise in clinical trials as a targeted therapy for patients with relapsed or refractory acute myeloid leukemia (AML). However, successful treatment of AML with Bcl-2-specific inhibitors is often followed by the rapid development of drug resistance. An emerging paradigm for overcoming drug resistance in cancer treatment is through the targeting of mitochondrial energetics and metabolism. In AML in particular, it was recently observed that inhibition of mitochondrial translation via administration of the antibiotic tedizolid significantly affects mitochondrial bioenergetics, activating the integrated stress response (ISR) and subsequently sensitizing drug-resistant AML cells to venetoclax. Here we develop an integrative systems biology approach to acquire a deeper understanding of the molecular mechanisms behind this process, and in particular, of the specific role of the ISR in the commitment of cells to apoptosis. Our multi-scale mathematical model couples the ISR to the intrinsic apoptosis pathway in venetoclax-resistant AML cells, includes the metabolic effects of treatment, and integrates RNA, protein level, and cellular viability data. Using the mathematical model, we identify the dominant mechanisms by which ISR activation helps to overcome venetoclax resistance, and we study the temporal sequencing of combination treatment to determine the most efficient and robust combination treatment protocol.

Dendrites of cortical pyramidal cells are densely populated by hyperpolarization-activated cyclic nucleotide-gated (HCN) channels, a.k.a. Ih channels. Ih channels are targeted by multiple neuromodulatory pathways, and thus are one of the key ion-channel populations regulating the pyramidal cell activity. Previous observations and theories attribute opposing effects of the Ih channels on neuronal excitability due to their mildly hyperpolarized reversal potential. These effects are difficult to measure experimentally due to the fine spatiotemporal landscape of the Ih activity in the dendrites, but computational models provide an efficient tool for studying this question in a reduced but generalizable setting. In this work, we build upon existing biophysically detailed models of thick-tufted layer V pyramidal cells and model the effects of over- and under-expression of Ih channels as well as their neuromodulation. We show that Ih channels facilitate the action potentials of layer V pyramidal cells in response to proximal dendritic stimulus while they hinder the action potentials in response to distal dendritic stimulus at the apical dendrite. We also show that the inhibitory action of the Ih channels in layer V pyramidal cells is due to the interactions between Ih channels and a hot zone of low voltage-activated Ca2+ channels at the apical dendrite. Our simulations suggest that a combination of Ih-enhancing neuromodulation at the proximal part of the apical dendrite and Ih-inhibiting modulation at the distal part of the apical dendrite can increase the layer V pyramidal excitability more than either of the two alone. Our analyses uncover the effects of Ih-channel neuromodulation of layer V pyramidal cells at a single-cell level and shed light on how these neurons integrate information and enable higher-order functions of the brain.

We perform a thorough analysis of RNA velocity methods, with a view towards understanding the suitability of the various assumptions underlying popular implementations. In addition to providing a self-contained exposition of the underlying mathematics, we undertake simulations and perform controlled experiments on biological datasets to assess workflow sensitivity to parameter choices and underlying biology. Finally, we argue for a more rigorous approach to RNA velocity, and present a framework for Markovian analysis that points to directions for improvement and mitigation of current problems.

Interactions across frontal cortex are critical for cognition. Animal studies suggest a role for mediodorsal thalamus (MD) in these interactions, but the computations performed and direct relevance to human decision making are unclear. Here, inspired by animal work, we extended a neural model of an executive frontal-MD network and trained it on a human decision-making task for which neuroimaging data were collected. Using a biologically-plausible learning rule, we found that the model MD thalamus compressed its cortical inputs (dorsolateral prefrontal cortex, dlPFC) underlying stimulus-response representations. Through direct feedback to dlPFC, this thalamic operation efficiently partitioned cortical activity patterns and enhanced task switching across different contingencies. To account for interactions with other frontal regions, we expanded the model to compute higher-order strategy signals outside dlPFC, and found that the MD offered a more efficient route for such signals to switch dlPFC activity patterns. Human fMRI data provided evidence that the MD engaged in feedback to dlPFC, and had a role in routing orbitofrontal cortex inputs when subjects switched behavioral strategy. Collectively, our findings contribute to the emerging evidence for thalamic regulation of frontal interactions in the human brain.

Studies comparing single cell RNA-Seq (scRNA-Seq) data between conditions mainly focus on differences in the proportion of cell types or on differentially expressed genes. In many cases these differences are driven by changes in cell interactions which are challenging to infer without spatial information. To determine cell-cell interactions that differ between conditions we developed the Cell Interaction Network Inference (CINS) pipeline. CINS combines Bayesian network analysis with regression-based modeling to identify differential cell type interactions and the proteins that underlie them. We tested CINS on a disease case control and on an aging mouse dataset. In both cases CINS correctly identifies cell type interactions and the ligands involved in these interactions improving on prior methods suggested for cell interaction predictions. We performed additional mouse aging scRNA-Seq experiments which further support the interactions identified by CINS.

Today, there is great interest in diets proposing new combinations of macronutrient compositions and fasting schedules. Unfortunately, there is little consensus regarding the impact of these different diets, since available studies measure different sets of variables in different populations, thus only providing partial, non-connected insights. We lack an approach for integrating all such partial insights into a useful and interconnected big picture. Herein, we present such an integrating tool. The tool uses a novel mathematical model that describes mechanisms regulating diet response and fasting metabolic fluxes, both for organ-organ crosstalk, and inside the liver. The tool can mechanistically explain and integrate data from several clinical studies, and correctly predict new independent data, including data from a new study. Using this model, we can predict non-measured variables, e.g. hepatic glycogen and gluconeogenesis, in response to fasting and different diets. Furthermore, we exemplify how such metabolic responses can be successfully adapted to a specific individual's sex, weight, height, as well as to the individual's historical data on metabolite dynamics. This tool enables an offline digital twin technology.

The great advances of sequencing technologies allow the in vivo measurement of nuclear processes-such as DNA repair after UV exposure-over entire cell populations. However, data sets usually contain only a few samples over several hours, missing possibly important information in between time points. We developed a data-driven approach to analyse CPD repair kinetics over time in Saccharomyces cerevisiae. In contrast to other studies that consider sequencing signals as an average behaviour, we understand them as the superposition of signals from independent cells. By motivating repair as a stochastic process, we derive a minimal model for which the parameters can be conveniently estimated. We correlate repair parameters to a variety of genomic features that are assumed to influence repair, including transcription rate and nucleosome density. The clearest link was found for the transcription unit length, which has been unreported for budding yeast to our knowledge. The framework hence allows a comprehensive analysis of nuclear processes on a population scale.

Prediction of therapeutic peptide is a significant step for the discovery of promising therapeutic drugs. Most of the existing studies have focused on the mono-functional therapeutic peptide prediction. However, the number of multi-functional therapeutic peptides (MFTP) is growing rapidly, which requires new computational schemes to be proposed to facilitate MFTP discovery. In this study, based on multi-head self-attention mechanism and class weight optimization algorithm, we propose a novel model called PrMFTP for MFTP prediction. PrMFTP exploits multi-scale convolutional neural network, bi-directional long short-term memory, and multi-head self-attention mechanisms to fully extract and learn informative features of peptide sequence to predict MFTP. In addition, we design a class weight optimization scheme to address the problem of label imbalanced data. Comprehensive evaluation demonstrate that PrMFTP is superior to other state-of-the-art computational methods for predicting MFTP. We provide a user-friendly web server of PrMFTP, which is available at http://bioinfo.ahu.edu.cn/PrMFTP.

Many genes have nycthemeral rhythms of expression, i.e. a 24-hours periodic variation, at either mRNA or protein level or both, and most rhythmic genes are tissue-specific. Here, we investigate and discuss the evolutionary origins of rhythms in gene expression. Our results suggest that rhythmicity of protein expression could have been favored by selection to minimize costs. Trends are consistent in bacteria, plants and animals, and are also supported by tissue-specific patterns in mouse. Unlike for protein level, cost cannot explain rhythm at the RNA level. We suggest that instead it allows to periodically reduce expression noise. Noise control had the strongest support in mouse, with limited evidence in other species. We have also found that genes under stronger purifying selection are rhythmically expressed at the mRNA level, and we propose that this is because they are noise sensitive genes. Finally, the adaptive role of rhythmic expression is supported by rhythmic genes being highly expressed yet tissue-specific. This provides a good evolutionary explanation for the observation that nycthemeral rhythms are often tissue-specific.

Determining transcriptional factor binding sites (TFBSs) is critical for understanding the molecular mechanisms regulating gene expression in different biological conditions. Biological assays designed to directly mapping TFBSs require large sample size and intensive resources. As an alternative, ATAC-seq assay is simple to conduct and provides genomic cleavage profiles that contain rich information for imputing TFBSs indirectly. Previous footprint-based tools are inheritably limited by the accuracy of their bias correction algorithms and the efficiency of their feature extraction models. Here we introduce TAMC (Transcriptional factor binding prediction from ATAC-seq profile at Motif-predicted binding sites using Convolutional neural networks), a deep-learning approach for predicting motif-centric TF binding activity from paired-end ATAC-seq data. TAMC does not require bias correction during signal processing. By leveraging a one-dimensional convolutional neural network (1D-CNN) model, TAMC make predictions based on both footprint and non-footprint features at binding sites for each TF and outperforms existing footprinting tools in TFBS prediction particularly for ATAC-seq data with limited sequencing depth.

The reproductive number is an important metric that has been widely used to quantify the infectiousness of communicable diseases. The time-varying instantaneous reproductive number is useful for monitoring the real-time dynamics of a disease to inform policy making for disease control. Local estimation of this metric, for instance at a county or city level, allows for more targeted interventions to curb transmission. However, simultaneous estimation of local reproductive numbers must account for potential sources of heterogeneity in these time-varying quantities-a key element of which is human mobility. We develop a statistical method that incorporates human mobility between multiple regions for estimating region-specific instantaneous reproductive numbers. The model also can account for exogenous cases imported from outside of the regions of interest. We propose two approaches to estimate the reproductive numbers, with mobility data used to adjust incidence in the first approach and to inform a formal priori distribution in the second (Bayesian) approach. Through a simulation study, we show that region-specific reproductive numbers can be well estimated if human mobility is reasonably well approximated by available data. We use this approach to estimate the instantaneous reproductive numbers of COVID-19 for 14 counties in Massachusetts using CDC case report data and the human mobility data collected by SafeGraph. We found that, accounting for mobility, our method produces estimates of reproductive numbers that are distinct across counties. In contrast, independent estimation of county-level reproductive numbers tends to produce similar values, as trends in county case-counts for the state are fairly concordant. These approaches can also be used to estimate any heterogeneity in transmission, for instance, age-dependent instantaneous reproductive number estimates. As people are more mobile and interact frequently in ways that permit transmission, it is important to account for this in the estimation of the reproductive number.