Computational modeling has been applied for data analysis in psychology, neuroscience, and psychiatry. One of its important uses is to infer the latent variables underlying behavior by which researchers can evaluate corresponding neural, physiological, or behavioral measures. This feature is especially crucial for computational psychiatry, in which altered computational processes underlying mental disorders are of interest. For instance, several studies employing model-based fMRI-a method for identifying brain regions correlated with latent variables-have shown that patients with mental disorders (e.g., depression) exhibit diminished neural responses to reward prediction errors (RPEs), which are the differences between experienced and predicted rewards. Such model-based analysis has the drawback that the parameter estimates and inference of latent variables are not necessarily correct-rather, they usually contain some errors. A previous study theoretically and empirically showed that the error in model-fitting does not necessarily cause a serious error in model-based fMRI. However, the study did not deal with certain situations relevant to psychiatry, such as group comparisons between patients and healthy controls. We developed a theoretical framework to explore such situations. We demonstrate that the parameter-misspecification can critically affect the results of group comparison. We demonstrate that even if the RPE response in patients is completely intact, a spurious difference to healthy controls is observable. Such a situation occurs when the ground-truth learning rate differs between groups but a common learning rate is used, as per previous studies. Furthermore, even if the parameters are appropriately fitted to individual participants, spurious group differences in RPE responses are observable when the model lacks a component that differs between groups. These results highlight the importance of appropriate model-fitting and the need for caution when interpreting the results of model-based fMRI.

Hepatitis B virus (HBV) infection contributes to hepatocellular carcinoma (HCC) initiation and is associated with worse outcomes. Many prior studies of HBV-related HCC have not accounted for potential heterogeneity among HBV-related tumors by assessing whether HBV activity is present in tumor tissue. Here, we measured tumor HBV RNA, a proxy for viral activity, and investigated the association between HBV RNA status and several clinicogenomic characteristics. We obtained clinical, mutation, RNA-Seq and survival data for 439 HCC tumors from The Cancer Genome Atlas and International Cancer Genome Consortium. Tumors were classified as HBV RNA positive if they harbored >1 HBV RNA read per million human reads. We investigated the association between HBV RNA status and nonsynonymous somatic mutations, gene set expression, homologous recombination deficiency (HRD) score and mutation-specific survival. HBV RNA positive status was associated with higher nonsynonymous mutation rates of multiple genes, including TP53 and CDKN2A, while HBV RNA negative status was associated with higher nonsynonymous BAP1 mutation rate. HBV RNA positive status was also associated with increased transcription of genes involved in multiple DNA damage repair pathways, genes upregulated by MYC and mTORC1, and genes overexpressed in several HCC subclasses associated with a proliferative phenotype. Further, HBV RNA positive status was associated with increased three-biomarker HRD score (22.2 for HBV RNA+ vs. 16.0 for HBV RNA-). Finally, HBV RNA status was associated with multiple mutation-specific survival differences, including decreased survival for HBV RNA positive patients with nonsynonymous KEAP1 mutations compared to those without (hazard ratio 4.26). HCC tumors harboring genomic evidence of HBV activity therefore constitute a distinct HCC subset characterized by specific differences in nonsynonymous mutations, gene set expression, three-biomarker HRD score and mutation-specific survival.

Hepatocellular carcinoma (HCC) is a significant health problem worldwide with poor prognosis. Drug repositioning represents a profitable strategy to accelerate drug discovery in the treatment of HCC. In this study, we developed a new approach for predicting therapeutic drugs for HCC based on tissue-specific pathways and identified three newly predicted drugs that are likely to be therapeutic drugs for the treatment of HCC. We validated these predicted drugs by analyzing their overlapping drug indications reported in PubMed literature. By using the cancer cell line data in the database, we constructed a Connectivity Map (CMap) profile similarity analysis and KEGG enrichment analysis on their related genes. By experimental validation, we found securinine and ajmaline significantly inhibited cell viability of HCC cells and induced apoptosis. Among them, securinine has lower toxicity to normal liver cell line, which is worthy of further research. Our results suggested that the proposed approach was effective and accurate for discovering novel therapeutic options for HCC. This method also could be used to indicate unmarked drug-disease associations in the Comparative Toxicogenomics Database. Meanwhile, our method could also be applied to predict the potential drugs for other types of tumors by changing the database.

Metagenomics has redefined many areas of microbiology. However, metagenome-assembled genomes (MAGs) are often fragmented, primarily when sequencing was performed with short reads. Recent long-read sequencing technologies promise to improve genome reconstruction. However, the integration of two different sequencing modalities makes downstream analyses complex. We, therefore, developed MUFFIN, a complete metagenomic workflow that uses short and long reads to produce high-quality bins and their annotations. The workflow is written by using Nextflow, a workflow orchestration software, to achieve high reproducibility and fast and straightforward use. This workflow also produces the taxonomic classification and KEGG pathways of the bins and can be further used for quantification and annotation by providing RNA-Seq data (optionally). We tested the workflow using twenty biogas reactor samples and assessed the capacity of MUFFIN to process and output relevant files needed to analyze the microbial community and their function. MUFFIN produces functional pathway predictions and, if provided de novo metatranscript annotations across the metagenomic sample and for each bin. MUFFIN is available on github under GNUv3 licence: https://github.com/RVanDamme/MUFFIN.

Traveling waves are commonly observed across the brain. While previous studies have suggested the role of traveling waves in learning, the mechanism remains unclear. We adopted a computational approach to investigate the effect of traveling waves on synaptic plasticity. Our results indicate that traveling waves facilitate the learning of poly-synaptic network paths when combined with a reward-dependent local synaptic plasticity rule. We also demonstrate that traveling waves expedite finding the shortest paths and learning nonlinear input/output mapping, such as exclusive or (XOR) function.

Many ecological studies employ general models that can feature an arbitrary number of populations. A critical requirement imposed on such models is clone consistency: If the individuals from two populations are indistinguishable, joining these populations into one shall not affect the outcome of the model. Otherwise a model produces different outcomes for the same scenario. Using functional analysis, we comprehensively characterize all clone-consistent models: We prove that they are necessarily composed from basic building blocks, namely linear combinations of parameters and abundances. These strong constraints enable a straightforward validation of model consistency. Although clone consistency can always be achieved with sufficient assumptions, we argue that it is important to explicitly name and consider the assumptions made: They may not be justified or limit the applicability of models and the generality of the results obtained with them. Moreover, our insights facilitate building new clone-consistent models, which we illustrate for a data-driven model of microbial communities. Finally, our insights point to new relevant forms of general models for theoretical ecology. Our framework thus provides a systematic way of comprehending ecological models, which can guide a wide range of studies.

Axonal connections are widely regarded as faithful transmitters of neuronal signals with fixed delays. The reasoning behind this is that extracellular potentials caused by spikes travelling along axons are too small to have an effect on other axons. Here we devise a computational framework that allows us to study the effect of extracellular potentials generated by spike volleys in axonal fibre bundles on axonal transmission delays. We demonstrate that, although the extracellular potentials generated by single spikes are of the order of microvolts, the collective extracellular potential generated by spike volleys can reach several millivolts. As a consequence, the resulting depolarisation of the axonal membranes increases the velocity of spikes, and therefore reduces axonal delays between brain areas. Driving a neural mass model with such spike volleys, we further demonstrate that only ephaptic coupling can explain the reduction of stimulus latencies with increased stimulus intensities, as observed in many psychological experiments.

Mechano-regulation during tendon healing, i.e. the relationship between mechanical stimuli and cellular response, has received more attention recently. However, the basic mechanobiological mechanisms governing tendon healing after a rupture are still not well-understood. Literature has reported spatial and temporal variations in the healing of ruptured tendon tissue. In this study, we explored a computational modeling approach to describe tendon healing. In particular, a novel 3D mechano-regulatory framework was developed to investigate spatio-temporal evolution of collagen content and orientation, and temporal evolution of tendon stiffness during early tendon healing. Based on an extensive literature search, two possible relationships were proposed to connect levels of mechanical stimuli to collagen production. Since literature remains unclear on strain-dependent collagen production at high levels of strain, the two investigated production laws explored the presence or absence of collagen production upon non-physiologically high levels of strain (>15%). Implementation in a finite element framework, pointed to large spatial variations in strain magnitudes within the callus tissue, which resulted in predictions of distinct spatial distributions of collagen over time. The simulations showed that the magnitude of strain was highest in the tendon core along the central axis, and decreased towards the outer periphery. Consequently, decreased levels of collagen production for high levels of tensile strain were shown to accurately predict the experimentally observed delayed collagen production in the tendon core. In addition, our healing framework predicted evolution of collagen orientation towards alignment with the tendon axis and the overall predicted tendon stiffness agreed well with experimental data. In this study, we explored the capability of a numerical model to describe spatial and temporal variations in tendon healing and we identified that understanding mechano-regulated collagen production can play a key role in explaining heterogeneities observed during tendon healing.

As cooperation incurs a cost to the cooperator for others to benefit, its evolution seems to contradict natural selection. How evolution has resolved this obstacle has been among the most intensely studied questions in evolutionary theory in recent decades. Here, we show that having a choice between different public resources provides a simple mechanism for cooperation to flourish. Such a mechanism can be at work in many biological or social contexts where individuals can form different groups or join different institutions to perform a collective action task, or when they can choose between collective actions with different profitability. As a simple evolutionary model suggests, defectors tend to join the highest quality resource in such a context. This allows cooperators to survive and out-compete defectors by sheltering in a lower quality resource. Cooperation is maximized, however, when the qualities of the two highest quality resources are similar, and thus, they are almost interchangeable.

Tyrosine and serine/threonine kinases are essential regulators of cell processes and are important targets for human therapies. Unfortunately, very little is known about specific kinase-substrate relationships, making it difficult to infer meaning from dysregulated phosphoproteomic datasets or for researchers to identify possible kinases that regulate specific or novel phosphorylation sites. The last two decades have seen an explosion in algorithms to extrapolate from what little is known into the larger unknown-predicting kinase relationships with site-specific substrates using a variety of approaches that include the sequence-specificity of kinase catalytic domains and various other factors, such as evolutionary relationships, co-expression, and protein-protein interaction networks. Unfortunately, a number of limitations prevent researchers from easily harnessing these resources, such as loss of resource accessibility, limited information in publishing that results in a poor mapping to a human reference, and not being updated to match the growth of the human phosphoproteome. Here, we propose a methodological framework for publishing predictions in a unified way, which entails ensuring predictions have been run on a current reference proteome, mapping the same substrates and kinases across resources to a common reference, filtering for the human phosphoproteome, and providing methods for updating the resource easily in the future. We applied this framework on three currently available resources, published in the last decade, which provide kinase-specific predictions in the human proteome. Using the unified datasets, we then explore the role of study bias, the emergent network properties of these predictive algorithms, and comparisons within and between predictive algorithms. The combination of the code for unification and analysis, as well as the unified predictions are available under the resource we named KinPred. We believe this resource will be useful for a wide range of applications and establishes best practices for long-term usability and sustainability for new and existing predictive algorithms.

Knowledge of the genetic basis for autotrophic metabolism is valuable since it relates to both the emergence of life and to the metabolic engineering challenge of incorporating CO2 as a potential substrate for biorefining. The most common CO2 fixation pathway is the Calvin cycle, which utilizes Rubisco and phosphoribulokinase enzymes. We searched thousands of microbial genomes and found that 6.0% contained the Calvin cycle. We then contrasted the genomes of Calvin cycle-positive, non-cyanobacterial microbes and their closest relatives by enrichment analysis, ancestral character estimation, and random forest machine learning, to explore genetic adaptations associated with acquisition of the Calvin cycle. The Calvin cycle overlaps with the pentose phosphate pathway and glycolysis, and we could confirm positive associations with fructose-1,6-bisphosphatase, aldolase, and transketolase, constituting a conserved operon, as well as ribulose-phosphate 3-epimerase, ribose-5-phosphate isomerase, and phosphoglycerate kinase. Additionally, carbohydrate storage enzymes, carboxysome proteins (that raise CO2 concentration around Rubisco), and Rubisco activases CbbQ and CbbX accompanied the Calvin cycle. Photorespiration did not appear to be adapted specifically for the Calvin cycle in the non-cyanobacterial microbes under study. Our results suggest that chemoautotrophy in Calvin cycle-positive organisms was commonly enabled by hydrogenase, and less commonly ammonia monooxygenase (nitrification). The enrichment of specific DNA-binding domains indicated Calvin-cycle associated genetic regulation. Metabolic regulatory adaptations were illustrated by negative correlation to AraC and the enzyme arabinose-5-phosphate isomerase, which suggests a downregulation of the metabolite arabinose-5-phosphate, which may interfere with the Calvin cycle through enzyme inhibition and substrate competition. Certain domains of unknown function that were found to be important in the analysis may indicate yet unknown regulatory mechanisms in Calvin cycle-utilizing microbes. Our gene ranking provides targets for experiments seeking to improve CO2 fixation, or engineer novel CO2-fixing organisms.

There is an emerging consensus that achieving global tuberculosis control targets will require more proactive case finding approaches than are currently used in high-incidence settings. Household contact tracing (HHCT), for which households of newly diagnosed cases are actively screened for additional infected individuals is a potentially efficient approach to finding new cases of tuberculosis, however randomized trials assessing the population-level effects of such interventions in settings with sustained community transmission have shown mixed results. One potential explanation for this is that household transmission is responsible for a variable proportion of population-level tuberculosis burden between settings. For example, transmission is more likely to occur in households in settings with a lower tuberculosis burden and where individuals mix preferentially in local areas, compared with settings with higher disease burden and more dispersed mixing. To better understand the relationship between endemic incidence levels, social mixing, and the impact of HHCT, we developed a spatially explicit model of coupled household and community transmission. We found that the impact of HHCT was robust across settings of varied incidence and community contact patterns. In contrast, we found that the effects of community contact tracing interventions were sensitive to community contact patterns. Our results suggest that the protective benefits of HHCT are robust and the benefits of this intervention are likely to be maintained across epidemiological settings.

The stimulation of a single neuron in the rat somatosensory cortex can elicit a behavioral response. The probability of a behavioral response does not depend appreciably on the duration or intensity of a constant stimulation, whereas the response probability increases significantly upon injection of an irregular current. Biological mechanisms that can potentially suppress a constant input signal are present in the dynamics of both neurons and synapses and seem ideal candidates to explain these experimental findings. Here, we study a large network of integrate-and-fire neurons with several salient features of neuronal populations in the rat barrel cortex. The model includes cellular spike-frequency adaptation, experimentally constrained numbers and types of chemical synapses endowed with short-term plasticity, and gap junctions. Numerical simulations of this model indicate that cellular and synaptic adaptation mechanisms alone may not suffice to account for the experimental results if the local network activity is read out by an integrator. However, a circuit that approximates a differentiator can detect the single-cell stimulation with a reliability that barely depends on the length or intensity of the stimulus, but that increases when an irregular signal is used. This finding is in accordance with the experimental results obtained for the stimulation of a regularly-spiking excitatory cell.

Since the seminal 1961 paper of Monod and Jacob, mathematical models of biomolecular circuits have guided our understanding of cell regulation. Model-based exploration of the functional capabilities of any given circuit requires systematic mapping of multidimensional spaces of model parameters. Despite significant advances in computational dynamical systems approaches, this analysis remains a nontrivial task. Here, we use a nonlinear system of ordinary differential equations to model oocyte selection in Drosophila, a robust symmetry-breaking event that relies on autoregulatory localization of oocyte-specification factors. By applying an algorithmic approach that implements symbolic computation and topological methods, we enumerate all phase portraits of stable steady states in the limit when nonlinear regulatory interactions become discrete switches. Leveraging this initial exact partitioning and further using numerical exploration, we locate parameter regions that are dense in purely asymmetric steady states when the nonlinearities are not infinitely sharp, enabling systematic identification of parameter regions that correspond to robust oocyte selection. This framework can be generalized to map the full parameter spaces in a broad class of models involving biological switches.

The novelty of new human coronavirus COVID-19/SARS-CoV-2 and the lack of effective drugs and vaccines gave rise to a wide variety of strategies employed to fight this worldwide pandemic. Many of these strategies rely on the repositioning of existing drugs that could shorten the time and reduce the cost compared to de novo drug discovery. In this study, we presented a new network-based algorithm for drug repositioning, called SAveRUNNER (Searching off-lAbel dRUg aNd NEtwoRk), which predicts drug-disease associations by quantifying the interplay between the drug targets and the disease-specific proteins in the human interactome via a novel network-based similarity measure that prioritizes associations between drugs and diseases locating in the same network neighborhoods. Specifically, we applied SAveRUNNER on a panel of 14 selected diseases with a consolidated knowledge about their disease-causing genes and that have been found to be related to COVID-19 for genetic similarity (i.e., SARS), comorbidity (e.g., cardiovascular diseases), or for their association to drugs tentatively repurposed to treat COVID-19 (e.g., malaria, HIV, rheumatoid arthritis). Focusing specifically on SARS subnetwork, we identified 282 repurposable drugs, including some the most rumored off-label drugs for COVID-19 treatments (e.g., chloroquine, hydroxychloroquine, tocilizumab, heparin), as well as a new combination therapy of 5 drugs (hydroxychloroquine, chloroquine, lopinavir, ritonavir, remdesivir), actually used in clinical practice. Furthermore, to maximize the efficiency of putative downstream validation experiments, we prioritized 24 potential anti-SARS-CoV repurposable drugs based on their network-based similarity values. These top-ranked drugs include ACE-inhibitors, monoclonal antibodies (e.g., anti-IFNγ, anti-TNFα, anti-IL12, anti-IL1β, anti-IL6), and thrombin inhibitors. Finally, our findings were in-silico validated by performing a gene set enrichment analysis, which confirmed that most of the network-predicted repurposable drugs may have a potential treatment effect against human coronavirus infections.

For many prevalent complex diseases, treatment regimens are frequently ineffective. For example, despite multiple available immunomodulators and immunosuppressants, inflammatory bowel disease (IBD) remains difficult to treat. Heterogeneity in the disease across patients makes it challenging to select the optimal treatment regimens, and some patients do not respond to any of the existing treatment choices. Drug repurposing strategies for IBD have had limited clinical success and have not typically offered individualized patient-level treatment recommendations. In this work, we present NetPTP, a Network-based Personalized Treatment Prediction framework which models measured drug effects from gene expression data and applies them to patient samples to generate personalized ranked treatment lists. To accomplish this, we combine publicly available network, drug target, and drug effect data to generate treatment rankings using patient data. These ranked lists can then be used to prioritize existing treatments and discover new therapies for individual patients. We demonstrate how NetPTP captures and models drug effects, and we apply our framework to individual IBD samples to provide novel insights into IBD treatment.

Neural activity is often low dimensional and dominated by only a few prominent neural covariation patterns. It has been hypothesised that these covariation patterns could form the building blocks used for fast and flexible motor control. Supporting this idea, recent experiments have shown that monkeys can learn to adapt their neural activity in motor cortex on a timescale of minutes, given that the change lies within the original low-dimensional subspace, also called neural manifold. However, the neural mechanism underlying this within-manifold adaptation remains unknown. Here, we show in a computational model that modification of recurrent weights, driven by a learned feedback signal, can account for the observed behavioural difference between within- and outside-manifold learning. Our findings give a new perspective, showing that recurrent weight changes do not necessarily lead to change in the neural manifold. On the contrary, successful learning is naturally constrained to a common subspace.

The visual cortex of the mouse brain can be divided into ten or more areas that each contain complete or partial retinotopic maps of the contralateral visual field. It is generally assumed that these areas represent discrete processing regions. In contrast to the conventional input-output characterizations of neuronal responses to standard visual stimuli, here we asked whether six of the core visual areas have responses that are functionally distinct from each other for a given visual stimulus set, by applying machine learning techniques to distinguish the areas based on their activity patterns. Visual areas defined by retinotopic mapping were examined using supervised classifiers applied to responses elicited by a range of stimuli. Using two distinct datasets obtained using wide-field and two-photon imaging, we show that the area labels predicted by the classifiers were highly consistent with the labels obtained using retinotopy. Furthermore, the classifiers were able to model the boundaries of visual areas using resting state cortical responses obtained without any overt stimulus, in both datasets. With the wide-field dataset, clustering neuronal responses using a constrained semi-supervised classifier showed graceful degradation of accuracy. The results suggest that responses from visual cortical areas can be classified effectively using data-driven models. These responses likely reflect unique circuits within each area that give rise to activity with stronger intra-areal than inter-areal correlations, and their responses to controlled visual stimuli across trials drive higher areal classification accuracy than resting state responses.

During developmental angiogenesis, endothelial cells respond to shear stress by migrating and remodelling the initially hyperbranched plexus, removing certain vessels whilst maintaining others. In this study, we argue that the key regulator of vessel preservation is cell decision behaviour at bifurcations. At flow-convergent bifurcations where migration paths diverge, cells must finely tune migration along both possible paths if the bifurcation is to persist. Experiments have demonstrated that disrupting the cells' ability to sense shear or the junction forces transmitted between cells impacts the preservation of bifurcations during the remodelling process. However, how these migratory cues integrate during cell decision making remains poorly understood. Therefore, we present the first agent-based model of endothelial cell flow-mediated migration suitable for interrogating the mechanisms behind bifurcation stability. The model simulates flow in a bifurcated vessel network composed of agents representing endothelial cells arranged into a lumen which migrate against flow. Upon approaching a bifurcation where more than one migration path exists, agents refer to a stochastic bifurcation rule which models the decision cells make as a combination of flow-based and collective-based migratory cues. With this rule, cells favour branches with relatively larger shear stress or cell number. We found that cells must integrate both cues nearly equally to maximise bifurcation stability. In simulations with stable bifurcations, we found competitive oscillations between flow and collective cues, and simulations that lost the bifurcation were unable to maintain these oscillations. The competition between these two cues is haemodynamic in origin, and demonstrates that a natural defence against bifurcation loss during remodelling exists: as vessel lumens narrow due to cell efflux, resistance to flow and shear stress increases, attracting new cells to enter and rescue the vessel from regression. Our work provides theoretical insight into the role of junction force transmission has in stabilising vasculature during remodelling and as an emergent mechanism to avoid functional shunting.

Neural information flow (NIF) provides a novel approach for system identification in neuroscience. It models the neural computations in multiple brain regions and can be trained end-to-end via stochastic gradient descent from noninvasive data. NIF models represent neural information processing via a network of coupled tensors, each encoding the representation of the sensory input contained in a brain region. The elements of these tensors can be interpreted as cortical columns whose activity encodes the presence of a specific feature in a spatiotemporal location. Each tensor is coupled to the measured data specific to a brain region via low-rank observation models that can be decomposed into the spatial, temporal and feature receptive fields of a localized neuronal population. Both these observation models and the convolutional weights defining the information processing within regions are learned end-to-end by predicting the neural signal during sensory stimulation. We trained a NIF model on the activity of early visual areas using a large-scale fMRI dataset recorded in a single participant. We show that we can recover plausible visual representations and population receptive fields that are consistent with empirical findings.

Individual cells show variability in their signaling dynamics that often correlates with phenotypic responses, indicating that cell-to-cell variability is not merely noise but can have functional consequences. Based on this observation, we reasoned that cell-to-cell variability under the same treatment condition could be explained in part by a single signaling motif that maps different upstream signals into a corresponding set of downstream responses. If this assumption holds, then repeated measurements of upstream and downstream signaling dynamics in a population of cells could provide information about the underlying signaling motif for a given pathway, even when no prior knowledge of that motif exists. To test these two hypotheses, we developed a computer algorithm called MISC (Motif Inference from Single Cells) that infers the underlying signaling motif from paired time-series measurements from individual cells. When applied to measurements of transcription factor and reporter gene expression in the yeast stress response, MISC predicted signaling motifs that were consistent with previous mechanistic models of transcription. The ability to detect the underlying mechanism became less certain when a cell's upstream signal was randomly paired with another cell's downstream response, demonstrating how averaging time-series measurements across a population obscures information about the underlying signaling mechanism. In some cases, motif predictions improved as more cells were added to the analysis. These results provide evidence that mechanistic information about cellular signaling networks can be systematically extracted from the dynamical patterns of single cells.

Spatio-temporal patterns of melanocytic proliferations observed in vivo are important for diagnosis but the mechanisms that produce them are poorly understood. Here we present an agent-based model for simulating the emergence of the main biologic patterns found in melanocytic proliferations. Our model portrays the extracellular matrix of the dermo-epidermal junction as a two-dimensional manifold and we simulate cellular migration in terms of geometric translations driven by adhesive, repulsive and random forces. Abstracted cellular functions and melanocyte-matrix interactions are modeled as stochastic events. For identification and validation we use visual renderings of simulated cell populations in a horizontal perspective that reproduce growth patterns observed in vivo by sequential dermatoscopy and corresponding vertical views that reproduce the arrangement of melanocytes observed in histopathologic sections. Our results show that a balanced interplay of proliferation and migration produces the typical reticular pattern of nevi, whereas the globular pattern involves additional cellular mechanisms. We further demonstrate that slight variations in the three basic cellular properties proliferation, migration, and adhesion are sufficient to produce a large variety of morphological appearances of nevi. We anticipate our model to be a starting point for the reproduction of more complex scenarios that will help to establish functional connections between abstracted microscopic behavior and macroscopic patterns in all types of melanocytic proliferations including melanoma.

Membrane proteins account for about one third of the cellular proteome, but it is still unclear how dynamic they are and how they establish functional contacts with cytoplasmic interaction partners. Here, we consider a membrane-integrated one-component receptor that also acts as a transcriptional activator, and analyze how it kinetically locates its specific binding site on the genome. We focus on the case of CadC, the pH receptor of the acid stress response Cad system in E. coli. CadC is a prime example of a one-component signaling protein that directly binds to its cognate target site on the chromosome to regulate transcription. We combined fluorescence microscopy experiments, mathematical analysis, and kinetic Monte Carlo simulations to probe this target search process. Using fluorescently labeled CadC, we measured the time from activation of the receptor until successful binding to the DNA in single cells, exploiting that stable receptor-DNA complexes are visible as fluorescent spots. Our experimental data indicate that CadC is highly mobile in the membrane and finds its target by a 2D diffusion and capture mechanism. DNA mobility is constrained due to the overall chromosome organization, but a labeled DNA locus in the vicinity of the target site appears sufficiently mobile to randomly come close to the membrane. Relocation of the DNA target site to a distant position on the chromosome had almost no effect on the mean search time, which was between four and five minutes in either case. However, a mutant strain with two binding sites displayed a mean search time that was reduced by about a factor of two. This behavior is consistent with simulations of a coarse-grained lattice model for the coupled dynamics of DNA within a cell volume and proteins on its surface. The model also rationalizes the experimentally determined distribution of search times. Overall our findings reveal that DNA target search does not present a much bigger kinetic challenge for membrane-integrated proteins than for cytoplasmic proteins. More generally, diffusion and capture mechanisms may be sufficient for bacterial membrane proteins to establish functional contacts with cytoplasmic targets.

Update of medRxiv. 2020 Jun 05;:

The SARS-CoV-2 pandemic of 2020 has mobilised scientists around the globe to research all aspects of the coronavirus virus and its infection. For fruitful and rapid investigation of viral pathomechanisms, a collaborative and interdisciplinary approach is required. Therefore, we have developed ViralLink: a systems biology workflow which reconstructs and analyses networks representing the effect of viruses on intracellular signalling. These networks trace the flow of signal from intracellular viral proteins through their human binding proteins and downstream signalling pathways, ending with transcription factors regulating genes differentially expressed upon viral exposure. In this way, the workflow provides a mechanistic insight from previously identified knowledge of virally infected cells. By default, the workflow is set up to analyse the intracellular effects of SARS-CoV-2, requiring only transcriptomics counts data as input from the user: thus, encouraging and enabling rapid multidisciplinary research. However, the wide-ranging applicability and modularity of the workflow facilitates customisation of viral context, a priori interactions and analysis methods. Through a case study of SARS-CoV-2 infected bronchial/tracheal epithelial cells, we evidence the functionality of the workflow and its ability to identify key pathways and proteins in the cellular response to infection. The application of ViralLink to different viral infections in a context specific manner using different available transcriptomics datasets will uncover key mechanisms in viral pathogenesis.

The aim of this study is to analyze patient movement patterns between hospital departments to derive the underlying intra-hospital movement network, and to assess if movement patterns differ between patients at high or low risk of colonization. For that purpose, we analyzed patient electronic medical record data from five hospitals to extract information on risk stratification and patient intra-hospital movements. Movement patterns were visualized as networks, and network centrality measures were calculated. Next, using an agent-based model where agents represent patients and intra-hospital patient movements were explicitly modeled, we simulated the spread of multidrug resistant enterobacteriacae (MDR-E) inside a hospital. Risk stratification of patients according to certain ICD-10 codes revealed that length of stay, patient age, and mean number of movements per admission were higher in the high-risk groups. Movement networks in all hospitals displayed a high variability among departments concerning their network centrality and connectedness with a few highly connected departments and many weakly connected peripheral departments. Simulating the spread of a pathogen in one hospital network showed positive correlation between department prevalence and network centrality measures. This study highlights the importance of intra-hospital patient movements and their possible impact on pathogen spread. Targeting interventions to departments of higher (weighted) degree may help to control the spread of MDR-E. Moreover, when the colonization status of patients coming from different departments is unknown, a ranking system based on department centralities may be used to design more effective interventions that mitigate pathogen spread.

Large-scale serological testing in the population is essential to determine the true extent of the current SARS-CoV-2 pandemic. Serological tests measure antibody responses against pathogens and use predefined cutoff levels that dichotomize the quantitative test measures into sero-positives and negatives and use this as a proxy for past infection. With the imperfect assays that are currently available to test for past SARS-CoV-2 infection, the fraction of seropositive individuals in serosurveys is a biased estimator of the cumulative incidence and is usually corrected to account for the sensitivity and specificity. Here we use an inference method-referred to as mixture-model approach-for the estimation of the cumulative incidence that does not require to define cutoffs by integrating the quantitative test measures directly into the statistical inference procedure. We confirm that the mixture model outperforms the methods based on cutoffs, leading to less bias and error in estimates of the cumulative incidence. We illustrate how the mixture model can be used to optimize the design of serosurveys with imperfect serological tests. We also provide guidance on the number of control and case sera that are required to quantify the test's ambiguity sufficiently to enable the reliable estimation of the cumulative incidence. Lastly, we show how this approach can be used to estimate the cumulative incidence of classes of infections with an unknown distribution of quantitative test measures. This is a very promising application of the mixture-model approach that could identify the elusive fraction of asymptomatic SARS-CoV-2 infections. An R-package implementing the inference methods used in this paper is provided. Our study advocates using serological tests without cutoffs, especially if they are used to determine parameters characterizing populations rather than individuals. This approach circumvents some of the shortcomings of cutoff-based methods at exactly the low cumulative incidence levels and test accuracies that we are currently facing in SARS-CoV-2 serosurveys.

The vast expansion of protein sequence databases provides an opportunity for new protein design approaches which seek to learn the sequence-function relationship directly from natural sequence variation. Deep generative models trained on protein sequence data have been shown to learn biologically meaningful representations helpful for a variety of downstream tasks, but their potential for direct use in the design of novel proteins remains largely unexplored. Here we show that variational autoencoders trained on a dataset of almost 70000 luciferase-like oxidoreductases can be used to generate novel, functional variants of the luxA bacterial luciferase. We propose separate VAE models to work with aligned sequence input (MSA VAE) and raw sequence input (AR-VAE), and offer evidence that while both are able to reproduce patterns of amino acid usage characteristic of the family, the MSA VAE is better able to capture long-distance dependencies reflecting the influence of 3D structure. To confirm the practical utility of the models, we used them to generate variants of luxA whose luminescence activity was validated experimentally. We further showed that conditional variants of both models could be used to increase the solubility of luxA without disrupting function. Altogether 6/12 of the variants generated using the unconditional AR-VAE and 9/11 generated using the unconditional MSA VAE retained measurable luminescence, together with all 23 of the less distant variants generated by conditional versions of the models; the most distant functional variant contained 35 differences relative to the nearest training set sequence. These results demonstrate the feasibility of using deep generative models to explore the space of possible protein sequences and generate useful variants, providing a method complementary to rational design and directed evolution approaches.

The physical basis of consciousness remains one of the most elusive concepts in current science. One influential conjecture is that consciousness is to do with some form of causality, measurable through information. The integrated information theory of consciousness (IIT) proposes that conscious experience, filled with rich and specific content, corresponds directly to a hierarchically organised, irreducible pattern of causal interactions; i.e. an integrated informational structure among elements of a system. Here, we tested this conjecture in a simple biological system (fruit flies), estimating the information structure of the system during wakefulness and general anesthesia. Consistent with this conjecture, we found that integrated interactions among populations of neurons during wakefulness collapsed to isolated clusters of interactions during anesthesia. We used classification analysis to quantify the accuracy of discrimination between wakeful and anesthetised states, and found that informational structures inferred conscious states with greater accuracy than a scalar summary of the structure, a measure which is generally championed as the main measure of IIT. In stark contrast to a view which assumes feedforward architecture for insect brains, especially fly visual systems, we found rich information structures, which cannot arise from purely feedforward systems, occurred across the fly brain. Further, these information structures collapsed uniformly across the brain during anesthesia. Our results speak to the potential utility of the novel concept of an "informational structure" as a measure for level of consciousness, above and beyond simple scalar values.

Low-cost, high-throughput sequencing has led to an enormous increase in the number of sequenced microbial genomes, with well over 100,000 genomes in public archives today. Automatic genome annotation tools are integral to understanding these organisms, yet older gene finding methods must be retrained on each new genome. We have developed a universal model of prokaryotic genes by fitting a temporal convolutional network to amino-acid sequences from a large, diverse set of microbial genomes. We incorporated the new model into a gene finding system, Balrog (Bacterial Annotation by Learned Representation Of Genes), which does not require genome-specific training and which matches or outperforms other state-of-the-art gene finding tools. Balrog is freely available under the MIT license at https://github.com/salzberg-lab/Balrog.

In this work, we show how the mechanical properties of the cellular microenvironment modulate the growth of tumour spheroids. Based on the composition of the extracellular matrix, its stiffness and architecture can significantly vary, subsequently influencing cell movement and tumour growth. However, it is still unclear exactly how both of these processes are regulated by the matrix composition. Here, we present a centre-based computational model that describes how collagen density, which modulates the steric hindrance properties of the matrix, governs individual cell migration and, consequently, leads to the formation of multicellular clusters of varying size. The model was calibrated using previously published experimental data, replicating a set of experiments in which cells were seeded in collagen matrices of different collagen densities, hence producing distinct mechanical properties. At an initial stage, we tracked individual cell trajectories and speeds. Subsequently, the formation of multicellular clusters was also analysed by quantifying their size. Overall, the results showed that our model could accurately replicate what was previously seen experimentally. Specifically, we showed that cells seeded in matrices with low collagen density tended to migrate more. Accordingly, cells strayed away from their original cluster and thus promoted the formation of small structures. In contrast, we also showed that high collagen densities hindered cell migration and produced multicellular clusters with increased volume. In conclusion, this model not only establishes a relation between matrix density and individual cell migration but also showcases how migration, or its inhibition, modulates tumour growth.

Intracranial electroencephalography is a standard tool in clinical evaluation of patients with focal epilepsy. Various early electrographic seizure patterns differing in frequency, amplitude, and waveform of the oscillations are observed. The pattern most common in the areas of seizure propagation is the so-called theta-alpha activity (TAA), whose defining features are oscillations in the θ - α range and gradually increasing amplitude. A deeper understanding of the mechanism underlying the generation of the TAA pattern is however lacking. In this work we evaluate the hypothesis that the TAA patterns are caused by seizures spreading across the cortex. To do so, we perform simulations of seizure dynamics on detailed patient-derived cortical surfaces using the spreading seizure model as well as reference models with one or two homogeneous sources. We then detect the occurrences of the TAA patterns both in the simulated stereo-electroencephalographic signals and in the signals of recorded epileptic seizures from a cohort of fifty patients, and we compare the features of the groups of detected TAA patterns to assess the plausibility of the different models. Our results show that spreading seizure hypothesis is qualitatively consistent with the evidence available in the seizure recordings, and it can explain the features of the detected TAA groups best among the examined models.

In this work we develop a novel algorithm for reconstructing the genomes of ancestral individuals, given genotype or sequence data from contemporary individuals and an extended pedigree of family relationships. A pedigree with complete genomes for every individual enables the study of allele frequency dynamics and haplotype diversity across generations, including deviations from neutrality such as transmission distortion. When studying heritable diseases, ancestral haplotypes can be used to augment genome-wide association studies and track disease inheritance patterns. The building blocks of our reconstruction algorithm are segments of Identity-By-Descent (IBD) shared between two or more genotyped individuals. The method alternates between identifying a source for each IBD segment and assembling IBD segments placed within each ancestral individual. Unlike previous approaches, our method is able to accommodate complex pedigree structures with hundreds of individuals genotyped at millions of SNPs. We apply our method to an Old Order Amish pedigree from Lancaster, Pennsylvania, whose founders came to North America from Europe during the early 18th century. The pedigree includes 1338 individuals from the past 12 generations, 394 with genotype data. The motivation for reconstruction is to understand the genetic basis of diseases segregating in the family through tracking haplotype transmission over time. Using our algorithm thread, we are able to reconstruct an average of 224 ancestral individuals per chromosome. For these ancestral individuals, on average we reconstruct 79% of their haplotypes. We also identify a region on chromosome 16 that is difficult to reconstruct-we find that this region harbors a short Amish-specific copy number variation and the gene HYDIN. thread was developed for endogamous populations, but can be applied to any extensive pedigree with the recent generations genotyped. We anticipate that this type of practical ancestral reconstruction will become more common and necessary to understand rare and complex heritable diseases in extended families.

With more microbiome studies being conducted by African-based research groups, there is an increasing demand for knowledge and skills in the design and analysis of microbiome studies and data. However, high-quality bioinformatics courses are often impeded by differences in computational environments, complicated software stacks, numerous dependencies, and versions of bioinformatics tools along with a lack of local computational infrastructure and expertise. To address this, H3ABioNet developed a 16S rRNA Microbiome Intermediate Bioinformatics Training course, extending its remote classroom model. The course was developed alongside experienced microbiome researchers, bioinformaticians, and systems administrators, who identified key topics to address. Development of containerised workflows has previously been undertaken by H3ABioNet, and Singularity containers were used here to enable the deployment of a standard replicable software stack across different hosting sites. The pilot ran successfully in 2019 across 23 sites registered in 11 African countries, with more than 200 participants formally enrolled and 106 volunteer staff for onsite support. The pulling, running, and testing of the containers, software, and analyses on various clusters were performed prior to the start of the course by hosting classrooms. The containers allowed the replication of analyses and results across all participating classrooms running a cluster and remained available posttraining ensuring analyses could be repeated on real data. Participants thus received the opportunity to analyse their own data, while local staff were trained and supported by experienced experts, increasing local capacity for ongoing research support. This provides a model for delivering topic-specific bioinformatics courses across Africa and other remote/low-resourced regions which overcomes barriers such as inadequate infrastructures, geographical distance, and access to expertise and educational materials.

Workflow management systems represent, manage, and execute multistep computational analyses and offer many benefits to bioinformaticians. They provide a common language for describing analysis workflows, contributing to reproducibility and to building libraries of reusable components. They can support both incremental build and re-entrancy-the ability to selectively re-execute parts of a workflow in the presence of additional inputs or changes in configuration and to resume execution from where a workflow previously stopped. Many workflow management systems enhance portability by supporting the use of containers, high-performance computing (HPC) systems, and clouds. Most importantly, workflow management systems allow bioinformaticians to delegate how their workflows are run to the workflow management system and its developers. This frees the bioinformaticians to focus on what these workflows should do, on their data analyses, and on their science. RiboViz is a package to extract biological insight from ribosome profiling data to help advance understanding of protein synthesis. At the heart of RiboViz is an analysis workflow, implemented in a Python script. To conform to best practices for scientific computing which recommend the use of build tools to automate workflows and to reuse code instead of rewriting it, the authors reimplemented this workflow within a workflow management system. To select a workflow management system, a rapid survey of available systems was undertaken, and candidates were shortlisted: Snakemake, cwltool, Toil, and Nextflow. Each candidate was evaluated by quickly prototyping a subset of the RiboViz workflow, and Nextflow was chosen. The selection process took 10 person-days, a small cost for the assurance that Nextflow satisfied the authors' requirements. The use of prototyping can offer a low-cost way of making a more informed selection of software to use within projects, rather than relying solely upon reviews and recommendations by others.

Increased availability of drug response and genomics data for many tumor cell lines has accelerated the development of pan-cancer prediction models of drug response. However, it is unclear how much between-tissue differences in drug response and molecular characteristics may contribute to pan-cancer predictions. Also unknown is whether the performance of pan-cancer models could vary by cancer type. Here, we built a series of pan-cancer models using two datasets containing 346 and 504 cell lines, each with MEK inhibitor (MEKi) response and mRNA expression, point mutation, and copy number variation data, and found that, while the tissue-level drug responses are accurately predicted (between-tissue ρ = 0.88-0.98), only 5 of 10 cancer types showed successful within-tissue prediction performance (within-tissue ρ = 0.11-0.64). Between-tissue differences make substantial contributions to the performance of pan-cancer MEKi response predictions, as exclusion of between-tissue signals leads to a decrease in Spearman's ρ from a range of 0.43-0.62 to 0.30-0.51. In practice, joint analysis of multiple cancer types usually has a larger sample size, hence greater power, than for one cancer type; and we observe that higher accuracy of pan-cancer prediction of MEKi response is almost entirely due to the sample size advantage. Success of pan-cancer prediction reveals how drug response in different cancers may invoke shared regulatory mechanisms despite tissue-specific routes of oncogenesis, yet predictions in different cancer types require flexible incorporation of between-cancer and within-cancer signals. As most datasets in genome sciences contain multiple levels of heterogeneity, careful parsing of group characteristics and within-group, individual variation is essential when making robust inference.

We live in an increasingly data-driven world, where high-throughput sequencing and mass spectrometry platforms are transforming biology into an information science. This has shifted major challenges in biological research from data generation and processing to interpretation and knowledge translation. However, postsecondary training in bioinformatics, or more generally data science for life scientists, lags behind current demand. In particular, development of accessible, undergraduate data science curricula has the potential to improve research and learning outcomes as well as better prepare students in the life sciences to thrive in public and private sector careers. Here, we describe the Experiential Data science for Undergraduate Cross-Disciplinary Education (EDUCE) initiative, which aims to progressively build data science competency across several years of integrated practice. Through EDUCE, students complete data science modules integrated into required and elective courses augmented with coordinated cocurricular activities. The EDUCE initiative draws on a community of practice consisting of teaching assistants (TAs), postdocs, instructors, and research faculty from multiple disciplines to overcome several reported barriers to data science for life scientists, including instructor capacity, student prior knowledge, and relevance to discipline-specific problems. Preliminary survey results indicate that even a single module improves student self-reported interest and/or experience in bioinformatics and computer science. Thus, EDUCE provides a flexible and extensible active learning framework for integration of data science curriculum into undergraduate courses and programs across the life sciences.

Erratum in PLoS Comput Biol. 2021 Dec 10;17(12):e1009686.

The collective behaviour of animal and human groups emerges from the individual decisions and actions of their constituent members. Recent research has revealed many ways in which the behaviour of groups can be influenced by differences amongst their constituent individuals. The existence of individual differences that have implications for collective behaviour raises important questions. How are these differences generated and maintained? Are individual differences driven by exogenous factors, or are they a response to the social dilemmas these groups face? Here I consider the classic case of patch selection by foraging agents under conditions of social competition. I introduce a multilevel model wherein the perceptual sensitivities of agents evolve in response to their foraging success or failure over repeated patch selections. This model reveals a bifurcation in the population, creating a class of agents with no perceptual sensitivity. These agents exploit the social environment to avoid the costs of accurate perception, relying on other agents to make fitness rewards insensitive to the choice of foraging patch. This provides a individual-based evolutionary basis for models incorporating perceptual limits that have been proposed to explain observed deviations from the Ideal Free Distribution (IFD) in empirical studies, while showing that the common assumption in such models that agents share identical sensory limits is likely false. Further analysis of the model shows how agents develop perceptual strategic niches in response to environmental variability. The emergence of agents insensitive to reward differences also has implications for societal resource allocation problems, including the use of financial and prediction markets as mechanisms for aggregating collective wisdom.

Crystallography and NMR system (CNS) is currently a widely used method for fragment-free ab initio protein folding from inter-residue distance or contact maps. Despite its widespread use in protein structure prediction, CNS is a decade-old macromolecular structure determination system that was originally developed for solving macromolecular geometry from experimental restraints as opposed to predictive modeling driven by interaction map data. As such, the adaptation of the CNS experimental structure determination protocol for ab initio protein folding is intrinsically anomalous that may undermine the folding accuracy of computational protein structure prediction. In this paper, we propose a new CNS-free hierarchical structure modeling method called DConStruct for folding both soluble and membrane proteins driven by distance and contact information. Rigorous experimental validation shows that DConStruct attains much better reconstruction accuracy than CNS when tested with the same input contact map at varying contact thresholds. The hierarchical modeling with iterative self-correction employed in DConStruct scales at a much higher degree of folding accuracy than CNS with the increase in contact thresholds, ultimately approaching near-optimal reconstruction accuracy at higher-thresholded contact maps. The folding accuracy of DConStruct can be further improved by exploiting distance-based hybrid interaction maps at tri-level thresholding, as demonstrated by the better performance of our method in folding free modeling targets from the 12th and 13th rounds of the Critical Assessment of techniques for protein Structure Prediction (CASP) experiments compared to popular CNS- and fragment-based approaches and energy-minimization protocols, some of which even using much finer-grained distance maps than ours. Additional large-scale benchmarking shows that DConStruct can significantly improve the folding accuracy of membrane proteins compared to a CNS-based approach. These results collectively demonstrate the feasibility of greatly improving the accuracy of ab initio protein folding by optimally exploiting the information encoded in inter-residue interaction maps beyond what is possible by CNS.

Neurons exhibit diverse intrinsic dynamics, which govern how they integrate synaptic inputs to produce spikes. Intrinsic dynamics are often plastic during development and learning, but the effects of these changes on stimulus encoding properties are not well known. To examine this relationship, we simulated auditory responses to zebra finch song using a linear-dynamical cascade model, which combines a linear spectrotemporal receptive field with a dynamical, conductance-based neuron model, then used generalized linear models to estimate encoding properties from the resulting spike trains. We focused on the effects of a low-threshold potassium current (KLT) that is present in a subset of cells in the zebra finch caudal mesopallium and is affected by early auditory experience. We found that KLT affects both spike adaptation and the temporal filtering properties of the receptive field. The direction of the effects depended on the temporal modulation tuning of the linear (input) stage of the cascade model, indicating a strongly nonlinear relationship. These results suggest that small changes in intrinsic dynamics in tandem with differences in synaptic connectivity can have dramatic effects on the tuning of auditory neurons.

Significant scientific and translational questions remain in auditory neuroscience surrounding the neural correlates of perception. Relating perceptual and neural data collected from humans can be useful; however, human-based neural data are typically limited to evoked far-field responses, which lack anatomical and physiological specificity. Laboratory-controlled preclinical animal models offer the advantage of comparing single-unit and evoked responses from the same animals. This ability provides opportunities to develop invaluable insight into proper interpretations of evoked responses, which benefits both basic-science studies of neural mechanisms and translational applications, e.g., diagnostic development. However, these comparisons have been limited by a disconnect between the types of spectrotemporal analyses used with single-unit spike trains and evoked responses, which results because these response types are fundamentally different (point-process versus continuous-valued signals) even though the responses themselves are related. Here, we describe a unifying framework to study temporal coding of complex sounds that allows spike-train and evoked-response data to be analyzed and compared using the same advanced signal-processing techniques. The framework uses a set of peristimulus-time histograms computed from single-unit spike trains in response to polarity-alternating stimuli to allow advanced spectral analyses of both slow (envelope) and rapid (temporal fine structure) response components. Demonstrated benefits include: (1) novel spectrally specific temporal-coding measures that are less confounded by distortions due to hair-cell transduction, synaptic rectification, and neural stochasticity compared to previous metrics, e.g., the correlogram peak-height, (2) spectrally specific analyses of spike-train modulation coding (magnitude and phase), which can be directly compared to modern perceptually based models of speech intelligibility (e.g., that depend on modulation filter banks), and (3) superior spectral resolution in analyzing the neural representation of nonstationary sounds, such as speech and music. This unifying framework significantly expands the potential of preclinical animal models to advance our understanding of the physiological correlates of perceptual deficits in real-world listening following sensorineural hearing loss.

Proteases are an important class of enzymes, whose activity is central to many physiologic and pathologic processes. Detailed knowledge of protease specificity is key to understanding their function. Although many methods have been developed to profile specificities of proteases, few have the diversity and quantitative grasp necessary to fully define specificity of a protease, both in terms of substrate numbers and their catalytic efficiencies. We have developed a concept of "selectome"; the set of substrate amino acid sequences that uniquely represent the specificity of a protease. We applied it to two closely related members of the Matrixin family-MMP-2 and MMP-9 by using substrate phage display coupled with Next Generation Sequencing and information theory-based data analysis. We have also derived a quantitative measure of substrate specificity, which accounts for both the number of substrates and their relative catalytic efficiencies. Using these advances greatly facilitates elucidation of substrate selectivity between closely related members of a protease family. The study also provides insight into the degree to which the catalytic cleft defines substrate recognition, thus providing basis for overcoming two of the major challenges in the field of proteolysis: 1) development of highly selective activity probes for studying proteases with overlapping specificities, and 2) distinguishing targeted proteolysis from bystander proteolytic events.

Phenotypic profiling of large three-dimensional microscopy data sets has not been widely adopted due to the challenges posed by cell segmentation and feature selection. The computational demands of automated processing further limit analysis of hard-to-segment images such as of neurons and organoids. Here we describe a comprehensive shallow-learning framework for automated quantitative phenotyping of three-dimensional (3D) image data using unsupervised data-driven voxel-based feature learning, which enables computationally facile classification, clustering and advanced data visualization. We demonstrate the analysis potential on complex 3D images by investigating the phenotypic alterations of: neurons in response to apoptosis-inducing treatments and morphogenesis for oncogene-expressing human mammary gland acinar organoids. Our novel implementation of image analysis algorithms called Phindr3D allowed rapid implementation of data-driven voxel-based feature learning into 3D high content analysis (HCA) operations and constitutes a major practical advance as the computed assignments represent the biology while preserving the heterogeneity of the underlying data. Phindr3D is provided as Matlab code and as a stand-alone program (https://github.com/DWALab/Phindr3D).

Biomineralization is the process by which organisms use minerals to harden their tissues and provide them with physical support. Biomineralizing cells concentrate the mineral in vesicles that they secret into a dedicated compartment where crystallization occurs. The dynamics of vesicle motion and the molecular mechanisms that control it, are not well understood. Sea urchin larval skeletogenesis provides an excellent platform for investigating the kinetics of mineral-bearing vesicles. Here we used lattice light-sheet microscopy to study the three-dimensional (3D) dynamics of calcium-bearing vesicles in the cells of normal sea urchin embryos and of embryos where skeletogenesis is blocked through the inhibition of Vascular Endothelial Growth Factor Receptor (VEGFR). We developed computational tools for displaying 3D-volumetric movies and for automatically quantifying vesicle dynamics. Our findings imply that calcium vesicles perform an active diffusion motion in both, calcifying (skeletogenic) and non-calcifying (ectodermal) cells of the embryo. The diffusion coefficient and vesicle speed are larger in the mesenchymal skeletogenic cells compared to the epithelial ectodermal cells. These differences are possibly due to the distinct mechanical properties of the two tissues, demonstrated by the enhanced f-actin accumulation and myosinII activity in the ectodermal cells compared to the skeletogenic cells. Vesicle motion is not directed toward the biomineralization compartment, but the vesicles slow down when they approach it, and probably bind for mineral deposition. VEGFR inhibition leads to an increase of vesicle volume but hardly changes vesicle kinetics and doesn't affect f-actin accumulation and myosinII activity. Thus, calcium vesicles perform an active diffusion motion in the cells of the sea urchin embryo, with diffusion length and speed that inversely correlate with the strength of the actomyosin network. Overall, our studies provide an unprecedented view of calcium vesicle 3D-dynamics and point toward cytoskeleton remodeling as an important effector of the motion of mineral-bearing vesicles.

Approximately 30% of patients who have Clostridioides difficile infection (CDI) will suffer at least one incident of reinfection. While the underlying causes of CDI recurrence are poorly understood, interactions between C. difficile and commensal gut bacteria are thought to play an important role. In this study, an in silico pipeline was used to process 16S rRNA gene amplicon sequence data of 225 stool samples from 93 CDI patients into sample-specific models of bacterial community metabolism. Clustered metabolite production rates generated from post-diagnosis samples generated a high Enterobacteriaceae abundance cluster containing disproportionately large numbers of recurrent samples and patients. This cluster was predicted to have significantly reduced capabilities for secondary bile acid synthesis but elevated capabilities for aromatic amino acid catabolism. When applied to 16S sequence data of 40 samples from fecal microbiota transplantation (FMT) patients suffering from recurrent CDI and their stool donors, the community modeling method generated a high Enterobacteriaceae abundance cluster with a disproportionate large number of pre-FMT samples. This cluster also was predicted to exhibit reduced secondary bile acid synthesis and elevated aromatic amino acid catabolism. Collectively, these in silico predictions suggest that Enterobacteriaceae may create a gut environment favorable for C. difficile spore germination and/or toxin synthesis.

Effective regulation of the sonic hedgehog (Shh) signalling pathway is essential for normal development in a wide variety of species. Correct Shh signalling requires the formation of Shh aggregates on the surface of producing cells. Shh aggregates subsequently diffuse away and are recognised in receiving cells located elsewhere in the developing embryo. Various mechanisms have been postulated regarding how these aggregates form and what their precise role is in the overall signalling process. To understand the role of these mechanisms in the overall signalling process, we formulate and analyse a mathematical model of Shh aggregation using nonlinear ordinary differential equations. We consider Shh aggregate formation to comprise of multimerisation, association with heparan sulfate proteoglycans (HSPG) and binding with lipoproteins. We show that the size distribution of the Shh aggregates formed on the producing cell surface resembles an exponential distribution, a result in agreement with experimental data. A detailed sensitivity analysis of our model reveals that this exponential distribution is robust to parameter changes, and subsequently, also to variations in the processes by which Shh is recruited by HSPGs and lipoproteins. The work demonstrates the time taken for different sized Shh aggregates to form and the important role this likely plays in Shh diffusion.

Tracking statistical regularities of the environment is important for shaping human behavior and perception. Evidence suggests that the brain learns environmental dependencies using Bayesian principles. However, much remains unknown about the employed algorithms, for somesthesis in particular. Here, we describe the cortical dynamics of the somatosensory learning system to investigate both the form of the generative model as well as its neural surprise signatures. Specifically, we recorded EEG data from 40 participants subjected to a somatosensory roving-stimulus paradigm and performed single-trial modeling across peri-stimulus time in both sensor and source space. Our Bayesian model selection procedure indicates that evoked potentials are best described by a non-hierarchical learning model that tracks transitions between observations using leaky integration. From around 70ms post-stimulus onset, secondary somatosensory cortices are found to represent confidence-corrected surprise as a measure of model inadequacy. Indications of Bayesian surprise encoding, reflecting model updating, are found in primary somatosensory cortex from around 140ms. This dissociation is compatible with the idea that early surprise signals may control subsequent model update rates. In sum, our findings support the hypothesis that early somatosensory processing reflects Bayesian perceptual learning and contribute to an understanding of its underlying mechanisms.

A population's spatial structure affects the rate of genetic change and the outcome of natural selection. These effects can be modeled mathematically using the Birth-death process on graphs. Individuals occupy the vertices of a weighted graph, and reproduce into neighboring vertices based on fitness. A key quantity is the probability that a mutant type will sweep to fixation, as a function of the mutant's fitness. Graphs that increase the fixation probability of beneficial mutations, and decrease that of deleterious mutations, are said to amplify selection. However, fixation probabilities are difficult to compute for an arbitrary graph. Here we derive an expression for the fixation probability, of a weakly-selected mutation, in terms of the time for two lineages to coalesce. This expression enables weak-selection fixation probabilities to be computed, for an arbitrary weighted graph, in polynomial time. Applying this method, we explore the range of possible effects of graph structure on natural selection, genetic drift, and the balance between the two. Using exhaustive analysis of small graphs and a genetic search algorithm, we identify families of graphs with striking effects on fixation probability, and we analyze these families mathematically. Our work reveals the nuanced effects of graph structure on natural selection and neutral drift. In particular, we show how these notions depend critically on the process by which mutations arise.

The availability of bacterial transcriptomes has dramatically increased in recent years. This data deluge could result in detailed inference of underlying regulatory networks, but the diversity of experimental platforms and protocols introduces critical biases that could hinder scalable analysis of existing data. Here, we show that the underlying structure of the E. coli transcriptome, as determined by Independent Component Analysis (ICA), is conserved across multiple independent datasets, including both RNA-seq and microarray datasets. We subsequently combined five transcriptomics datasets into a large compendium containing over 800 expression profiles and discovered that its underlying ICA-based structure was still comparable to that of the individual datasets. With this understanding, we expanded our analysis to over 3,000 E. coli expression profiles and predicted three high-impact regulons that respond to oxidative stress, anaerobiosis, and antibiotic treatment. ICA thus enables deep analysis of disparate data to uncover new insights that were not visible in the individual datasets.

Finding non-standard or new metabolic pathways has important applications in metabolic engineering, synthetic biology and the analysis and reconstruction of metabolic networks. Branched metabolic pathways dominate in metabolic networks and depict a more comprehensive picture of metabolism compared to linear pathways. Although progress has been developed to find branched metabolic pathways, few efforts have been made in identifying branched metabolic pathways via atom group tracking. In this paper, we present a pathfinding method called BPFinder for finding branched metabolic pathways by atom group tracking, which aims to guide the synthetic design of metabolic pathways. BPFinder enumerates linear metabolic pathways by tracking the movements of atom groups in metabolic network and merges the linear atom group conserving pathways into branched pathways. Two merging rules based on the structure of conserved atom groups are proposed to accurately merge the branched compounds of linear pathways to identify branched pathways. Furthermore, the integrated information of compound similarity, thermodynamic feasibility and conserved atom groups is also used to rank the pathfinding results for feasible branched pathways. Experimental results show that BPFinder is more capable of recovering known branched metabolic pathways as compared to other existing methods, and is able to return biologically relevant branched pathways and discover alternative branched pathways of biochemical interest. The online server of BPFinder is available at http://114.215.129.245:8080/atomic/. The program, source code and data can be downloaded from https://github.com/hyr0771/BPFinder.

Relaxed clock models enable estimation of molecular substitution rates across lineages and are widely used in phylogenetics for dating evolutionary divergence times. Under the (uncorrelated) relaxed clock model, tree branches are associated with molecular substitution rates which are independently and identically distributed. In this article we delved into the internal complexities of the relaxed clock model in order to develop efficient MCMC operators for Bayesian phylogenetic inference. We compared three substitution rate parameterisations, introduced an adaptive operator which learns the weights of other operators during MCMC, and we explored how relaxed clock model estimation can benefit from two cutting-edge proposal kernels: the AVMVN and Bactrian kernels. This work has produced an operator scheme that is up to 65 times more efficient at exploring continuous relaxed clock parameters compared with previous setups, depending on the dataset. Finally, we explored variants of the standard narrow exchange operator which are specifically designed for the relaxed clock model. In the most extreme case, this new operator traversed tree space 40% more efficiently than narrow exchange. The methodologies introduced are adaptive and highly effective on short as well as long alignments. The results are available via the open source optimised relaxed clock (ORC) package for BEAST 2 under a GNU licence (https://github.com/jordandouglas/ORC).

In tumor metastasis, the margination and adhesion of tumor cells are two critical and closely related steps, which may determine the destination where the tumor cells extravasate to. We performed a direct three-dimensional simulation on the behaviors of the tumor cells in a real microvascular network, by a hybrid method of the smoothed dissipative particle dynamics and immersed boundary method (SDPD-IBM). The tumor cells are found to adhere at the microvascular bifurcations more frequently, and there is a positive correlation between the adhesion of the tumor cells and the wall-directed force from the surrounding red blood cells (RBCs). The larger the wall-directed force is, the closer the tumor cells are marginated towards the wall, and the higher the probability of adhesion behavior happen is. A relatively low or high hematocrit can help to prevent the adhesion of tumor cells, and similarly, increasing the shear rate of blood flow can serve the same purpose. These results suggest that the tumor cells may be more likely to extravasate at the microvascular bifurcations if the blood flow is slow and the hematocrit is moderate.

Rare variants are thought to play an important role in the etiology of complex diseases and may explain a significant fraction of the missing heritability in genetic disease studies. Next-generation sequencing facilitates the association of rare variants in coding or regulatory regions with complex diseases in large cohorts at genome-wide scale. However, rare variant association studies (RVAS) still lack power when cohorts are small to medium-sized and if genetic variation explains a small fraction of phenotypic variance. Here we present a novel Bayesian rare variant Association Test using Integrated Nested Laplace Approximation (BATI). Unlike existing RVAS tests, BATI allows integration of individual or variant-specific features as covariates, while efficiently performing inference based on full model estimation. We demonstrate that BATI outperforms established RVAS methods on realistic, semi-synthetic whole-exome sequencing cohorts, especially when using meaningful biological context, such as functional annotation. We show that BATI achieves power above 70% in scenarios in which competing tests fail to identify risk genes, e.g. when risk variants in sum explain less than 0.5% of phenotypic variance. We have integrated BATI, together with five existing RVAS tests in the 'Rare Variant Genome Wide Association Study' (rvGWAS) framework for data analyzed by whole-exome or whole genome sequencing. rvGWAS supports rare variant association for genes or any other biological unit such as promoters, while allowing the analysis of essential functionalities like quality control or filtering. Applying rvGWAS to a Chronic Lymphocytic Leukemia study we identified eight candidate predisposition genes, including EHMT2 and COPS7A.

Erratum in PLoS Comput Biol. 2021 Apr 21;17(4):e1008945.

Positive interactions, including intraspecies cooperation and interspecies mutualisms, play crucial roles in shaping the structure and function of many ecosystems, ranging from plant communities to the human microbiome. While the evolutionary forces that form and maintain positive interactions have been investigated extensively, the influence of positive interactions on the ability of species to adapt to new environments is still poorly understood. Here, we use numerical simulations and theoretical analyses to study how positive interactions impact the likelihood that populations survive after an environment deteriorates, such that survival in the new environment requires quick adaptation via the rise of new mutants-a scenario known as evolutionary rescue. We find that the probability of evolutionary rescue in populations engaged in positive interactions is reduced significantly. In cooperating populations, this reduction is largely due to the fact that survival may require at least a minimal number of individuals, meaning that adapted mutants must arise and spread before the population declines below this threshold. In mutualistic populations, the rescue probability is decreased further due to two additional effects-the need for both mutualistic partners to adapt to the new environment, and competition between the two species. Finally, we show that the presence of cheaters reduces the likelihood of evolutionary rescue even further, making it extremely unlikely. These results indicate that while positive interactions may be beneficial in stable environments, they can hinder adaptation to changing environments and thereby elevate the risk of population collapse. Furthermore, these results may hint at the selective pressures that drove co-dependent unicellular species to form more adaptable organisms able to differentiate into multiple phenotypes, including multicellular life.

Remote acquisition of information on ecosystem dynamics is essential for conservation management, especially for the deep ocean. Soundscape offers unique opportunities to study the behavior of soniferous marine animals and their interactions with various noise-generating activities at a fine temporal resolution. However, the retrieval of soundscape information remains challenging owing to limitations in audio analysis techniques that are effective in the face of highly variable interfering sources. This study investigated the application of a seafloor acoustic observatory as a long-term platform for observing marine ecosystem dynamics through audio source separation. A source separation model based on the assumption of source-specific periodicity was used to factorize time-frequency representations of long-duration underwater recordings. With minimal supervision, the model learned to discriminate source-specific spectral features and prove to be effective in the separation of sounds made by cetaceans, soniferous fish, and abiotic sources from the deep-water soundscapes off northeastern Taiwan. Results revealed phenological differences among the sound sources and identified diurnal and seasonal interactions between cetaceans and soniferous fish. The application of clustering to source separation results generated a database featuring the diversity of soundscapes and revealed a compositional shift in clusters of cetacean vocalizations and fish choruses during diurnal and seasonal cycles. The source separation model enables the transformation of single-channel audio into multiple channels encoding the dynamics of biophony, geophony, and anthropophony, which are essential for characterizing the community of soniferous animals, quality of acoustic habitat, and their interactions. Our results demonstrated the application of source separation could facilitate acoustic diversity assessment, which is a crucial task in soundscape-based ecosystem monitoring. Future implementation of soundscape information retrieval in long-term marine observation networks will lead to the use of soundscapes as a new tool for conservation management in an increasingly noisy ocean.

N6-methyladenine (6mA) is an important DNA modification form associated with a wide range of biological processes. Identifying accurately 6mA sites on a genomic scale is crucial for under-standing of 6mA's biological functions. However, the existing experimental techniques for detecting 6mA sites are cost-ineffective, which implies the great need of developing new computational methods for this problem. In this paper, we developed, without requiring any prior knowledge of 6mA and manually crafted sequence features, a deep learning framework named Deep6mA to identify DNA 6mA sites, and its performance is superior to other DNA 6mA prediction tools. Specifically, the 5-fold cross-validation on a benchmark dataset of rice gives the sensitivity and specificity of Deep6mA as 92.96% and 95.06%, respectively, and the overall prediction accuracy is 94%. Importantly, we find that the sequences with 6mA sites share similar patterns across different species. The model trained with rice data predicts well the 6mA sites of other three species: Arabidopsis thaliana, Fragaria vesca and Rosa chinensis with a prediction accuracy over 90%. In addition, we find that (1) 6mA tends to occur at GAGG motifs, which means the sequence near the 6mA site may be conservative; (2) 6mA is enriched in the TATA box of the promoter, which may be the main source of its regulating downstream gene expression.

We propose a developmental model inspired by the cortico-basal system (CX-BG) for vocal learning in babies and for solving the correspondence mismatch problem they face when they hear unfamiliar voices, with different tones and pitches. This model is based on the neural architecture INFERNO standing for Iterative Free-Energy Optimization of Recurrent Neural Networks. Free-energy minimization is used for rapidly exploring, selecting and learning the optimal choices of actions to perform (eg sound production) in order to reproduce and control as accurately as possible the spike trains representing desired perceptions (eg sound categories). We detail in this paper the CX-BG system responsible for linking causally the sound and motor primitives at the order of a few milliseconds. Two experiments performed with a small and a large audio database show the capabilities of exploration, generalization and robustness to noise of our neural architecture in retrieving audio primitives during vocal learning and during acoustic matching with unheared voices (different genders and tones).

Gene function annotation is important for a variety of downstream analyses of genetic data. But experimental characterization of function remains costly and slow, making computational prediction an important endeavor. Phylogenetic approaches to prediction have been developed, but implementation of a practical Bayesian framework for parameter estimation remains an outstanding challenge. We have developed a computationally efficient model of evolution of gene annotations using phylogenies based on a Bayesian framework using Markov Chain Monte Carlo for parameter estimation. Unlike previous approaches, our method is able to estimate parameters over many different phylogenetic trees and functions. The resulting parameters agree with biological intuition, such as the increased probability of function change following gene duplication. The method performs well on leave-one-out cross-validation, and we further validated some of the predictions in the experimental scientific literature.

Segregation and integration are two fundamental principles of brain structural and functional organization. Neuroimaging studies have shown that the brain transits between different functionally segregated and integrated states, and neuromodulatory systems have been proposed as key to facilitate these transitions. Although whole-brain computational models have reproduced this neuromodulatory effect, the role of local inhibitory circuits and their cholinergic modulation has not been studied. In this article, we consider a Jansen & Rit whole-brain model in a network interconnected using a human connectome, and study the influence of the cholinergic and noradrenergic neuromodulatory systems on the segregation/integration balance. In our model, we introduce a local inhibitory feedback as a plausible biophysical mechanism that enables the integration of whole-brain activity, and that interacts with the other neuromodulatory influences to facilitate the transition between different functional segregation/integration regimes in the brain.

Lung cancer is one of the leading causes of cancer-related deaths worldwide and is characterized by hijacking immune system for active growth and aggressive metastasis. Neutrophils, which in their original form should establish immune activities to the tumor as a first line of defense, are undermined by tumor cells to promote tumor invasion in several ways. In this study, we investigate the mutual interactions between the tumor cells and the neutrophils that facilitate tumor invasion by developing a mathematical model that involves taxis-reaction-diffusion equations for the critical components in the interaction. These include the densities of tumor and neutrophils, and the concentrations of signaling molecules and structure such as neutrophil extracellular traps (NETs). We apply the mathematical model to a Boyden invasion assay used in the experiments to demonstrate that the tumor-associated neutrophils can enhance tumor cell invasion by secreting the neutrophil elastase. We show that the model can both reproduce the major experimental observation on NET-mediated cancer invasion and make several important predictions to guide future experiments with the goal of the development of new anti-tumor strategies. Moreover, using this model, we investigate the fundamental mechanism of NET-mediated invasion of cancer cells and the impact of internal and external heterogeneity on the migration patterning of tumour cells and their response to different treatment schedules.

Surgical interventions in epileptic patients aimed at the removal of the epileptogenic zone have success rates at only 60-70%. This failure can be partly attributed to the insufficient spatial sampling by the implanted intracranial electrodes during the clinical evaluation, leading to an incomplete picture of spatio-temporal seizure organization in the regions that are not directly observed. Utilizing the partial observations of the seizure spreading through the brain network, complemented by the assumption that the epileptic seizures spread along the structural connections, we infer if and when are the unobserved regions recruited in the seizure. To this end we introduce a data-driven model of seizure recruitment and propagation across a weighted network, which we invert using the Bayesian inference framework. Using a leave-one-out cross-validation scheme on a cohort of 45 patients we demonstrate that the method can improve the predictions of the states of the unobserved regions compared to an empirical estimate that does not use the structural information, yet it is on the same level as the estimate that takes the structure into account. Furthermore, a comparison with the performed surgical resection and the surgery outcome indicates a link between the inferred excitable regions and the actual epileptogenic zone. The results emphasize the importance of the structural connectome in the large-scale spatio-temporal organization of epileptic seizures and introduce a novel way to integrate the patient-specific connectome and intracranial seizure recordings in a whole-brain computational model of seizure spread.

Computational models of animal biosonar seek to identify critical aspects of echo processing responsible for the superior, real-time performance of echolocating bats and dolphins in target tracking and clutter rejection. The Spectrogram Correlation and Transformation (SCAT) model replicates aspects of biosonar imaging in both species by processing wideband biosonar sounds and echoes with auditory mechanisms identified from experiments with bats. The model acquires broadband biosonar broadcasts and echoes, represents them as time-frequency spectrograms using parallel bandpass filters, translates the filtered signals into ten parallel amplitude threshold levels, and then operates on the resulting time-of-occurrence values at each frequency to estimate overall echo range delay. It uses the structure of the echo spectrum by depicting it as a series of local frequency nulls arranged regularly along the frequency axis of the spectrograms after dechirping them relative to the broadcast. Computations take place entirely on the timing of threshold-crossing events for each echo relative to threshold-events for the broadcast. Threshold-crossing times take into account amplitude-latency trading, a physiological feature absent from conventional digital signal processing. Amplitude-latency trading transposes the profile of amplitudes across frequencies into a profile of time-registrations across frequencies. Target shape is extracted from the spacing of the object's individual acoustic reflecting points, or glints, using the mutual interference pattern of peaks and nulls in the echo spectrum. These are merged with the overall range-delay estimate to produce a delay-based reconstruction of the object's distance as well as its glints. Clutter echoes indiscriminately activate multiple parts in the null-detecting system, which then produces the equivalent glint-delay spacings in images, thus blurring the overall echo-delay estimates by adding spurious glint delays to the image. Blurring acts as an anticorrelation process that rejects clutter intrusion into perceptions.

Evolutionary branching occurs when a population with a unimodal phenotype distribution diversifies into a multimodally distributed population consisting of two or more strains. Branching results from frequency-dependent selection, which is caused by interactions between individuals. For example, a population performing a social task may diversify into a cooperator strain and a defector strain. Branching can also occur in multi-dimensional phenotype spaces, such as when two tasks are performed simultaneously. In such cases, the strains may diverge in different directions: possible outcomes include division of labor (with each population performing one of the tasks) or the diversification into a strain that performs both tasks and another that performs neither. Here we show that the shape of the population's phenotypic distribution plays a role in determining the direction of branching. Furthermore, we show that the shape of the distribution is, in turn, contingent on the direction of approach to the evolutionary branching point. This results in a distribution-selection feedback that is not captured in analytical models of evolutionary branching, which assume monomorphic populations. Finally, we show that this feedback can influence long-term evolutionary dynamics and promote the evolution of division of labor.

Spectral similarity is used as a proxy for structural similarity in many tandem mass spectrometry (MS/MS) based metabolomics analyses such as library matching and molecular networking. Although weaknesses in the relationship between spectral similarity scores and the true structural similarities have been described, little development of alternative scores has been undertaken. Here, we introduce Spec2Vec, a novel spectral similarity score inspired by a natural language processing algorithm-Word2Vec. Spec2Vec learns fragmental relationships within a large set of spectral data to derive abstract spectral embeddings that can be used to assess spectral similarities. Using data derived from GNPS MS/MS libraries including spectra for nearly 13,000 unique molecules, we show how Spec2Vec scores correlate better with structural similarity than cosine-based scores. We demonstrate the advantages of Spec2Vec in library matching and molecular networking. Spec2Vec is computationally more scalable allowing structural analogue searches in large databases within seconds.

Developmental mechanisms that canalize or compensate perturbations of organismal development (targeted or compensatory growth) are widely considered a prerequisite of individual health and the evolution of complex life, but little is known about the nature of these mechanisms. It is even unclear if and how a "target trajectory" of individual development is encoded in the organism's genetic-developmental system or, instead, emerges as an epiphenomenon. Here we develop a statistical model of developmental canalization based on an extended autoregressive model. We show that under certain assumptions the strength of canalization and the amount of canalized variance in a population can be estimated, or at least approximated, from longitudinal phenotypic measurements, even if the target trajectories are unobserved. We extend this model to multivariate measures and discuss reifications of the ensuing parameter matrix. We apply these approaches to longitudinal geometric morphometric data on human postnatal craniofacial size and shape as well as to the size of the frontal sinuses. Craniofacial size showed strong developmental canalization during the first 5 years of life, leading to a 50% reduction of cross-sectional size variance, followed by a continual increase in variance during puberty. Frontal sinus size, by contrast, did not show any signs of canalization. Total variance of craniofacial shape decreased slightly until about 5 years of age and increased thereafter. However, different features of craniofacial shape showed very different developmental dynamics. Whereas the relative dimensions of the nasopharynx showed strong canalization and a reduction of variance throughout postnatal development, facial orientation continually increased in variance. Some of the signals of canalization may owe to independent variation in developmental timing of cranial components, but our results indicate evolved, partly mechanically induced mechanisms of canalization that ensure properly sized upper airways and facial dimensions.

Contemporary accounts of the initiation of cardiac arrhythmias typically rely on after-depolarizations as the trigger for reentrant activity. The after-depolarizations are usually triggered by calcium entry or spontaneous release within the cells of the myocardium or the conduction system. Here we propose an alternative mechanism whereby arrhythmias are triggered autonomously by cardiac cells that fail to repolarize after a normal heartbeat. We investigated the proposal by representing the heart as an excitable medium of FitzHugh-Nagumo cells where a proportion of cells were capable of remaining depolarized indefinitely. As such, those cells exhibit bistable membrane dynamics. We found that heterogeneous media can tolerate a surprisingly large number of bistable cells and still support normal rhythmic activity. Yet there is a critical limit beyond which the medium is persistently arrhythmogenic. Numerical analysis revealed that the critical threshold for arrhythmogenesis depends on both the strength of the coupling between cells and the extent to which the abnormal cells resist repolarization. Moreover, arrhythmogenesis was found to emerge preferentially at tissue boundaries where cells naturally have fewer neighbors to influence their behavior. These findings may explain why atrial fibrillation typically originates from tissue boundaries such as the cuff of the pulmonary vein.

During sporogony, malaria-causing parasites infect a mosquito, reproduce and migrate to the mosquito salivary glands where they can be transmitted the next time blood feeding occurs. The time required for sporogony, known as the extrinsic incubation period (EIP), is an important determinant of malaria transmission intensity. The EIP is typically estimated as the time for a given percentile, x, of infected mosquitoes to develop salivary gland sporozoites (the infectious parasite life stage), which is denoted by EIPx. Many mechanisms, however, affect the observed sporozoite prevalence including the human-to-mosquito transmission probability and possibly differences in mosquito mortality according to infection status. To account for these various mechanisms, we present a mechanistic mathematical model, which explicitly models key processes at the parasite, mosquito and observational scales. Fitting this model to experimental data, we find greater variation in the EIP than previously thought: we estimated the range between EIP10 and EIP90 (at 27°C) as 4.5 days compared to 0.9 days using existing statistical methods. This pattern holds over the range of study temperatures included in the dataset. Increasing temperature from 21°C to 34°C decreased the EIP50 from 16.1 to 8.8 days. Our work highlights the importance of mechanistic modelling of sporogony to (1) improve estimates of malaria transmission under different environmental conditions or disease control programs and (2) evaluate novel interventions that target the mosquito life stages of the parasite.

Short QT (SQT) syndrome is a genetic cardiac disorder characterized by an abbreviated QT interval of the patient's electrocardiogram. The syndrome is associated with increased risk of arrhythmia and sudden cardiac death and can arise from a number of ion channel mutations. Cardiomyocytes derived from induced pluripotent stem cells generated from SQT patients (SQT hiPSC-CMs) provide promising platforms for testing pharmacological treatments directly in human cardiac cells exhibiting mutations specific for the syndrome. However, a difficulty is posed by the relative immaturity of hiPSC-CMs, with the possibility that drug effects observed in SQT hiPSC-CMs could be very different from the corresponding drug effect in vivo. In this paper, we apply a multistep computational procedure for translating measured drug effects from these cells to human QT response. This process first detects drug effects on individual ion channels based on measurements of SQT hiPSC-CMs and then uses these results to estimate the drug effects on ventricular action potentials and QT intervals of adult SQT patients. We find that the procedure is able to identify IC50 values in line with measured values for the four drugs quinidine, ivabradine, ajmaline and mexiletine. In addition, the predicted effect of quinidine on the adult QT interval is in good agreement with measured effects of quinidine for adult patients. Consequently, the computational procedure appears to be a useful tool for helping predicting adult drug responses from pure in vitro measurements of patient derived cell lines.

The voltage-dependent anion channel (VDAC) is a critical β-barrel membrane protein of the mitochondrial outer membrane, which regulates the transport of ions and ATP between mitochondria and the cytoplasm. In addition, VDAC plays a central role in the control of apoptosis and is therefore of great interest in both cancer and neurodegenerative diseases. Although not fully understood, it is presumed that the gating mechanism of VDAC is governed by its N-terminal region which, in the open state of the channel, exhibits an α-helical structure positioned midway inside the pore and strongly interacting with the β-barrel wall. In the present work, we performed molecular simulations with a recently developed force field for disordered systems to shed new light on known experimental results, showing that the N-terminus of VDAC is an intrinsically disordered region (IDR). First, simulation of the N-terminal segment as a free peptide highlighted its disordered nature and the importance of using an IDR-specific force field to properly sample its conformational landscape. Secondly, accelerated dynamics simulation of a double cysteine VDAC mutant under applied voltage revealed metastable low conducting states of the channel representative of closed states observed experimentally. Related structures were characterized by partial unfolding and rearrangement of the N-terminal tail, that led to steric hindrance of the pore. Our results indicate that the disordered properties of the N-terminus are crucial to properly account for the gating mechanism of VDAC.

Collective behavior is an emergent property of numerous complex systems, from financial markets to cancer cells to predator-prey ecological systems. Characterizing modes of collective behavior is often done through human observation, training generative models, or other supervised learning techniques. Each of these cases requires knowledge of and a method for characterizing the macro-state(s) of the system. This presents a challenge for studying novel systems where there may be little prior knowledge. Here, we present a new unsupervised method of detecting emergent behavior in complex systems, and discerning between distinct collective behaviors. We require only metrics, d(1), d(2), defined on the set of agents, X, which measure agents' nearness in variables of interest. We apply the method of diffusion maps to the systems (X, d(i)) to recover efficient embeddings of their interaction networks. Comparing these geometries, we formulate a measure of similarity between two networks, called the map alignment statistic (MAS). A large MAS is evidence that the two networks are codetermined in some fashion, indicating an emergent relationship between the metrics d(1) and d(2). Additionally, the form of the macro-scale organization is encoded in the covariances among the two sets of diffusion map components. Using these covariances we discern between different modes of collective behavior in a data-driven, unsupervised manner. This method is demonstrated on a synthetic flocking model as well as empirical fish schooling data. We show that our state classification subdivides the known behaviors of the school in a meaningful manner, leading to a finer description of the system's behavior.

Skilled behavior often displays signatures of Bayesian inference. In order for the brain to implement the required computations, neuronal activity must carry accurate information about the uncertainty of sensory inputs. Two major approaches have been proposed to study neuronal representations of uncertainty. The first one, the Bayesian decoding approach, aims primarily at decoding the posterior probability distribution of the stimulus from population activity using Bayes' rule, and indirectly yields uncertainty estimates as a by-product. The second one, which we call the correlational approach, searches for specific features of neuronal activity (such as tuning-curve width and maximum firing-rate) which correlate with uncertainty. To compare these two approaches, we derived a new normative model of sound source localization by Interaural Time Difference (ITD), that reproduces a wealth of behavioral and neural observations. We found that several features of neuronal activity correlated with uncertainty on average, but none provided an accurate estimate of uncertainty on a trial-by-trial basis, indicating that the correlational approach may not reliably identify which aspects of neuronal responses represent uncertainty. In contrast, the Bayesian decoding approach reveals that the activity pattern of the entire population was required to reconstruct the trial-to-trial posterior distribution with Bayes' rule. These results suggest that uncertainty is unlikely to be represented in a single feature of neuronal activity, and highlight the importance of using a Bayesian decoding approach when exploring the neural basis of uncertainty.

In this work, we introduce an entirely data-driven and automated approach to reveal disease-associated biomarker and risk factor networks from heterogeneous and high-dimensional healthcare data. Our workflow is based on Bayesian networks, which are a popular tool for analyzing the interplay of biomarkers. Usually, data require extensive manual preprocessing and dimension reduction to allow for effective learning of Bayesian networks. For heterogeneous data, this preprocessing is hard to automatize and typically requires domain-specific prior knowledge. We here combine Bayesian network learning with hierarchical variable clustering in order to detect groups of similar features and learn interactions between them entirely automated. We present an optimization algorithm for the adaptive refinement of such group Bayesian networks to account for a specific target variable, like a disease. The combination of Bayesian networks, clustering, and refinement yields low-dimensional but disease-specific interaction networks. These networks provide easily interpretable, yet accurate models of biomarker interdependencies. We test our method extensively on simulated data, as well as on data from the Study of Health in Pomerania (SHIP-TREND), and demonstrate its effectiveness using non-alcoholic fatty liver disease and hypertension as examples. We show that the group network models outperform available biomarker scores, while at the same time, they provide an easily interpretable interaction network.

Intratumour heterogeneity is increasingly recognized as a frequent problem for cancer treatment as it allows for the evolution of resistance against treatment. While cancer genotyping becomes more and more established and allows to determine the genetic heterogeneity, less is known about the phenotypic heterogeneity among cancer cells. We investigate how phenotypic differences can impact the efficiency of therapy options that select on this diversity, compared to therapy options that are independent of the phenotype. We employ the ecological concept of trait distributions and characterize the cancer cell population as a collection of subpopulations that differ in their growth rate. We show in a deterministic model that growth rate-dependent treatment types alter the trait distribution of the cell population, resulting in a delayed relapse compared to a growth rate-independent treatment. Whether the cancer cell population goes extinct or relapse occurs is determined by stochastic dynamics, which we investigate using a stochastic model. Again, we find that relapse is delayed for the growth rate-dependent treatment type, albeit an increased relapse probability, suggesting that slowly growing subpopulations are shielded from extinction. Sequential application of growth rate-dependent and growth rate-independent treatment types can largely increase treatment efficiency and delay relapse. Interestingly, even longer intervals between decisions to change the treatment type may achieve close-to-optimal efficiencies and relapse times. Monitoring patients at regular check-ups may thus provide the temporally resolved guidance to tailor treatments to the changing cancer cell trait distribution and allow clinicians to cope with this dynamic heterogeneity.

Drug combinations have demonstrated great potential in cancer treatments. They alleviate drug resistance and improve therapeutic efficacy. The fast-growing number of anti-cancer drugs has caused the experimental investigation of all drug combinations to become costly and time-consuming. Computational techniques can improve the efficiency of drug combination screening. Despite recent advances in applying machine learning to synergistic drug combination prediction, several challenges remain. First, the performance of existing methods is suboptimal. There is still much space for improvement. Second, biological knowledge has not been fully incorporated into the model. Finally, many models are lack interpretability, limiting their clinical applications. To address these challenges, we have developed a knowledge-enabled and self-attention transformer boosted deep learning model, TranSynergy, which improves the performance and interpretability of synergistic drug combination prediction. TranSynergy is designed so that the cellular effect of drug actions can be explicitly modeled through cell-line gene dependency, gene-gene interaction, and genome-wide drug-target interaction. A novel Shapley Additive Gene Set Enrichment Analysis (SA-GSEA) method has been developed to deconvolute genes that contribute to the synergistic drug combination and improve model interpretability. Extensive benchmark studies demonstrate that TranSynergy outperforms the state-of-the-art method, suggesting the potential of mechanism-driven machine learning. Novel pathways that are associated with the synergistic combinations are revealed and supported by experimental evidences. They may provide new insights into identifying biomarkers for precision medicine and discovering new anti-cancer therapies. Several new synergistic drug combinations have been predicted with high confidence for ovarian cancer which has few treatment options. The code is available at https://github.com/qiaoliuhub/drug_combination.

We present OpenAWSEM and Open3SPN2, new cross-compatible implementations of coarse-grained models for protein (AWSEM) and DNA (3SPN2) molecular dynamics simulations within the OpenMM framework. These new implementations retain the chemical accuracy and intrinsic efficiency of the original models while adding GPU acceleration and the ease of forcefield modification provided by OpenMM's Custom Forces software framework. By utilizing GPUs, we achieve around a 30-fold speedup in protein and protein-DNA simulations over the existing LAMMPS-based implementations running on a single CPU core. We showcase the benefits of OpenMM's Custom Forces framework by devising and implementing two new potentials that allow us to address important aspects of protein folding and structure prediction and by testing the ability of the combined OpenAWSEM and Open3SPN2 to model protein-DNA binding. The first potential is used to describe the changes in effective interactions that occur as a protein becomes partially buried in a membrane. We also introduced an interaction to describe proteins with multiple disulfide bonds. Using simple pairwise disulfide bonding terms results in unphysical clustering of cysteine residues, posing a problem when simulating the folding of proteins with many cysteines. We now can computationally reproduce Anfinsen's early Nobel prize winning experiments by using OpenMM's Custom Forces framework to introduce a multi-body disulfide bonding term that prevents unphysical clustering. Our protein-DNA simulations show that the binding landscape is funneled towards structures that are quite similar to those found using experiments. In summary, this paper provides a simulation tool for the molecular biophysics community that is both easy to use and sufficiently efficient to simulate large proteins and large protein-DNA systems that are central to many cellular processes. These codes should facilitate the interplay between molecular simulations and cellular studies, which have been hampered by the large mismatch between the time and length scales accessible to molecular simulations and those relevant to cell biology.

Erratum in PLoS Comput Biol. 2022 Oct 5;18(10):e1010592.

Analyzing the dynamical properties of mobile objects requires to extract trajectories from recordings, which is often done by tracking movies. We compiled a database of two-dimensional movies for very different biological and physical systems spanning a wide range of length scales and developed a general-purpose, optimized, open-source, cross-platform, easy to install and use, self-updating software called FastTrack. It can handle a changing number of deformable objects in a region of interest, and is particularly suitable for animal and cell tracking in two-dimensions. Furthermore, we introduce the probability of incursions as a new measure of a movie's trackability that doesn't require the knowledge of ground truth trajectories, since it is resilient to small amounts of errors and can be computed on the basis of an ad hoc tracking. We also leveraged the versatility and speed of FastTrack to implement an iterative algorithm determining a set of nearly-optimized tracking parameters-yet further reducing the amount of human intervention-and demonstrate that FastTrack can be used to explore the space of tracking parameters to optimize the number of swaps for a batch of similar movies. A benchmark shows that FastTrack is orders of magnitude faster than state-of-the-art tracking algorithms, with a comparable tracking accuracy. The source code is available under the GNU GPLv3 at https://github.com/FastTrackOrg/FastTrack and pre-compiled binaries for Windows, Mac and Linux are available at http://www.fasttrack.sh.

Prolonged lag time can be induced by starvation contributing to the antibiotic tolerance of bacteria. We analyze the optimal lag time to survive and grow the iterative and stochastic application of antibiotics. A simple model shows that the optimal lag time can exhibit a discontinuous transition when the severeness of the antibiotic application, such as the probability to be exposed the antibiotic, the death rate under the exposure, and the duration of the exposure, is increased. This suggests the possibility of reducing tolerant bacteria by controlled usage of antibiotics application. When the bacterial populations are able to have two phenotypes with different lag times, the fraction of the second phenotype that has different lag time shows a continuous transition. We then present a generic framework to investigate the optimal lag time distribution for total population fitness for a given distribution of the antibiotic application duration. The obtained optimal distributions have multiple peaks for a wide range of the antibiotic application duration distributions, including the case where the latter is monotonically decreasing. The analysis supports the advantage in evolving multiple, possibly discrete phenotypes in lag time for bacterial long-term fitness.

MIGNON is a workflow for the analysis of RNA-Seq experiments, which not only efficiently manages the estimation of gene expression levels from raw sequencing reads, but also calls genomic variants present in the transcripts analyzed. Moreover, this is the first workflow that provides a framework for the integration of transcriptomic and genomic data based on a mechanistic model of signaling pathway activities that allows a detailed biological interpretation of the results, including a comprehensive functional profiling of cell activity. MIGNON covers the whole process, from reads to signaling circuit activity estimations, using state-of-the-art tools, it is easy to use and it is deployable in different computational environments, allowing an optimized use of the resources available.

Pentameric ligand-gated ion channels (pLGICs) are receptor proteins that are sensitive to their membrane environment, but the mechanism for how lipids modulate function under physiological conditions in a state dependent manner is not known. The glycine receptor is a pLGIC whose structure has been resolved in different functional states. Using a realistic model of a neuronal membrane coupled with coarse-grained molecular dynamics simulations, we demonstrate that some key lipid-protein interactions are dependent on the receptor state, suggesting that lipids may regulate the receptor's conformational dynamics. Comparison with existing structural data confirms known lipid binding sites, but we also predict further protein-lipid interactions including a site at the communication interface between the extracellular and transmembrane domain. Moreover, in the active state, cholesterol can bind to the binding site of the positive allosteric modulator ivermectin. These protein-lipid interaction sites could in future be exploited for the rational design of lipid-like allosteric drugs.

One of the significant unanswered questions about COVID-19 epidemiology relates to the role of children in transmission. This study uses data on infections within households in order to estimate the susceptibility and infectivity of children compared to those of adults. The data were collected from households in the city of Bnei Brak, Israel, in which all household members were tested for COVID-19 using PCR (637 households, average household size of 5.3). In addition, serological tests were performed on a subset of the individuals in the study. Inspection of the PCR data shows that children are less likely to be tested positive compared to adults (25% of children positive over all households, 44% of adults positive over all households, excluding index cases), and the chance of being positive increases with age. Analysis of joint PCR/serological data shows that there is under-detection of infections in the PCR testing, which is more substantial in children. However, the differences in detection rates are not sufficient to account for the differences in PCR positive rates in the two age groups. To estimate relative transmission parameters, we employ a discrete stochastic model of the spread of infection within a household, allowing for susceptibility and infectivity parameters to differ among children and adults. The model is fitted to the household data using a simulated maximum likelihood approach. To adjust parameter estimates for under-detection of infections in the PCR results, we employ a multiple imputation procedure using estimates of under-detection in children and adults, based on the available serological data. We estimate that the susceptibility of children (under 20 years old) is 43% (95% CI: [31%, 55%]) of the susceptibility of adults. The infectivity of children was estimated to be 63% (95% CI: [37%, 88%]) relative to that of adults.

The spatial dynamics of epidemics are fundamentally affected by patterns of human mobility. Mobile phone call detail records (CDRs) are a rich source of mobility data, and allow semi-mechanistic models of movement to be parameterised even for resource-poor settings. While the gravity model typically reproduces human movement reasonably well at the administrative level spatial scale, past studies suggest that parameter estimates vary with the level of spatial discretisation at which models are fitted. Given that privacy concerns usually preclude public release of very fine-scale movement data, such variation would be problematic for individual-based simulations of epidemic spread parametrised at a fine spatial scale. We therefore present new methods to fit fine-scale mathematical mobility models (here we implement variants of the gravity and radiation models) to spatially aggregated movement data and investigate how model parameter estimates vary with spatial resolution. We use gridded population data at 1km resolution to derive population counts at different spatial scales (down to ∼ 5km grids) and implement mobility models at each scale. Parameters are estimated from administrative-level flow data between overnight locations in Kenya and Namibia derived from CDRs: where the model spatial resolution exceeds that of the mobility data, we compare the flow data between a particular origin and destination with the sum of all model flows between cells that lie within those particular origin and destination administrative units. Clear evidence of over-dispersion supports the use of negative binomial instead of Poisson likelihood for count data with high values. Radiation models use fewer parameters than the gravity model and better predict trips between overnight locations for both considered countries. Results show that estimates for some parameters change between countries and with spatial resolution and highlight how imperfect flow data and spatial population distribution can influence model fit.

The correct identification of metabolic activity in tissues or cells under different conditions can be extremely elusive due to mechanisms such as post-transcriptional modification of enzymes or different rates in protein degradation, making difficult to perform predictions on the basis of gene expression alone. Context-specific metabolic network reconstruction can overcome some of these limitations by leveraging the integration of multi-omics data into genome-scale metabolic networks (GSMN). Using the experimental information, context-specific models are reconstructed by extracting from the generic GSMN the sub-network most consistent with the data, subject to biochemical constraints. One advantage is that these context-specific models have more predictive power since they are tailored to the specific tissue, cell or condition, containing only the reactions predicted to be active in such context. However, an important limitation is that there are usually many different sub-networks that optimally fit the experimental data. This set of optimal networks represent alternative explanations of the possible metabolic state. Ignoring the set of possible solutions reduces the ability to obtain relevant information about the metabolism and may bias the interpretation of the true metabolic states. In this work we formalize the problem of enumerating optimal metabolic networks and we introduce DEXOM, an unified approach for diversity-based enumeration of context-specific metabolic networks. We developed different strategies for this purpose and we performed an exhaustive analysis using simulated and real data. In order to analyze the extent to which these results are biologically meaningful, we used the alternative solutions obtained with the different methods to measure: 1) the improvement of in silico predictions of essential genes in Saccharomyces cerevisiae using ensembles of metabolic network; and 2) the detection of alternative enriched pathways in different human cancer cell lines. We also provide DEXOM as an open-source library compatible with COBRA Toolbox 3.0, available at https://github.com/MetExplore/dexom.

Recent work of Sottoriva, Graham, and collaborators have led to the controversial claim that exponentially growing tumors have a site frequency spectrum that follows the 1/f law consistent with neutral evolution. This conclusion has been criticized based on data quality issues, statistical considerations, and simulation results. Here, we use rigorous mathematical arguments to investigate the site frequency spectrum in the two-type model of clonal evolution. If the fitnesses of the two types are λ0 < λ1, then the site frequency spectrum is c/fα where α = λ0/λ1. This is due to the advantageous mutations that produce the founders of the type 1 population. Mutations within the growing type 0 and type 1 populations follow the 1/f law. Our results show that, in contrast to published criticisms, neutral evolution in an exponentially growing tumor can be distinguished from the two-type model using the site frequency spectrum.

Optogenetic targeting of astrocytes provides a robust experimental model to differentially induce Ca2+ signals in astrocytes in vivo. However, a systematic study quantifying the response of optogenetically modified astrocytes to light is yet to be performed. Here, we propose a novel stochastic model of Ca2+ dynamics in astrocytes that incorporates a light sensitive component-channelrhodopsin 2 (ChR2). Utilizing this model, we investigated the effect of different light stimulation paradigms on cells expressing select variants of ChR2 (wild type, ChETA, and ChRET/TC). Results predict that depending on paradigm specification, astrocytes might undergo drastic changes in their basal Ca2+ level and spiking probability. Furthermore, we performed a global sensitivity analysis to assess the effect of variation in parameters pertinent to the shape of the ChR2 photocurrent on astrocytic Ca2+ dynamics. Results suggest that directing variants towards the first open state of the ChR2 photocycle (o1) enhances spiking activity in astrocytes during optical stimulation. Evaluation of the effect of Ca2+ buffering and coupling coefficient in a network of ChR2-expressing astrocytes demonstrated basal level elevations in the stimulated region and propagation of calcium activity to unstimulated cells. Buffering reduced the diffusion range of Ca2+ within the network, thereby limiting propagation and influencing the activity of astrocytes. Collectively, the framework presented in this study provides valuable information for the selection of light stimulation paradigms that elicit desired astrocytic activity using existing ChR2 constructs, as well as aids in the engineering of future application-oriented optogenetic variants.

Pavlovian associations drive approach towards reward-predictive cues, and avoidance of punishment-predictive cues. These associations "misbehave" when they conflict with correct instrumental behavior. This raises the question of how Pavlovian and instrumental influences on behavior are arbitrated. We test a computational theory according to which Pavlovian influence will be stronger when inferred controllability of outcomes is low. Using a model-based analysis of a Go/NoGo task with human subjects, we show that theta-band oscillatory power in frontal cortex tracks inferred controllability, and that these inferences predict Pavlovian action biases. Functional MRI data revealed an inferior frontal gyrus correlate of action probability and a ventromedial prefrontal correlate of outcome valence, both of which were modulated by inferred controllability.

Increasingly complex in silico modeling approaches offer a way to simultaneously access cancerous processes at different spatio-temporal scales. High-level models, such as those based on partial differential equations, are computationally affordable and allow large tumor sizes and long temporal windows to be studied, but miss the discrete nature of many key underlying cellular processes. Individual-based approaches provide a much more detailed description of tumors, but have difficulties when trying to handle full-sized real cancers. Thus, there exists a trade-off between the integration of macroscopic and microscopic information, now widely available, and the ability to attain clinical tumor sizes. In this paper we put forward a stochastic mesoscopic simulation framework that incorporates key cellular processes during tumor progression while keeping computational costs to a minimum. Our framework captures a physical scale that allows both the incorporation of microscopic information, tracking the spatio-temporal emergence of tumor heterogeneity and the underlying evolutionary dynamics, and the reconstruction of clinically sized tumors from high-resolution medical imaging data, with the additional benefit of low computational cost. We illustrate the functionality of our modeling approach for the case of glioblastoma, a paradigm of tumor heterogeneity that remains extremely challenging in the clinical setting.

The dynamics of cerebellar neuronal networks is controlled by the underlying building blocks of neurons and synapses between them. For which, the computation of Purkinje cells (PCs), the only output cells of the cerebellar cortex, is implemented through various types of neural pathways interactively routing excitation and inhibition converged to PCs. Such tuning of excitation and inhibition, coming from the gating of specific pathways as well as short-term plasticity (STP) of the synapses, plays a dominant role in controlling the PC dynamics in terms of firing rate and spike timing. PCs receive cascade feedforward inputs from two major neural pathways: the first one is the feedforward excitatory pathway from granule cells (GCs) to PCs; the second one is the feedforward inhibition pathway from GCs, via molecular layer interneurons (MLIs), to PCs. The GC-PC pathway, together with short-term dynamics of excitatory synapses, has been a focus over past decades, whereas recent experimental evidence shows that MLIs also greatly contribute to controlling PC activity. Therefore, it is expected that the diversity of excitation gated by STP of GC-PC synapses, modulated by strong inhibition from MLI-PC synapses, can promote the computation performed by PCs. However, it remains unclear how these two neural pathways are interacted to modulate PC dynamics. Here using a computational model of PC network installed with these two neural pathways, we addressed this question to investigate the change of PC firing dynamics at the level of single cell and network. We show that the nonlinear characteristics of excitatory STP dynamics can significantly modulate PC spiking dynamics mediated by inhibition. The changes in PC firing rate, firing phase, and temporal spike pattern, are strongly modulated by these two factors in different ways. MLIs mainly contribute to variable delays in the postsynaptic action potentials of PCs while modulated by excitation STP. Notably, the diversity of synchronization and pause response in the PC network is governed not only by the balance of excitation and inhibition, but also by the synaptic STP, depending on input burst patterns. Especially, the pause response shown in the PC network can only emerge with the interaction of both pathways. Together with other recent findings, our results show that the interaction of feedforward pathways of excitation and inhibition, incorporated with synaptic short-term dynamics, can dramatically regulate the PC activities that consequently change the network dynamics of the cerebellar circuit.

Different miRNA profiling protocols and technologies introduce differences in the resulting quantitative expression profiles. These include differences in the presence (and measurability) of certain miRNAs. We present and examine a method based on quantile normalization, Adjusted Quantile Normalization (AQuN), to combine miRNA expression data from multiple studies in breast cancer into a single joint dataset for integrative analysis. By pooling multiple datasets, we obtain increased statistical power, surfacing patterns that do not emerge as statistically significant when separately analyzing these datasets. To merge several datasets, as we do here, one needs to overcome both technical and batch differences between these datasets. We compare several approaches for merging and jointly analyzing miRNA datasets. We investigate the statistical confidence for known results and highlight potential new findings that resulted from the joint analysis using AQuN. In particular, we detect several miRNAs to be differentially expressed in estrogen receptor (ER) positive versus ER negative samples. In addition, we identify new potential biomarkers and therapeutic targets for both clinical groups. As a specific example, using the AQuN-derived dataset we detect hsa-miR-193b-5p to have a statistically significant over-expression in the ER positive group, a phenomenon that was not previously reported. Furthermore, as demonstrated by functional assays in breast cancer cell lines, overexpression of hsa-miR-193b-5p in breast cancer cell lines resulted in decreased cell viability in addition to inducing apoptosis. Together, these observations suggest a novel functional role for this miRNA in breast cancer. Packages implementing AQuN are provided for Python and Matlab: https://github.com/YakhiniGroup/PyAQN.

Mathematical disease modelling has long operated under the assumption that any one infectious disease is caused by one transmissible pathogen spreading among a population. This paradigm has been useful in simplifying the biological reality of epidemics and has allowed the modelling community to focus on the complexity of other factors such as population structure and interventions. However, there is an increasing amount of evidence that the strain diversity of pathogens, and their interplay with the host immune system, can play a large role in shaping the dynamics of epidemics. Here, we introduce a disease model with an underlying genotype network to account for two important mechanisms. One, the disease can mutate along network pathways as it spreads in a host population. Two, the genotype network allows us to define a genetic distance between strains and therefore to model the transcendence of immunity often observed in real world pathogens. We study the emergence of epidemics in this model, through its epidemic phase transitions, and highlight the role of the genotype network in driving cyclicity of diseases, large scale fluctuations, sequential epidemic transitions, as well as localization around specific strains of the associated pathogen. More generally, our model illustrates the richness of behaviours that are possible even in well-mixed host populations once we consider strain diversity and go beyond the "one disease equals one pathogen" paradigm.

Epidemics may pose a significant dilemma for governments and individuals. The personal or public health consequences of inaction may be catastrophic; but the economic consequences of drastic response may likewise be catastrophic. In the face of these trade-offs, governments and individuals must therefore strike a balance between the economic and personal health costs of reducing social contacts and the public health costs of neglecting to do so. As risk of infection increases, potentially infectious contact between people is deliberately reduced either individually or by decree. This must be balanced against the social and economic costs of having fewer people in contact, and therefore active in the labor force or enrolled in school. Although the importance of adaptive social contact on epidemic outcomes has become increasingly recognized, the most important properties of coupled human-natural epidemic systems are still not well understood. We develop a theoretical model for adaptive, optimal control of the effective social contact rate using traditional epidemic modeling tools and a utility function with delayed information. This utility function trades off the population-wide contact rate with the expected cost and risk of increasing infections. Our analytical and computational analysis of this simple discrete-time deterministic strategic model reveals the existence of an endemic equilibrium, oscillatory dynamics around this equilibrium under some parametric conditions, and complex dynamic regimes that shift under small parameter perturbations. These results support the supposition that infectious disease dynamics under adaptive behavior change may have an indifference point, may produce oscillatory dynamics without other forcing, and constitute complex adaptive systems with associated dynamics. Implications for any epidemic in which adaptive behavior influences infectious disease dynamics include an expectation of fluctuations, for a considerable time, around a quasi-equilibrium that balances public health and economic priorities, that shows multiple peaks and surges in some scenarios, and that implies a high degree of uncertainty in mathematical projections.

Many viral infections can be prevented by immunizing with live, attenuated vaccines. Early methods of attenuation were hit-and-miss, now much improved by genetic engineering. However, even current methods operate on the principle of genetic harm, reducing the virus's ability to grow. Reduced viral growth has the undesired side-effect of reducing the host immune response below that of infection with wild-type. Might some methods of attenuation instead lead to an increased immune response? We use mathematical models of the dynamics of virus with innate and adaptive immunity to explore the tradeoff between attenuation of virus pathology and immunity. We find that modification of some virus immune-evasion pathways can indeed reduce pathology yet enhance immunity. Thus, attenuated vaccines can, in principle, be directed to be safe yet create better immunity than is elicited by the wild-type virus.

Newly-translated glycoproteins in the endoplasmic reticulum (ER) often undergo cycles of chaperone binding and release in order to assist in folding. Quality control is required to distinguish between proteins that have completed native folding, those that have yet to fold, and those that have misfolded. Using quantitative modeling, we explore how the design of the quality-control pathway modulates its efficiency. Our results show that an energy-consuming cyclic quality-control process, similar to the observed physiological system, outperforms alternative designs. The kinetic parameters that optimize the performance of this system drastically change with protein production levels, while remaining relatively insensitive to the protein folding rate. Adjusting only the degradation rate, while fixing other parameters, allows the pathway to adapt across a range of protein production levels, aligning with in vivo measurements that implicate the release of degradation-associated enzymes as a rapid-response system for perturbations in protein homeostasis. The quantitative models developed here elucidate design principles for effective glycoprotein quality control in the ER, improving our mechanistic understanding of a system crucial to maintaining cellular health.

Nuclear Magnetic Resonance (NMR) spectroscopy is one of the three primary experimental means of characterizing macromolecular structures, including protein structures. Structure determination by solution NMR spectroscopy has traditionally relied heavily on distance restraints derived from nuclear Overhauser effect (NOE) measurements. While structure determination of proteins from NOE-based restraints is well understood and broadly used, structure determination from Residual Dipolar Couplings (RDCs) is relatively less well developed. Here, we describe the new features of the protein structure modeling program REDCRAFT and focus on the new Adaptive Decimation (AD) feature. The AD plays a critical role in improving the robustness of REDCRAFT to missing or noisy data, while allowing structure determination of larger proteins from less data. In this report we demonstrate the successful application of REDCRAFT in structure determination of proteins ranging in size from 50 to 145 residues using experimentally collected data, and of larger proteins (145 to 573 residues) using simulated RDC data. In both cases, REDCRAFT uses only RDC data that can be collected from perdeuterated proteins. Finally, we compare the accuracy of structure determination from RDCs alone with traditional NOE-based methods for the structurally novel PF.2048.1 protein. The RDC-based structure of PF.2048.1 exhibited 1.0 Å BB-RMSD with respect to a high-quality NOE-based structure. Although optimal strategies would include using RDC data together with chemical shift, NOE, and other NMR data, these studies provide proof-of-principle for robust structure determination of largely-perdeuterated proteins from RDC data alone using REDCRAFT.

RNA is considered as an attractive target for new small molecule drugs. Designing active compounds can be facilitated by computational modeling. Most of the available tools developed for these prediction purposes, such as molecular docking or scoring functions, are parametrized for protein targets. The performance of these methods, when applied to RNA-ligand systems, is insufficient. To overcome these problems, we developed AnnapuRNA, a new knowledge-based scoring function designed to evaluate RNA-ligand complex structures, generated by any computational docking method. We also evaluated three main factors that may influence the structure prediction, i.e., the starting conformer of a ligand, the docking program, and the scoring function used. We applied the AnnapuRNA method for a post-hoc study of the recently published structures of the FMN riboswitch. Software is available at https://github.com/filipspl/AnnapuRNA.

While large datasets of HIV-1 sequences are increasingly being generated, many studies rely on a single gene or fragment of the genome and few comparative studies across genes have been done. We performed genome-based and gene-specific Bayesian phylogenetic analyses to investigate how certain factors impact estimates of the infection dates in an acute HIV-1 infection cohort, RV217. In this cohort, HIV-1 diagnosis corresponded to the first RNA positive test and occurred a median of four days after the last negative test, allowing us to compare timing estimates using BEAST to a narrow window of infection. We analyzed HIV-1 sequences sampled one week, one month and six months after HIV-1 diagnosis in 39 individuals. We found that shared diversity and temporal signal was limited in acute infection, and insufficient to allow timing inferences in the shortest HIV-1 genes, thus dated phylogenies were primarily analyzed for env, gag, pol and near full-length genomes. There was no one best-fitting model across participants and genes, though relaxed molecular clocks (73% of best-fitting models) and the Bayesian skyline (49%) tended to be favored. For infections with single founders, the infection date was estimated to be around one week pre-diagnosis for env (IQR: 3-9 days) and gag (IQR: 5-9 days), whilst the genome placed it at a median of 10 days (IQR: 4-19). Multiply-founded infections proved problematic to date. Our ability to compare timing inferences to precise estimates of HIV-1 infection (within a week) highlights that molecular dating methods can be applied to within-host datasets from early infection. Nonetheless, our results also suggest caution when using uniform clock and population models or short genes with limited information content.