Clinical trial data-sharing is seen as an imperative for research integrity and is becoming increasingly encouraged or even required by funders, journals, and other stakeholders. However, early experiences with data-sharing have been disappointing because they are not always conducted properly. Health data is indeed sensitive and not always easy to share in a responsible way. We propose 10 rules for researchers wishing to share their data. These rules cover the majority of elements to be considered in order to start the commendable process of clinical trial data-sharing: Rule 1: Abide by local legal and regulatory data protection requirementsRule 2: Anticipate the possibility of clinical trial data-sharing before obtaining fundingRule 3: Declare your intent to share data in the registration stepRule 4: Involve research participantsRule 5: Determine the method of data accessRule 6: Remember there are several other elements to shareRule 7: Do not proceed aloneRule 8: Deploy optimal data management to ensure that the data shared is usefulRule 9: Minimize risksRule 10: Strive for excellence.

Injuries to the skin heal through coordinated action of fibroblast-mediated extracellular matrix (ECM) deposition, ECM remodeling, and wound contraction. Defects involving the dermis result in fibrotic scars featuring increased stiffness and altered collagen content and organization. Although computational models are crucial to unravel the underlying biochemical and biophysical mechanisms, simulations of the evolving wound biomechanics are seldom benchmarked against measurements. Here, we leverage recent quantifications of local tissue stiffness in murine wounds to refine a previously-proposed systems-mechanobiological finite-element model. Fibroblasts are considered as the main cell type involved in ECM remodeling and wound contraction. Tissue rebuilding is coordinated by the release and diffusion of a cytokine wave, e.g. TGF-β, itself developed in response to an earlier inflammatory signal triggered by platelet aggregation. We calibrate a model of the evolving wound biomechanics through a custom-developed hierarchical Bayesian inverse analysis procedure. Further calibration is based on published biochemical and morphological murine wound healing data over a 21-day healing period. The calibrated model recapitulates the temporal evolution of: inflammatory signal, fibroblast infiltration, collagen buildup, and wound contraction. Moreover, it enables in silico hypothesis testing, which we explore by: (i) quantifying the alteration of wound contraction profiles corresponding to the measured variability in local wound stiffness; (ii) proposing alternative constitutive links connecting the dynamics of the biochemical fields to the evolving mechanical properties; (iii) discussing the plausibility of a stretch- vs. stiffness-mediated mechanobiological coupling. Ultimately, our model challenges the current understanding of wound biomechanics and mechanobiology, beside offering a versatile tool to explore and eventually control scar fibrosis after injury.

The majority of gene expression studies focus on the search for genes whose mean expression is different between two or more populations of samples in the so-called "differential expression analysis" approach. However, a difference in variance in gene expression may also be biologically and physiologically relevant. In the classical statistical model used to analyze RNA-sequencing (RNA-seq) data, the dispersion, which defines the variance, is only considered as a parameter to be estimated prior to identifying a difference in mean expression between conditions of interest. Here, we propose to evaluate four recently published methods, which detect differences in both the mean and dispersion in RNA-seq data. We thoroughly investigated the performance of these methods on simulated datasets and characterized parameter settings to reliably detect genes with a differential expression dispersion. We applied these methods to The Cancer Genome Atlas datasets. Interestingly, among the genes with an increased expression dispersion in tumors and without a change in mean expression, we identified some key cellular functions, most of which were related to catabolism and were overrepresented in most of the analyzed cancers. In particular, our results highlight autophagy, whose role in cancerogenesis is context-dependent, illustrating the potential of the differential dispersion approach to gain new insights into biological processes and to discover new biomarkers.

[This corrects the article DOI: 10.1371/journal.pcbi.1009492.].

When exposed to increased mechanical resistance from the soil, plant roots display non-linear growth responses that cannot be solely explained by mechanical principles. Here, we aim to investigate how changes in tissue mechanical properties are biologically regulated in response to soil strength. A particle-based model was developed to solve root-soil mechanical interactions at the cellular scale, and a detailed numerical study explored factors that affect root responses to soil resistance. Results showed how softening of root tissues at the tip may contribute to root responses to soil impedance, a mechanism likely linked to soil cavity expansion. The model also predicted the shortening and decreased anisotropy of the zone where growth occurs, which may improve the mechanical stability of the root against axial forces. The study demonstrates the potential of advanced modeling tools to help identify traits that confer plant resistance to abiotic stress.

[This corrects the article DOI: 10.1371/journal.pcbi.1010220.].

The availability of patient cohorts with several types of omics data opens new perspectives for exploring the disease's underlying biological processes and developing predictive models. It also comes with new challenges in computational biology in terms of integrating high-dimensional and heterogeneous data in a fashion that captures the interrelationships between multiple genes and their functions. Deep learning methods offer promising perspectives for integrating multi-omics data. In this paper, we review the existing integration strategies based on autoencoders and propose a new customizable one whose principle relies on a two-phase approach. In the first phase, we adapt the training to each data source independently before learning cross-modality interactions in the second phase. By taking into account each source's singularity, we show that this approach succeeds at taking advantage of all the sources more efficiently than other strategies. Moreover, by adapting our architecture to the computation of Shapley additive explanations, our model can provide interpretable results in a multi-source setting. Using multiple omics sources from different TCGA cohorts, we demonstrate the performance of the proposed method for cancer on test cases for several tasks, such as the classification of tumor types and breast cancer subtypes, as well as survival outcome prediction. We show through our experiments the great performances of our architecture on seven different datasets with various sizes and provide some interpretations of the results obtained. Our code is available on (https://github.com/HakimBenkirane/CustOmics).

Missense variants in genes encoding ion channels are associated with a spectrum of severe diseases. Variant effects on biophysical function correlate with clinical features and can be categorized as gain- or loss-of-function. This information enables a timely diagnosis, facilitates precision therapy, and guides prognosis. Functional characterization presents a bottleneck in translational medicine. Machine learning models may be able to rapidly generate supporting evidence by predicting variant functional effects. Here, we describe a multi-task multi-kernel learning framework capable of harmonizing functional results and structural information with clinical phenotypes. This novel approach extends the human phenotype ontology towards kernel-based supervised machine learning. Our gain- or loss-of-function classifier achieves high performance (mean accuracy 0.853 SD 0.016, mean AU-ROC 0.912 SD 0.025), outperforming both conventional baseline and state-of-the-art methods. Performance is robust across different phenotypic similarity measures and largely insensitive to phenotypic noise or sparsity. Localized multi-kernel learning offered biological insight and interpretability by highlighting channels with implicit genotype-phenotype correlations or latent task similarity for downstream analysis.

Repetitive transcranial magnetic stimulation (rTMS) is a promising alternative therapy for treatment-resistant depression, although its limited remission rate indicates room for improvement. As depression is a phenomenological construction, the biological heterogeneity within this syndrome needs to be considered to improve the existing therapies. Whole-brain modeling provides an integrative multi-modal framework for capturing disease heterogeneity in a holistic manner. Computational modelling combined with probabilistic nonparametric fitting was applied to the resting-state fMRI data from 42 patients (21 women), to parametrize baseline brain dynamics in depression. All patients were randomly assigned to two treatment groups, namely active (i.e., rTMS, n = 22) or sham (n = 20). The active treatment group received rTMS treatment with an accelerated intermittent theta burst protocol over the dorsomedial prefrontal cortex. The sham treatment group underwent the identical procedure but with the magnetically shielded side of the coil. We stratified the depression sample into distinct covert subtypes based on their baseline attractor dynamics captured by different model parameters. Notably, the two detected depression subtypes exhibited different phenotypic behaviors at baseline. Our stratification could predict the diverse response to the active treatment that could not be explained by the sham treatment. Critically, we further found that one group exhibited more distinct improvement in certain affective and negative symptoms. The subgroup of patients with higher responsiveness to treatment exhibited blunted frequency dynamics for intrinsic activity at baseline, as indexed by lower global metastability and synchrony. Our findings suggested that whole-brain modeling of intrinsic dynamics may constitute a determinant for stratifying patients into treatment groups and bringing us closer towards precision medicine.

We present a multiagent-based model that captures the interactions between different types of cells with their microenvironment, and enables the analysis of the emergent global behavior during tissue regeneration and tumor development. Using this model, we are able to reproduce the temporal dynamics of regular healthy cells and cancer cells, as well as the evolution of their three-dimensional spatial distributions. By tuning the system with the characteristics of the individual patients, our model reproduces a variety of spatial patterns of tissue regeneration and tumor growth, resembling those found in clinical imaging or biopsies. In order to calibrate and validate our model we study the process of liver regeneration after surgical hepatectomy in different degrees. In the clinical context, our model is able to predict the recurrence of a hepatocellular carcinoma after a 70% partial hepatectomy. The outcomes of our simulations are in agreement with experimental and clinical observations. By fitting the model parameters to specific patient factors, it might well become a useful platform for hypotheses testing in treatments protocols.

Diffuse low grade gliomas are invasive and incurable brain tumors that inevitably transform into higher grade ones. A classical treatment to delay this transition is radiotherapy (RT). Following RT, the tumor gradually shrinks during a period of typically 6 months to 4 years before regrowing. To improve the patient's health-related quality of life and help clinicians build personalized follow-ups, one would benefit from predictions of the time during which the tumor is expected to decrease. The challenge is to provide a reliable estimate of this regrowth time shortly after RT (i.e. with few data), although patients react differently to the treatment. To this end, we analyze the tumor size dynamics from a batch of 20 high-quality longitudinal data, and propose a simple and robust analytical model, with just 4 parameters. From the study of their correlations, we build a statistical constraint that helps determine the regrowth time even for patients for which we have only a few measurements of the tumor size. We validate the procedure on the data and predict the regrowth time at the moment of the first MRI after RT, with precision of, typically, 6 months. Using virtual patients, we study whether some forecast is still possible just three months after RT. We obtain some reliable estimates of the regrowth time in 75% of the cases, in particular for all "fast-responders". The remaining 25% represent cases where the actual regrowth time is large and can be safely estimated with another measurement a year later. These results show the feasibility of making personalized predictions of the tumor regrowth time shortly after RT.

Here we introduce chiLife, a Python package for site-directed spin label (SDSL) modeling for electron paramagnetic resonance (EPR) spectroscopy, in particular double electron-electron resonance (DEER). It is based on in silico attachment of rotamer ensemble representations of spin labels to protein structures. chiLife enables the development of custom protein analysis and modeling pipelines using SDSL EPR experimental data. It allows the user to add custom spin labels, scoring functions and spin label modeling methods. chiLife is designed with integration into third-party software in mind, to take advantage of the diverse and rapidly expanding set of molecular modeling tools available with a Python interface. This article describes the main design principles of chiLife and presents a series of examples.

Advances in long-read sequencing technologies have dramatically improved the contiguity and completeness of genome assemblies. Using the latest nanopore-based sequencers, we can generate enough data for the assembly of a human genome from a single flow cell. With the long-read data from these sequences, we can now routinely produce de novo genome assemblies in which half or more of a genome is contained in megabase-scale contigs. Assemblies produced from nanopore data alone, though, have relatively high error rates and can benefit from a process called polishing, in which more-accurate reads are used to correct errors in the consensus sequence. In this manuscript, we present a novel tool for genome polishing called JASPER (Jellyfish-based Assembly Sequence Polisher for Error Reduction). In contrast to many other polishing methods, JASPER gains efficiency by avoiding the alignment of reads to the assembly. Instead, JASPER uses a database of k-mer counts that it creates from the reads to detect and correct errors in the consensus. Our experiments demonstrate that JASPER is faster than alignment-based polishers, and both faster and more accurate than other k-mer based polishing methods. We also introduce the idea of using a polishing tool to create population-specific reference genomes, and illustrate this idea using sequence data from multiple individuals from Tokyo, Japan.

Drug-target binding affinity prediction plays a key role in the early stage of drug discovery. Numerous experimental and data-driven approaches have been developed for predicting drug-target binding affinity. However, experimental methods highly rely on the limited structural-related information from drug-target pairs, domain knowledge, and time-consuming assays. On the other hand, learning-based methods have shown an acceptable prediction performance. However, most of them utilize several simple and complex types of proteins and drug compounds data, ranging from the protein sequences to the topology of a graph representation of drug compounds, employing multiple deep neural networks for encoding and feature extraction, and so, leads to the computational overheads. In this study, we propose a unified measure for protein sequence encoding, named BiComp, which provides compression-based and evolutionary-related features from the protein sequences. Specifically, we employ Normalized Compression Distance and Smith-Waterman measures for capturing complementary information from the algorithmic information theory and biological domains, respectively. We utilize the proposed measure to encode the input proteins feeding a new deep neural network-based method for drug-target binding affinity prediction, named BiComp-DTA. BiComp-DTA is evaluated utilizing four benchmark datasets for drug-target binding affinity prediction. Compared to the state-of-the-art methods, which employ complex models for protein encoding and feature extraction, BiComp-DTA provides superior efficiency in terms of accuracy, runtime, and the number of trainable parameters. The latter achievement facilitates execution of BiComp-DTA on a normal desktop computer in a fast fashion. As a comparative study, we evaluate BiComp's efficiency against its components for drug-target binding affinity prediction. The results have shown superior accuracy of BiComp due to the orthogonality and complementary nature of Smith-Waterman and Normalized Compression Distance measures for protein sequences. Such a protein sequence encoding provides efficient representation with no need for multiple sources of information, deep domain knowledge, and complex neural networks.

Although some methods for estimating the instantaneous reproductive number during epidemics have been developed, the existing frameworks usually require information on the distribution of the serial interval and/or additional contact tracing data. However, in the case of outbreaks of emerging infectious diseases with an unknown natural history or undetermined characteristics, the serial interval and/or contact tracing data are often not available, resulting in inaccurate estimates for this quantity. In the present study, a new framework was specifically designed for joint estimates of the instantaneous reproductive number and serial interval. Concretely, a likelihood function for the two quantities was first introduced. Then, the instantaneous reproductive number and the serial interval were modeled parametrically as a function of time using the interpolation method and a known traditional distribution, respectively. Using the Bayesian information criterion and the Markov Chain Monte Carlo method, we ultimately obtained their estimates and distribution. The simulation study revealed that our estimates of the two quantities were consistent with the ground truth. Seven data sets of historical epidemics were considered and further verified the robust performance of our method. Therefore, to some extent, even if we know only the daily incidence, our method can accurately estimate the instantaneous reproductive number and serial interval to provide crucial information for policymakers to design appropriate prevention and control interventions during epidemics.

Predicting protein-protein interactions from sequences is an important goal of computational biology. Various sources of information can be used to this end. Starting from the sequences of two interacting protein families, one can use phylogeny or residue coevolution to infer which paralogs are specific interaction partners within each species. We show that these two signals can be combined to improve the performance of the inference of interaction partners among paralogs. For this, we first align the sequence-similarity graphs of the two families through simulated annealing, yielding a robust partial pairing. We next use this partial pairing to seed a coevolution-based iterative pairing algorithm. This combined method improves performance over either separate method. The improvement obtained is striking in the difficult cases where the average number of paralogs per species is large or where the total number of sequences is modest.

Progress continues in the field of cancer biology, yet much remains to be unveiled regarding the mechanisms of cancer invasion. In particular, complex biophysical mechanisms enable a tumor to remodel the surrounding extracellular matrix (ECM), allowing cells to invade alone or collectively. Tumor spheroids cultured in collagen represent a simplified, reproducible 3D model system, which is sufficiently complex to recapitulate the evolving organization of cells and interaction with the ECM that occur during invasion. Recent experimental approaches enable high resolution imaging and quantification of the internal structure of invading tumor spheroids. Concurrently, computational modeling enables simulations of complex multicellular aggregates based on first principles. The comparison between real and simulated spheroids represents a way to fully exploit both data sources, but remains a challenge. We hypothesize that comparing any two spheroids requires first the extraction of basic features from the raw data, and second the definition of key metrics to match such features. Here, we present a novel method to compare spatial features of spheroids in 3D. To do so, we define and extract features from spheroid point cloud data, which we simulated using Cells in Silico (CiS), a high-performance framework for large-scale tissue modeling previously developed by us. We then define metrics to compare features between individual spheroids, and combine all metrics into an overall deviation score. Finally, we use our features to compare experimental data on invading spheroids in increasing collagen densities. We propose that our approach represents the basis for defining improved metrics to compare large 3D data sets. Moving forward, this approach will enable the detailed analysis of spheroids of any origin, one application of which is informing in silico spheroids based on their in vitro counterparts. This will enable both basic and applied researchers to close the loop between modeling and experiments in cancer research.

We analyzed large-scale post-translational modification (PTM) data to outline cell signaling pathways affected by tyrosine kinase inhibitors (TKIs) in ten lung cancer cell lines. Tyrosine phosphorylated, lysine ubiquitinated, and lysine acetylated proteins were concomitantly identified using sequential enrichment of post translational modification (SEPTM) proteomics. Machine learning was used to identify PTM clusters that represent functional modules that respond to TKIs. To model lung cancer signaling at the protein level, PTM clusters were used to create a co-cluster correlation network (CCCN) and select protein-protein interactions (PPIs) from a large network of curated PPIs to create a cluster-filtered network (CFN). Next, we constructed a Pathway Crosstalk Network (PCN) by connecting pathways from NCATS BioPlanet whose member proteins have PTMs that co-cluster. Interrogating the CCCN, CFN, and PCN individually and in combination yields insights into the response of lung cancer cells to TKIs. We highlight examples where cell signaling pathways involving EGFR and ALK exhibit crosstalk with BioPlanet pathways: Transmembrane transport of small molecules; and Glycolysis and gluconeogenesis. These data identify known and previously unappreciated connections between receptor tyrosine kinase (RTK) signal transduction and oncogenic metabolic reprogramming in lung cancer. Comparison to a CFN generated from a previous multi-PTM analysis of lung cancer cell lines reveals a common core of PPIs involving heat shock/chaperone proteins, metabolic enzymes, cytoskeletal components, and RNA-binding proteins. Elucidation of points of crosstalk among signaling pathways employing different PTMs reveals new potential drug targets and candidates for synergistic attack through combination drug therapy.

As researchers develop computational models of neural systems with increasing sophistication and scale, it is often the case that fully de novo model development is impractical and inefficient. Thus arises a critical need to quickly find, evaluate, re-use, and build upon models and model components developed by other researchers. We introduce the NeuroML Database (NeuroML-DB.org), which has been developed to address this need and to complement other model sharing resources. NeuroML-DB stores over 1,500 previously published models of ion channels, cells, and networks that have been translated to the modular NeuroML model description language. The database also provides reciprocal links to other neuroscience model databases (ModelDB, Open Source Brain) as well as access to the original model publications (PubMed). These links along with Neuroscience Information Framework (NIF) search functionality provide deep integration with other neuroscience community modeling resources and greatly facilitate the task of finding suitable models for reuse. Serving as an intermediate language, NeuroML and its tooling ecosystem enable efficient translation of models to other popular simulator formats. The modular nature also enables efficient analysis of a large number of models and inspection of their properties. Search capabilities of the database, together with web-based, programmable online interfaces, allow the community of researchers to rapidly assess stored model electrophysiology, morphology, and computational complexity properties. We use these capabilities to perform a database-scale analysis of neuron and ion channel models and describe a novel tetrahedral structure formed by cell model clusters in the space of model properties and features. This analysis provides further information about model similarity to enrich database search.

Understanding how neural populations encode sensory stimuli remains a central problem in neuroscience. Here we performed multi-unit recordings from sensory neural populations in the electrosensory system of the weakly electric fish Apteronotus leptorhynchus in response to stimuli located at different positions along the rostro-caudal axis. Our results reveal that the spatial dependence of correlated activity along receptive fields can help mitigate the deleterious effects that these correlations would otherwise have if they were spatially independent. Moreover, using mathematical modeling, we show that experimentally observed heterogeneities in the receptive fields of neurons help optimize information transmission as to object location. Taken together, our results have important implications for understanding how sensory neurons whose receptive fields display antagonistic center-surround organization encode location. Important similarities between the electrosensory system and other sensory systems suggest that our results will be applicable elsewhere.

Accurate prediction of synergistic effects of drug combinations can reduce the experimental costs for drug development and facilitate the discovery of novel efficacious combination therapies for clinical studies. The drug combinations with high synergy scores are regarded as synergistic ones, while those with moderate or low synergy scores are additive or antagonistic ones. The existing methods usually exploit the synergy data from the aspect of synergistic drug combinations, paying little attention to the additive or antagonistic ones. Also, they usually do not leverage the common patterns of drug combinations across different cell lines. In this paper, we propose a multi-channel graph autoencoder (MGAE)-based method for predicting the synergistic effects of drug combinations (DC), and shortly denote it as MGAE-DC. A MGAE model is built to learn the drug embeddings by considering not only synergistic combinations but also additive and antagonistic ones as three input channels. The later two channels guide the model to explicitly characterize the features of non-synergistic combinations through an encoder-decoder learning process, and thus the drug embeddings become more discriminative between synergistic and non-synergistic combinations. In addition, an attention mechanism is incorporated to fuse each cell-line's drug embeddings across various cell lines, and a common drug embedding is extracted to capture the invariant patterns by developing a set of cell-line shared decoders. The generalization performance of our model is further improved with the invariant patterns. With the cell-line specific and common drug embeddings, our method is extended to predict the synergy scores of drug combinations by a neural network module. Experiments on four benchmark datasets demonstrate that MGAE-DC consistently outperforms the state-of-the-art methods. In-depth literature survey is conducted to find that many drug combinations predicted by MGAE-DC are supported by previous experimental studies. The source code and data are available at https://github.com/yushenshashen/MGAE-DC.

Mathematical models have been an important tool during the COVID-19 pandemic, for example to predict demand of critical resources such as medical devices, personal protective equipment and diagnostic tests. Many COVID-19 models have been developed. However, there is relatively little information available regarding reliability of model predictions. Here we present a general model validation framework for epidemiological models focused around predictive capability for questions relevant to decision-making end-users. COVID-19 models are typically comprised of multiple releases, and provide predictions for multiple localities, and these characteristics are systematically accounted for in the framework, which is based around a set of validation scores or metrics that quantify model accuracy of specific quantities of interest including: date of peak, magnitude of peak, rate of recovery, and monthly cumulative counts. We applied the framework to retrospectively assess accuracy of death predictions for four COVID-19 models, and accuracy of hospitalization predictions for one COVID-19 model (models for which sufficient data was publicly available). When predicting date of peak deaths, the most accurate model had errors of approximately 15 days or less, for releases 3-6 weeks in advance of the peak. Death peak magnitude relative errors were generally in the 50% range 3-6 weeks before peak. Hospitalization predictions were less accurate than death predictions. All models were highly variable in predictive accuracy across regions. Overall, our framework provides a wealth of information on the predictive accuracy of epidemiological models and could be used in future epidemics to evaluate new models or support existing modeling methodologies, and thereby aid in informed model-based public health decision making. The code for the validation framework is available at https://doi.org/10.5281/zenodo.7102854.

[This corrects the article DOI: 10.1371/journal.pcbi.1009398.].

Plant pathogens respond to selection pressures exerted by disease management strategies. This can lead to fungicide resistance and/or the breakdown of disease-resistant cultivars, each of which significantly threaten food security. Both fungicide resistance and cultivar breakdown can be characterised as qualitative or quantitative. Qualitative (monogenic) resistance/breakdown involves a step change in the characteristics of the pathogen population with respect to disease control, often caused by a single genetic change. Quantitative (polygenic) resistance/breakdown instead involves multiple genetic changes, each causing a smaller shift in pathogen characteristics, leading to a gradual alteration in the effectiveness of disease control over time. Although resistance/breakdown to many fungicides/cultivars currently in use is quantitative, the overwhelming majority of modelling studies focus on the much simpler case of qualitative resistance. Further, those very few models of quantitative resistance/breakdown which do exist are not fitted to field data. Here we present a model of quantitative resistance/breakdown applied to Zymoseptoria tritici, which causes Septoria leaf blotch, the most prevalent disease of wheat worldwide. Our model is fitted to data from field trials in the UK and Denmark. For fungicide resistance, we show that the optimal disease management strategy depends on the timescale of interest. Greater numbers of fungicide applications per year lead to greater selection for resistant strains, although over short timescales this can be oset by the increased control oered by more sprays. However, over longer timescales higher yields are attained using fewer fungicide applications per year. Deployment of disease-resistant cultivars is not only a valuable disease management strategy, but also oers the secondary benefit of protecting fungicide effectiveness by delaying the development of fungicide resistance. However, disease-resistant cultivars themselves erode over time. We show how an integrated disease management strategy with frequent replacement of disease-resistant cultivars can give a large improvement in fungicide durability and yields.

Lipid molecules such as cholesterol interact with the surface of integral membrane proteins (IMP) in a mode different from drug-like molecules in a protein binding pocket. These differences are due to the lipid molecule's shape, the membrane's hydrophobic environment, and the lipid's orientation in the membrane. We can use the recent increase in experimental structures in complex with cholesterol to understand protein-cholesterol interactions. We developed the RosettaCholesterol protocol consisting of (1) a prediction phase using an energy grid to sample and score native-like binding poses and (2) a specificity filter to calculate the likelihood that a cholesterol interaction site may be specific. We used a multi-pronged benchmark (self-dock, flip-dock, cross-dock, and global-dock) of protein-cholesterol complexes to validate our method. RosettaCholesterol improved sampling and scoring of native poses over the standard RosettaLigand baseline method in 91% of cases and performs better regardless of benchmark complexity. On the β2AR, our method found one likely-specific site, which is described in the literature. The RosettaCholesterol protocol quantifies cholesterol binding site specificity. Our approach provides a starting point for high-throughput modeling and prediction of cholesterol binding sites for further experimental validation.

A model for DNA and nucleosomes is introduced with the goal of studying chromosomes from a single base level all the way to higher-order chromatin structures. This model, dubbed the Widely Editable Chromatin Model (WEChroM), reproduces the complex mechanics of the double helix including its bending persistence length and twisting persistence length, and the temperature dependence of the former. The WEChroM Hamiltonian is composed of chain connectivity, steric interactions, and associative memory terms representing all remaining interactions leading to the structure, dynamics, and mechanical characteristics of the B-DNA. Several applications of this model are discussed to demonstrate its applicability. WEChroM is used to investigate the behavior of circular DNA in the presence of positive and negative supercoiling. We show that it recapitulates the formation of plectonemes and of structural defects that relax mechanical stress. The model spontaneously manifests an asymmetric behavior with respect to positive or negative supercoiling, similar to what was previously observed in experiments. Additionally, we show that the associative memory Hamiltonian is also capable of reproducing the free energy of partial DNA unwrapping from nucleosomes. WEChroM is designed to emulate the continuously variable mechanical properties of the 10nm fiber and, by virtue of its simplicity, is ready to be scaled up to molecular systems large enough to investigate the structural ensembles of genes. WEChroM is implemented in the OpenMM simulation toolkits and is freely available for public use.

Displaced communication, whereby individuals communicate regarding a subject that is not immediately present (spatially or temporally), is one of the key features of human language. It also occurs in a few animal species, most notably the honeybee, where the waggle dance is used to communicate the location and quality of a patch of flowers. However, it is difficult to study how it emerged given the paucity of species displaying this capacity and the fact that it often occurs via complex multimodal signals. To address this issue, we developed a novel paradigm in which we conducted experimental evolution with foraging agents endowed with neural networks that regulate their movement and the production of signals. Displaced communication readily evolved but, surprisingly, agents did not use signal amplitude to convey information on food location. Instead, they used signal onset-delay and duration-based mode of communication, which depends on the motion of the agent within a communication area. When agents were experimentally prevented from using these modes of communication, they evolved to use signal amplitude instead. Interestingly, this mode of communication was more efficient and led to higher performance. Subsequent controlled experiments suggested that this more efficient mode of communication failed to evolve because it took more generations to emerge than communication grounded on the onset-delay and length of signaling. These results reveal that displaced communication is likely to initially evolve from non-communicative behavioral cues providing incidental information with evolution later leading to more efficient communication systems through a ritualization process.

Those building predictive models from transcriptomic data are faced with two conflicting perspectives. The first, based on the inherent high dimensionality of biological systems, supposes that complex non-linear models such as neural networks will better match complex biological systems. The second, imagining that complex systems will still be well predicted by simple dividing lines prefers linear models that are easier to interpret. We compare multi-layer neural networks and logistic regression across multiple prediction tasks on GTEx and Recount3 datasets and find evidence in favor of both possibilities. We verified the presence of non-linear signal when predicting tissue and metadata sex labels from expression data by removing the predictive linear signal with Limma, and showed the removal ablated the performance of linear methods but not non-linear ones. However, we also found that the presence of non-linear signal was not necessarily sufficient for neural networks to outperform logistic regression. Our results demonstrate that while multi-layer neural networks may be useful for making predictions from gene expression data, including a linear baseline model is critical because while biological systems are high-dimensional, effective dividing lines for predictive models may not be.

The rise of single-cell data highlights the need for a nondeterministic view of gene expression, while offering new opportunities regarding gene regulatory network inference. We recently introduced two strategies that specifically exploit time-course data, where single-cell profiling is performed after a stimulus: HARISSA, a mechanistic network model with a highly efficient simulation procedure, and CARDAMOM, a scalable inference method seen as model calibration. Here, we combine the two approaches and show that the same model driven by transcriptional bursting can be used simultaneously as an inference tool, to reconstruct biologically relevant networks, and as a simulation tool, to generate realistic transcriptional profiles emerging from gene interactions. We verify that CARDAMOM quantitatively reconstructs causal links when the data is simulated from HARISSA, and demonstrate its performance on experimental data collected on in vitro differentiating mouse embryonic stem cells. Overall, this integrated strategy largely overcomes the limitations of disconnected inference and simulation.

Machine learning models have difficulty generalizing to data outside of the distribution they were trained on. In particular, vision models are usually vulnerable to adversarial attacks or common corruptions, to which the human visual system is robust. Recent studies have found that regularizing machine learning models to favor brain-like representations can improve model robustness, but it is unclear why. We hypothesize that the increased model robustness is partly due to the low spatial frequency preference inherited from the neural representation. We tested this simple hypothesis with several frequency-oriented analyses, including the design and use of hybrid images to probe model frequency sensitivity directly. We also examined many other publicly available robust models that were trained on adversarial images or with data augmentation, and found that all these robust models showed a greater preference to low spatial frequency information. We show that preprocessing by blurring can serve as a defense mechanism against both adversarial attacks and common corruptions, further confirming our hypothesis and demonstrating the utility of low spatial frequency information in robust object recognition.

We introduce a new spatial statistic, the weighted pair correlation function (wPCF). The wPCF extends the existing pair correlation function (PCF) and cross-PCF to describe spatial relationships between points marked with combinations of discrete and continuous labels. We validate its use through application to a new agent-based model (ABM) which simulates interactions between macrophages and tumour cells. These interactions are influenced by the spatial positions of the cells and by macrophage phenotype, a continuous variable that ranges from anti-tumour to pro-tumour. By varying model parameters that regulate macrophage phenotype, we show that the ABM exhibits behaviours which resemble the 'three Es of cancer immunoediting': Equilibrium, Escape, and Elimination. We use the wPCF to analyse synthetic images generated by the ABM. We show that the wPCF generates a 'human readable' statistical summary of where macrophages with different phenotypes are located relative to both blood vessels and tumour cells. We also define a distinct 'PCF signature' that characterises each of the three Es of immunoediting, by combining wPCF measurements with the cross-PCF describing interactions between vessels and tumour cells. By applying dimension reduction techniques to this signature, we identify its key features and train a support vector machine classifier to distinguish between simulation outputs based on their PCF signature. This proof-of-concept study shows how multiple spatial statistics can be combined to analyse the complex spatial features that the ABM generates, and to partition them into interpretable groups. The intricate spatial features produced by the ABM are similar to those generated by state-of-the-art multiplex imaging techniques which distinguish the spatial distribution and intensity of multiple biomarkers in biological tissue regions. Applying methods such as the wPCF to multiplex imaging data would exploit the continuous variation in biomarker intensities and generate more detailed characterisation of the spatial and phenotypic heterogeneity in tissue samples.

Surface antigens of pathogens are commonly targeted by vaccine-elicited antibodies but antigenic variability, notably in RNA viruses such as influenza, HIV and SARS-CoV-2, pose challenges for control by vaccination. For example, influenza A(H3N2) entered the human population in 1968 causing a pandemic and has since been monitored, along with other seasonal influenza viruses, for the emergence of antigenic drift variants through intensive global surveillance and laboratory characterisation. Statistical models of the relationship between genetic differences among viruses and their antigenic similarity provide useful information to inform vaccine development, though accurate identification of causative mutations is complicated by highly correlated genetic signals that arise due to the evolutionary process. Here, using a sparse hierarchical Bayesian analogue of an experimentally validated model for integrating genetic and antigenic data, we identify the genetic changes in influenza A(H3N2) virus that underpin antigenic drift. We show that incorporating protein structural data into variable selection helps resolve ambiguities arising due to correlated signals, with the proportion of variables representing haemagglutinin positions decisively included, or excluded, increased from 59.8% to 72.4%. The accuracy of variable selection judged by proximity to experimentally determined antigenic sites was improved simultaneously. Structure-guided variable selection thus improves confidence in the identification of genetic explanations of antigenic variation and we also show that prioritising the identification of causative mutations is not detrimental to the predictive capability of the analysis. Indeed, incorporating structural information into variable selection resulted in a model that could more accurately predict antigenic assay titres for phenotypically-uncharacterised virus from genetic sequence. Combined, these analyses have the potential to inform choices of reference viruses, the targeting of laboratory assays, and predictions of the evolutionary success of different genotypes, and can therefore be used to inform vaccine selection processes.

Neural mass models (NMMs) are important for helping us interpret observations of brain dynamics. They provide a means to understand data in terms of mechanisms such as synaptic interactions between excitatory and inhibitory neuronal populations. To interpret data using NMMs we need to quantitatively compare the output of NMMs with data, and thereby find parameter values for which the model can produce the observed dynamics. Mapping dynamics to NMM parameter values in this way has the potential to improve our understanding of the brain in health and disease. Though abstract, NMMs still comprise of many parameters that are difficult to constrain a priori. This makes it challenging to explore the dynamics of NMMs and elucidate regions of parameter space in which their dynamics best approximate data. Existing approaches to overcome this challenge use a combination of linearising models, constraining the values they can take and exploring restricted subspaces by fixing the values of many parameters a priori. As such, we have little knowledge of the extent to which different regions of parameter space of NMMs can yield dynamics that approximate data, how nonlinearities in models can affect parameter mapping or how best to quantify similarities between model output and data. These issues need to be addressed in order to fully understand the potential and limitations of NMMs, and to aid the development of new models of brain dynamics in the future. To begin to overcome these issues, we present a global nonlinear approach to recovering parameters of NMMs from data. We use global optimisation to explore all parameters of nonlinear NMMs simultaneously, in a minimally constrained way. We do this using multi-objective optimisation (multi-objective evolutionary algorithm, MOEA) so that multiple data features can be quantified. In particular, we use the weighted horizontal visibility graph (wHVG), which is a flexible framework for quantifying different aspects of time series, by converting them into networks. We study EEG alpha activity recorded during the eyes closed resting state from 20 healthy individuals and demonstrate that the MOEA performs favourably compared to single objective approaches. The addition of the wHVG objective allows us to better constrain the model output, which leads to the recovered parameter values being restricted to smaller regions of parameter space, thus improving the practical identifiability of the model. We then use the MOEA to study differences in the alpha rhythm observed in EEG recorded from 20 people with epilepsy. We find that a small number of parameters can explain this difference and that, counterintuitively, the mean excitatory synaptic gain parameter is reduced in people with epilepsy compared to control. In addition, we propose that the MOEA could be used to mine for the presence of pathological rhythms, and demonstrate the application of this to epileptiform spike-wave discharges.

Social insects demonstrate adaptive behaviour for a given colony size. Remarkably, most species do this even without visual information in a dark environment. However, how they achieve this is yet unknown. Based on individual trait expression, an agent-based simulation was used to identify an explicit mechanism for understanding colony size dependent behaviour. Through repeated physical contact between the queen and individual workers, individual colony members monitor the physiological states of others, reflecting such contact information in their physiology and behaviour. Feedback between the sensing of physiological states and the corresponding behaviour patterns leads to self-organisation with colonies shifting according to their size. We showed (1) the queen can exhibit adaptive behaviour patterns for the increase in colony size while density per space remains unchanged, and (2) such physical constraints can underlie the adaptive switching of colony stages from successful patrol behaviour to unsuccessful patrol behaviour, which leads to constant ovary development (production of reproductive castes). The feedback loops embedded in the queen between the perception of internal states of the workers and behavioural patterns can explain the adaptive behaviour as a function of colony size.

Understanding how pollinators move across space is key to understanding plant mating patterns. Bees are typically assumed to search for flowers randomly or using simple movement rules, so that the probability of discovering a flower should primarily depend on its distance to the nest. However, experimental work shows this is not always the case. Here, we explored the influence of flower size and density on their probability of being discovered by bees by developing a movement model of central place foraging bees, based on experimental data collected on bumblebees. Our model produces realistic bee trajectories by taking into account the autocorrelation of the bee's angular speed, the attraction to the nest (homing), and a gaussian noise. Simulations revealed a « masking effect » that reduces the detection of flowers close to another, with potential far reaching consequences on plant-pollinator interactions. At the plant level, flowers distant to the nest were more often discovered by bees in low density environments. At the bee colony level, foragers found more flowers when they were small and at medium densities. Our results indicate that the processes of search and discovery of resources are potentially more complex than usually assumed, and question the importance of resource distribution and abundance on bee foraging success and plant pollination.

Thermodynamic stability is a crucial fitness constraint in protein evolution and is a central factor in shaping the sequence landscapes of proteins. The correlation between stability and molecular fitness depends on the mechanism that relates the biophysical property with biological function. In the simplest case, stability and fitness are related by the amount of folded protein. However, when proteins are toxic in the unfolded state, the fitness function shifts, resulting in higher stability under mutation-selection balance. Likewise, a higher population size results in a similar change in protein stability, as it magnifies the effect of the selection pressure in evolutionary dynamics. This study investigates how such factors affect the evolution of protein stability, site-specific mutation rates, and residue-residue covariation. To simulate evolutionary trajectories with realistic modeling of protein energetics, we develop an all-atom simulator of protein evolution, RosettaEvolve. By evolving proteins under different fitness functions, we can study how the fitness function affects the distribution of proposed and accepted mutations, site-specific rates, and the prevalence of correlated amino acid substitutions. We demonstrate that fitness pressure affects the proposal distribution of mutational effects, that changes in stability can largely explain variations in site-specific substitution rates in evolutionary trajectories, and that increased fitness pressure results in a stronger covariation signal. Our results give mechanistic insight into the evolutionary consequences of variation in protein stability and provide a basis to rationalize the strong covariation signal observed in natural sequence alignments.

Switch-like responses arising from bistability have been linked to cell signaling processes and memory. Revealing the shape and properties of the set of parameters that lead to bistability is necessary to understand the underlying biological mechanisms, but is a complex mathematical problem. We present an efficient approach to address a basic topological property of the parameter region of multistationary, namely whether it is connected. The connectivity of this region can be interpreted in terms of the biological mechanisms underlying bistability and the switch-like patterns that the system can create. We provide an algorithm to assert that the parameter region of multistationarity is connected, targeting reaction networks with mass-action kinetics. We show that this is the case for numerous relevant cell signaling motifs, previously described to exhibit bistability. The method relies on linear programming and bypasses the expensive computational cost of direct and generic approaches to study parametric polynomial systems. This characteristic makes it suitable for mass-screening of reaction networks. Although the algorithm can only be used to certify connectivity, we illustrate that the ideas behind the algorithm can be adapted on a case-by-case basis to also decide that the region is not connected. In particular, we show that for a motif displaying a phosphorylation cycle with allosteric enzyme regulation, the region of multistationarity has two distinct connected components, corresponding to two different, but symmetric, biological mechanisms.

Guidelines concerning the potentially harmful effects of scientific studies have historically focused on ethical considerations for minimizing risk for participants. However, studies can also indirectly inflict harm on individuals and social groups through how they are designed, reported, and disseminated. As evidenced by recent criticisms and retractions of high-profile studies dealing with a wide variety of social issues, there is a scarcity of resources and guidance on how one can conduct research in a socially responsible manner. As such, even motivated researchers might publish work that has negative social impacts due to a lack of awareness. To address this, we propose 10 simple rules for researchers who wish to conduct socially responsible science. These rules, which cover major considerations throughout the life cycle of a study from inception to dissemination, are not aimed as a prescriptive list or a deterministic code of conduct. Rather, they are meant to help motivated scientists to reflect on their social responsibility as researchers and actively engage with the potential social impact of their research.

Phase amplitude coupling (PAC) between slow and fast oscillations is found throughout the brain and plays important functional roles. Its neural origin remains unclear. Experimental findings are often puzzling and sometimes contradictory. Most computational models rely on pairs of pacemaker neurons or neural populations tuned at different frequencies to produce PAC. Here, using a data-driven model of a hippocampal microcircuit, we demonstrate that PAC can naturally emerge from a single feedback mechanism involving an inhibitory and excitatory neuron population, which interplay to generate theta frequency periodic bursts of higher frequency gamma. The model suggests the conditions under which a CA1 microcircuit can operate to elicit theta-gamma PAC, and highlights the modulatory role of OLM and PVBC cells, recurrent connectivity, and short term synaptic plasticity. Surprisingly, the results suggest the experimentally testable prediction that the generation of the slow population oscillation requires the fast one and cannot occur without it.

Despite the vast phenotypic differences observed across primates, their protein products are largely similar to each other at the sequence level. We hypothesized that, since proteins accomplish all their functions via interactions with other molecules, alterations in the sites that participate in these interactions may be of critical importance. To uncover the extent to which these sites evolve across primates, we built a structurally-derived dataset of ~4,200 one-to-one orthologous sequence groups across 18 primate species, consisting of ~68,000 ligand-binding sites that interact with DNA, RNA, small molecules, ions, or peptides. Using this dataset, we identify functionally important patterns of conservation and variation within the amino acid residues that facilitate protein-ligand interactions across the primate phylogeny. We uncover that interaction sites are significantly more conserved than other sites, and that sites binding DNA and RNA further exhibit the lowest levels of variation. We also show that the subset of ligand-binding sites that do vary are enriched in components of gene regulatory pathways and uncover several instances of human-specific ligand-binding site changes within transcription factors. Altogether, our results suggest that ligand-binding sites have experienced selective pressure in primates and propose that variation in these sites may have an outsized effect on phenotypic variation in primates through pleiotropic effects on gene regulation.

Mitochondria are highly dynamic organelles, containing vital populations of mitochondrial DNA (mtDNA) distributed throughout the cell. Mitochondria form diverse physical structures in different cells, from cell-wide reticulated networks to fragmented individual organelles. These physical structures are known to influence the genetic makeup of mtDNA populations between cell divisions, but their influence on the inheritance of mtDNA at divisions remains less understood. Here, we use statistical and computational models of mtDNA content inside and outside the reticulated network to quantify how mitochondrial network structure can control the variances of inherited mtDNA copy number and mutant load. We assess the use of moment-based approximations to describe heteroplasmy variance and identify several cases where such an approach has shortcomings. We show that biased inclusion of one mtDNA type in the network can substantially increase heteroplasmy variance (acting as a genetic bottleneck), and controlled distribution of network mass and mtDNA through the cell can conversely reduce heteroplasmy variance below a binomial inheritance picture. Network structure also allows the generation of heteroplasmy variance while controlling copy number inheritance to sub-binomial levels, reconciling several observations from the experimental literature. Overall, different network structures and mtDNA arrangements within them can control the variances of key variables to suit a palette of different inheritance priorities.

There exists an ongoing need to improve the validity and accuracy of computational fluid dynamics (CFD) simulations of turbulent airflows in the extra-thoracic and upper airways. Yet, a knowledge gap remains in providing experimentally-resolved 3D flow benchmarks with sufficient data density and completeness for useful comparison with widely-employed numerical schemes. Motivated by such shortcomings, the present work details to the best of our knowledge the first attempt to deliver in vitro-in silico correlations of 3D respiratory airflows in a generalized mouth-throat model and thereby assess the performance of Large Eddy Simulations (LES) and Reynolds-Averaged Numerical Simulations (RANS). Numerical predictions are compared against 3D volumetric flow measurements using Tomographic Particle Image Velocimetry (TPIV) at three steady inhalation flowrates varying from shallow to deep inhalation conditions. We find that a RANS k-ω SST model adequately predicts velocity flow patterns for Reynolds numbers spanning 1'500 to 7'000, supporting results in close proximity to a more computationally-expensive LES model. Yet, RANS significantly underestimates turbulent kinetic energy (TKE), thus underlining the advantages of LES as a higher-order turbulence modeling scheme. In an effort to bridge future endevours across respiratory research disciplines, we provide end users with the present in vitro-in silico correlation data for improved predictive CFD models towards inhalation therapy and therapeutic or toxic dosimetry endpoints.

COVID-19 is one of the deadliest respiratory diseases, and its emergence caught the pharmaceutical industry off guard. While vaccines have been rapidly developed, treatment options for infected people remain scarce, and COVID-19 poses a substantial global threat. This study presents a novel workflow to predict robust druggable targets against emerging RNA viruses using metabolic networks and information of the viral structure and its genome sequence. For this purpose, we implemented pymCADRE and PREDICATE to create tissue-specific metabolic models, construct viral biomass functions and predict host-based antiviral targets from more than one genome. We observed that pymCADRE reduces the computational time of flux variability analysis for internal optimizations. We applied these tools to create a new metabolic network of primary bronchial epithelial cells infected with SARS-CoV-2 and identified enzymatic reactions with inhibitory effects. The most promising reported targets were from the purine metabolism, while targeting the pyrimidine and carbohydrate metabolisms seemed to be promising approaches to enhance viral inhibition. Finally, we computationally tested the robustness of our targets in all known variants of concern, verifying our targets' inhibitory effects. Since laboratory tests are time-consuming and involve complex readouts to track processes, our workflow focuses on metabolic fluxes within infected cells and is applicable for rapid hypothesis-driven identification of potentially exploitable antivirals concerning various viruses and host cell types.

One of the main obstacles to the successful treatment of cancer is the phenomenon of drug resistance. A common strategy to overcome resistance is the use of combination therapies. However, the space of possibilities is huge and efficient search strategies are required. Machine Learning (ML) can be a useful tool for the discovery of novel, clinically relevant anti-cancer drug combinations. In particular, deep learning (DL) has become a popular choice for modeling drug combination effects. Here, we set out to examine the impact of different methodological choices on the performance of multimodal DL-based drug synergy prediction methods, including the use of different input data types, preprocessing steps and model architectures. Focusing on the NCI ALMANAC dataset, we found that feature selection based on prior biological knowledge has a positive impact-limiting gene expression data to cancer or drug response-specific genes improved performance. Drug features appeared to be more predictive of drug response, with a 41% increase in coefficient of determination (R2) and 26% increase in Spearman correlation relative to a baseline model that used only cell line and drug identifiers. Molecular fingerprint-based drug representations performed slightly better than learned representations-ECFP4 fingerprints increased R2 by 5.3% and Spearman correlation by 2.8% w.r.t the best learned representations. In general, fully connected feature-encoding subnetworks outperformed other architectures. DL outperformed other ML methods by more than 35% (R2) and 14% (Spearman). Additionally, an ensemble combining the top DL and ML models improved performance by about 6.5% (R2) and 4% (Spearman). Using a state-of-the-art interpretability method, we showed that DL models can learn to associate drug and cell line features with drug response in a biologically meaningful way. The strategies explored in this study will help to improve the development of computational methods for the rational design of effective drug combinations for cancer therapy.

Missing observations in trait datasets pose an obstacle for analyses in myriad biological disciplines. Considering the mixed results of imputation, the wide variety of available methods, and the varied structure of real trait datasets, a framework for selecting a suitable imputation method is advantageous. We invoked a real data-driven simulation strategy to select an imputation method for a given mixed-type (categorical, count, continuous) target dataset. Candidate methods included mean/mode imputation, k-nearest neighbour, random forests, and multivariate imputation by chained equations (MICE). Using a trait dataset of squamates (lizards and amphisbaenians; order: Squamata) as a target dataset, a complete-case dataset consisting of species with nearly complete information was formed for the imputation method selection. Missing data were induced by removing values from this dataset under different missingness mechanisms: missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). For each method, combinations with and without phylogenetic information from single gene (nuclear and mitochondrial) or multigene trees were used to impute the missing values for five numerical and two categorical traits. The performances of the methods were evaluated under each missing mechanism by determining the mean squared error and proportion falsely classified rates for numerical and categorical traits, respectively. A random forest method supplemented with a nuclear-derived phylogeny resulted in the lowest error rates for the majority of traits, and this method was used to impute missing values in the original dataset. Data with imputed values better reflected the characteristics and distributions of the original data compared to complete-case data. However, caution should be taken when imputing trait data as phylogeny did not always improve performance for every trait and in every scenario. Ultimately, these results support the use of a real data-driven simulation strategy for selecting a suitable imputation method for a given mixed-type trait dataset.

Catalysis and fidelity of multisubunit RNA polymerases rely on a highly conserved active site domain called the trigger loop (TL), which achieves roles in transcription through conformational changes and interaction with NTP substrates. The mutations of TL residues cause distinct effects on catalysis including hypo- and hyperactivity and altered fidelity. We applied molecular dynamics simulation (MD) and machine learning (ML) techniques to characterize TL mutations in the Saccharomyces cerevisiae RNA Polymerase II (Pol II) system. We did so to determine relationships between individual mutations and phenotypes and to associate phenotypes with MD simulated structural alterations. Using fitness values of mutants under various stress conditions, we modeled phenotypes along a spectrum of continual values. We found that ML could predict the phenotypes with 0.68 R2 correlation from amino acid sequences alone. It was more difficult to incorporate MD data to improve predictions from machine learning, presumably because MD data is too noisy and possibly incomplete to directly infer functional phenotypes. However, a variational auto-encoder model based on the MD data allowed the clustering of mutants with different phenotypes based on structural details. Overall, we found that a subset of loss-of-function (LOF) and lethal mutations tended to increase distances of TL residues to the NTP substrate, while another subset of LOF and lethal substitutions tended to confer an increase in distances between TL and bridge helix (BH). In contrast, some of the gain-of-function (GOF) mutants appear to cause disruption of hydrophobic contacts among TL and nearby helices.

The coalescent is a powerful statistical framework that allows us to infer past population dynamics leveraging the ancestral relationships reconstructed from sampled molecular sequence data. In many biomedical applications, such as in the study of infectious diseases, cell development, and tumorgenesis, several distinct populations share evolutionary history and therefore become dependent. The inference of such dependence is a highly important, yet a challenging problem. With advances in sequencing technologies, we are well positioned to exploit the wealth of high-resolution biological data for tackling this problem. Here, we present adaPop, a probabilistic model to estimate past population dynamics of dependent populations and to quantify their degree of dependence. An essential feature of our approach is the ability to track the time-varying association between the populations while making minimal assumptions on their functional shapes via Markov random field priors. We provide nonparametric estimators, extensions of our base model that integrate multiple data sources, and fast scalable inference algorithms. We test our method using simulated data under various dependent population histories and demonstrate the utility of our model in shedding light on evolutionary histories of different variants of SARS-CoV-2.

Phased DNA methylation states within bisulfite sequencing reads are valuable source of information that can be used to estimate epigenetic diversity across cells as well as epigenomic instability in individual cells. Various measures capturing the heterogeneity of DNA methylation states have been proposed for a decade. However, in routine analyses on DNA methylation, this heterogeneity is often ignored by computing average methylation levels at CpG sites, even though such information exists in bisulfite sequencing data in the form of phased methylation states, or methylation patterns. In this study, to facilitate the application of the DNA methylation heterogeneity measures in downstream epigenomic analyses, we present a Rust-based, extremely fast and lightweight bioinformatics toolkit called Metheor. As the analysis of DNA methylation heterogeneity requires the examination of pairs or groups of CpGs throughout the genome, existing softwares suffer from high computational burden, which almost make a large-scale DNA methylation heterogeneity studies intractable for researchers with limited resources. In this study, we benchmark the performance of Metheor against existing code implementations for DNA methylation heterogeneity measures in three different scenarios of simulated bisulfite sequencing datasets. Metheor was shown to dramatically reduce the execution time up to 300-fold and memory footprint up to 60-fold, while producing identical results with the original implementation, thereby facilitating a large-scale study of DNA methylation heterogeneity profiles. To demonstrate the utility of the low computational burden of Metheor, we show that the methylation heterogeneity profiles of 928 cancer cell lines can be computed with standard computing resources. With those profiles, we reveal the association between DNA methylation heterogeneity and various omics features. Source code for Metheor is at https://github.com/dohlee/metheor and is freely available under the GPL-3.0 license.

As evidenced by the worldwide pandemic, respiratory infectious diseases and their airborne transmission must be studied to safeguard public health. This study focuses on the emission and transport of speech-generated droplets, which can pose risk of infection depending on the loudness of the speech, its duration and the initial angle of exhalation. We have numerically investigated the transport of these droplets into the human respiratory tract by way of a natural breathing cycle in order to predict the infection probability of three strains of SARS-CoV-2 on a person who is listening at a one-meter distance. Numerical methods were used to set the boundary conditions of the speaking and breathing models and large eddy simulation (LES) was used for the unsteady simulation of approximately 10 breathing cycles. Four different mouth angles when speaking were contrasted to evaluate real conditions of human communication and the possibility of infection. Breathed virions were counted using two different approaches: the breathing zone of influence and direction deposition on the tissue. Our results show that infection probability drastically changes based on the mouth angle and the breathing zone of influence overpredicts the inhalation risk in all cases. We conclude that to portray real conditions, the probability of infection should be based on direct tissue deposition results to avoid overprediction and that several mouth angles must be considered in future analyses.

A perfect bacterial genome assembly is one where the assembled sequence is an exact match for the organism's genome-each replicon sequence is complete and contains no errors. While this has been difficult to achieve in the past, improvements in long-read sequencing, assemblers, and polishers have brought perfect assemblies within reach. Here, we describe our recommended approach for assembling a bacterial genome to perfection using a combination of Oxford Nanopore Technologies long reads and Illumina short reads: Trycycler long-read assembly, Medaka long-read polishing, Polypolish short-read polishing, followed by other short-read polishing tools and manual curation. We also discuss potential pitfalls one might encounter when assembling challenging genomes, and we provide an online tutorial with sample data (github.com/rrwick/perfect-bacterial-genome-tutorial).

The movement of groups can be heavily influenced by 'leader' individuals who differ from the others in some way. A major source of differences between individuals is the repeatability and consistency of their behaviour, commonly considered as their 'personality', which can influence both position within a group as well as the tendency to lead. However, links between personality and behaviour may also depend upon the immediate social environment of the individual; individuals who behave consistently in one way when alone may not express the same behaviour socially, when they may be conforming with the behaviour of others. Experimental evidence shows that personality differences can be eroded in social situations, but there is currently a lack of theory to identify the conditions where we would expect personality to be suppressed. Here, we develop a simple individual-based framework considering a small group of individuals with differing tendencies to perform risky behaviours when travelling away from a safe home site towards a foraging site, and compare the group behaviours when the individuals follow differing rules for aggregation behaviour determining how much attention they pay to the actions of their fellow group-members. We find that if individuals pay attention to the other members of the group, the group will tend to remain at the safe site for longer, but then travel faster towards the foraging site. This demonstrates that simple social behaviours can result in the repression of consistent inter-individual differences in behaviour, giving the first theoretical consideration of the social mechanisms behind personality suppression.

The demand for well-trained bioinformaticians to support genomics research continues to rise. Unfortunately, undergraduate training in Kenya does not prepare students for specialization in bioinformatics. Graduates are often unaware of the career opportunities in bioinformatics, and those who are may lack mentors to help them choose a specialization. The Bioinformatics Mentorship and Incubation Program seeks to bridge the gap by laying the foundation for a bioinformatics training pipeline using project-based learning. The program selects six participants through an intensive open recruitment exercise for highly competitive students to join the program for four months. The six interns undergo intensive training within the first one and a half months before being assigned to mini-projects. We track the progress of the interns weekly through code review sessions and a final presentation at the end of the four months. We have trained five cohorts, most of whom have secured master's scholarships within and outside the country and job opportunities. We demonstrate the benefit of structured mentorship using project-based learning in filling the training gap after undergraduate programs to generate well-trained bioinformaticians who are competitive in graduate programs and bioinformatics jobs.

In mass spectrometry (MS)-based proteomics, protein inference from identified peptides (protein fragments) is a critical step. We present ProInfer (Protein Inference), a novel protein assembly method that takes advantage of information in biological networks. ProInfer assists recovery of proteins supported only by ambiguous peptides (a peptide which maps to more than one candidate protein) and enhances the statistical confidence for proteins supported by both unique and ambiguous peptides. Consequently, ProInfer rescues weakly supported proteins thereby improving proteome coverage. Evaluated across THP1 cell line, lung cancer and RAW267.4 datasets, ProInfer always infers the most numbers of true positives, in comparison to mainstream protein inference tools Fido, EPIFANY and PIA. ProInfer is also adept at retrieving differentially expressed proteins, signifying its usefulness for functional analysis and phenotype profiling. Source codes of ProInfer are available at https://github.com/PennHui2016/ProInfer.

During breast cancer metastasis, the developmental process epithelial-mesenchymal (EM) transition is abnormally activated. Transcriptional regulatory networks controlling EM transition are well-studied; however, alternative RNA splicing also plays a critical regulatory role during this process. Alternative splicing was proved to control the EM transition process, and RNA-binding proteins were determined to regulate alternative splicing. A comprehensive understanding of alternative splicing and the RNA-binding proteins that regulate it during EM transition and their dynamic impact on breast cancer remains largely unknown. To accurately study the dynamic regulatory relationships, time-series data of the EM transition process are essential. However, only cross-sectional data of epithelial and mesenchymal specimens are available. Therefore, we developed a pseudotemporal causality-based Bayesian (PCB) approach to infer the dynamic regulatory relationships between alternative splicing events and RNA-binding proteins. Our study sheds light on facilitating the regulatory network-based approach to identify key RNA-binding proteins or target alternative splicing events for the diagnosis or treatment of cancers. The data and code for PCB are available at: http://hkumath.hku.hk/~wkc/PCB(data+code).zip.

[This corrects the article DOI: 10.1371/journal.pcbi.1009813.].

Conducting a replication study is a valuable way for undergraduate students to learn about the scientific process and gain research experience. By promoting the evaluation of existing studies to confirm their reliability, replications play a unique, though often underappreciated, role in the scientific enterprise. Involving students early in this process can help make replication mainstream among the new generation of scientists. Beyond their benefit to science, replications also provide an invaluable learning ground for students, from encouraging the development of critical thinking to emphasizing the importance of details and honing research skills. In this piece, we outline 10 simple rules for designing and conducting undergraduate replication projects, from conceptualization to implementation and dissemination. We hope that these guidelines can help educators provide students with a meaningful and constructive pedagogical experience, without compromising the scientific value of the replication project, therefore ensuring robust, valuable contributions to our understanding of the world.

Computational models of infectious diseases have become valuable tools for research and the public health response against epidemic threats. The reproducibility of computational models has been limited, undermining the scientific process and possibly trust in modeling results and related response strategies, such as vaccination. We translated published reproducibility guidelines from a wide range of scientific disciplines into an implementation framework for improving reproducibility of infectious disease computational models. The framework comprises 22 elements that should be described, grouped into 6 categories: computational environment, analytical software, model description, model implementation, data, and experimental protocol. The framework can be used by scientific communities to develop actionable tools for sharing computational models in a reproducible way.

Estimating feature importance, which is the contribution of a prediction or several predictions due to a feature, is an essential aspect of explaining data-based models. Besides explaining the model itself, an equally relevant question is which features are important in the underlying data generating process. We present a Shapley-value-based framework for inferring the importance of individual features, including uncertainty in the estimator. We build upon the recently published model-agnostic feature importance score of SAGE (Shapley additive global importance) and introduce Sub-SAGE. For tree-based models, it has the advantage that it can be estimated without computationally expensive resampling. We argue that for all model types the uncertainties in our Sub-SAGE estimator can be estimated using bootstrapping and demonstrate the approach for tree ensemble methods. The framework is exemplified on synthetic data as well as large genotype data for predicting feature importance with respect to obesity.

We introduce a self-describing serialized format for bulk biomedical data called the Portable Format for Biomedical (PFB) data. The Portable Format for Biomedical data is based upon Avro and encapsulates a data model, a data dictionary, the data itself, and pointers to third party controlled vocabularies. In general, each data element in the data dictionary is associated with a third party controlled vocabulary to make it easier for applications to harmonize two or more PFB files. We also introduce an open source software development kit (SDK) called PyPFB for creating, exploring and modifying PFB files. We describe experimental studies showing the performance improvements when importing and exporting bulk biomedical data in the PFB format versus using JSON and SQL formats.

Deaths by suicide, as well as suicidal ideations, plans and attempts, have been increasing in the US for the past two decades. Deployment of effective interventions would require timely, geographically well-resolved estimates of suicide activity. In this study, we evaluated the feasibility of a two-step process for predicting suicide mortality: a) generation of hindcasts, mortality estimates for past months for which observational data would not have been available if forecasts were generated in real-time; and b) generation of forecasts with observational data augmented with hindcasts. Calls to crisis hotline services and online queries to the Google search engine for suicide-related terms were used as proxy data sources to generate hindcasts. The primary hindcast model (auto) is an Autoregressive Integrated Moving average model (ARIMA), trained on suicide mortality rates alone. Three regression models augment hindcast estimates from auto with call rates (calls), GHT search rates (ght) and both datasets together (calls_ght). The 4 forecast models used are ARIMA models trained with corresponding hindcast estimates. All models were evaluated against a baseline random walk with drift model. Rolling monthly 6-month ahead forecasts for all 50 states between 2012 and 2020 were generated. Quantile score (QS) was used to assess the quality of the forecast distributions. Median QS for auto was better than baseline (0.114 vs. 0.21. Median QS of augmented models were lower than auto, but not significantly different from each other (Wilcoxon signed-rank test, p > .05). Forecasts from augmented models were also better calibrated. Together, these results provide evidence that proxy data can address delays in release of suicide mortality data and improve forecast quality. An operational forecast system of state-level suicide risk may be feasible with sustained engagement between modelers and public health departments to appraise data sources and methods as well as to continuously evaluate forecast accuracy.

BACKGROUND: Pneumonia remains a leading cause of hospitalization and death among young children worldwide, and the diagnostic challenge of differentiating bacterial from non-bacterial pneumonia is the main driver of antibiotic use for treating pneumonia in children. Causal Bayesian networks (BNs) serve as powerful tools for this problem as they provide clear maps of probabilistic relationships between variables and produce results in an explainable way by incorporating both domain expert knowledge and numerical data. METHODS: We used domain expert knowledge and data in combination and iteratively, to construct, parameterise and validate a causal BN to predict causative pathogens for childhood pneumonia. Expert knowledge elicitation occurred through a series of group workshops, surveys and one-on-one meetings involving 6-8 experts from diverse domain areas. The model performance was evaluated based on both quantitative metrics and qualitative expert validation. Sensitivity analyses were conducted to investigate how the target output is influenced by varying key assumptions of a particularly high degree of uncertainty around data or domain expert knowledge. RESULTS: Designed to apply to a cohort of children with X-ray confirmed pneumonia who presented to a tertiary paediatric hospital in Australia, the resulting BN offers explainable and quantitative predictions on a range of variables of interest, including the diagnosis of bacterial pneumonia, detection of respiratory pathogens in the nasopharynx, and the clinical phenotype of a pneumonia episode. Satisfactory numeric performance has been achieved including an area under the receiver operating characteristic curve of 0.8 in predicting clinically-confirmed bacterial pneumonia with sensitivity 88% and specificity 66% given certain input scenarios (i.e., information that is available and entered into the model) and trade-off preferences (i.e., relative weightings of the consequences of false positive versus false negative predictions). We specifically highlight that a desirable model output threshold for practical use is very dependent upon different input scenarios and trade-off preferences. Three commonly encountered scenarios were presented to demonstrate the potential usefulness of the BN outputs in various clinical pictures. CONCLUSIONS: To our knowledge, this is the first causal model developed to help determine the causative pathogen for paediatric pneumonia. We have shown how the method works and how it would help decision making on the use of antibiotics, providing insight into how computational model predictions may be translated to actionable decisions in practice. We discussed key next steps including external validation, adaptation and implementation. Our model framework and the methodological approach can be adapted beyond our context to broad respiratory infections and geographical and healthcare settings.

The signature of early cancer dynamics on the spatial arrangement of tumour cells is poorly understood, and yet could encode information about how sub-clones grew within the expanding tumour. Novel methods of quantifying spatial tumour data at the cellular scale are required to link evolutionary dynamics to the resulting spatial architecture of the tumour. Here, we propose a framework using first passage times of random walks to quantify the complex spatial patterns of tumour cell population mixing. First, using a simple model of cell mixing we demonstrate how first passage time statistics can distinguish between different pattern structures. We then apply our method to simulated patterns of mutated and non-mutated tumour cell population mixing, generated using an agent-based model of expanding tumours, to explore how first passage times reflect mutant cell replicative advantage, time of emergence and strength of cell pushing. Finally, we explore applications to experimentally measured human colorectal cancer, and estimate parameters of early sub-clonal dynamics using our spatial computational model. We infer a wide range of sub-clonal dynamics, with mutant cell division rates varying between 1 and 4 times the rate of non-mutated cells across our sample set. Some mutated sub-clones emerged after as few as 100 non-mutant cell divisions, and others only after 50,000 divisions. The majority were consistent with boundary driven growth or short-range cell pushing. By analysing multiple sub-sampled regions in a small number of samples, we explore how the distribution of inferred dynamics could inform about the initial mutational event. Our results demonstrate the efficacy of first passage time analysis as a new methodology in spatial analysis of solid tumour tissue, and suggest that patterns of sub-clonal mixing can provide insights into early cancer dynamics.

G-quadruplexes are non-B-DNA structures that form in the genome facilitated by Hoogsteen bonds between guanines in single or multiple strands of DNA. The functions of G-quadruplexes are linked to various molecular and disease phenotypes, and thus researchers are interested in measuring G-quadruplex formation genome-wide. Experimentally measuring G-quadruplexes is a long and laborious process. Computational prediction of G-quadruplex propensity from a given DNA sequence is thus a long-standing challenge. Unfortunately, despite the availability of high-throughput datasets measuring G-quadruplex propensity in the form of mismatch scores, extant methods to predict G-quadruplex formation either rely on small datasets or are based on domain-knowledge rules. We developed G4mismatch, a novel algorithm to accurately and efficiently predict G-quadruplex propensity for any genomic sequence. G4mismatch is based on a convolutional neural network trained on almost 400 millions human genomic loci measured in a single G4-seq experiment. When tested on sequences from a held-out chromosome, G4mismatch, the first method to predict mismatch scores genome-wide, achieved a Pearson correlation of over 0.8. When benchmarked on independent datasets derived from various animal species, G4mismatch trained on human data predicted G-quadruplex propensity genome-wide with high accuracy (Pearson correlations greater than 0.7). Moreover, when tested in detecting G-quadruplexes genome-wide using the predicted mismatch scores, G4mismatch achieved superior performance compared to extant methods. Last, we demonstrate the ability to deduce the mechanism behind G-quadruplex formation by unique visualization of the principles learned by the model.

In many real-world applications, such as those based on electronic health records, prognostic prediction of patient survival is based on heterogeneous sets of clinical laboratory measurements. To address the trade-off between the predictive accuracy of a prognostic model and the costs related to its clinical implementation, we propose an optimized L0-pseudonorm approach to learn sparse solutions in multivariable regression. The model sparsity is maintained by restricting the number of nonzero coefficients in the model with a cardinality constraint, which makes the optimization problem NP-hard. In addition, we generalize the cardinality constraint for grouped feature selection, which makes it possible to identify key sets of predictors that may be measured together in a kit in clinical practice. We demonstrate the operation of our cardinality constraint-based feature subset selection method, named OSCAR, in the context of prognostic prediction of prostate cancer patients, where it enables one to determine the key explanatory predictors at different levels of model sparsity. We further explore how the model sparsity affects the model accuracy and implementation cost. Lastly, we demonstrate generalization of the presented methodology to high-dimensional transcriptomics data.

A quantitative understanding of the dynamics of bee colonies is important to support global efforts to improve bee health and enhance pollination services. Traditional approaches focus either on theoretical models or data-centred statistical analyses. Here we argue that the combination of these two approaches is essential to obtain interpretable information on the state of bee colonies and show how this can be achieved in the case of time series of intra-day weight variation. We model how the foraging and food processing activities of bees affect global hive weight through a set of ordinary differential equations and show how to estimate the parameters of this model from measurements on a single day. Our analysis of 10 hives at different times shows that the estimation of crucial indicators of the health of honey bee colonies are statistically reliable and fall in ranges compatible with previously reported results. The crucial indicators, which include the amount of food collected (foraging success) and the number of active foragers, may be used to develop early warning indicators of colony failure.

Directed laboratory evolution applies iterative rounds of mutation and selection to explore the protein fitness landscape and provides rich information regarding the underlying relationships between protein sequence, structure, and function. Laboratory evolution data consist of protein sequences sampled from evolving populations over multiple generations and this data type does not fit into established supervised and unsupervised machine learning approaches. We develop a statistical learning framework that models the evolutionary process and can infer the protein fitness landscape from multiple snapshots along an evolutionary trajectory. We apply our modeling approach to dihydrofolate reductase (DHFR) laboratory evolution data and the resulting landscape parameters capture important aspects of DHFR structure and function. We use the resulting model to understand the structure of the fitness landscape and find numerous examples of epistasis but an overall global peak that is evolutionarily accessible from most starting sequences. Finally, we use the model to perform an in silico extrapolation of the DHFR laboratory evolution trajectory and computationally design proteins from future evolutionary rounds.