Search

Computational modeling has revealed that human research participants use both rapid working memory (WM) and incremental reinforcement learning (RL) (RL+WM) to solve a simple instrumental learning task, relying on WM when the number of stimuli is small and supplementing with RL when the number of stimuli exceeds WM capacity. Inspired by this work, we examined which learning systems and strategies are used by adolescent and adult mice when they first acquire a conditional associative learning task. In a version of the human RL+WM task translated for rodents, mice were required to associate odor stimuli (from a set of 2 or 4 odors) with a left or right port to receive reward. Using logistic regression and computational models to analyze the first 200 trials per odor, we determined that mice used both incremental RL and stimulus-insensitive, one-back strategies to solve the task. While these one-back strategies may be a simple form of short-term or working memory, they did not approximate the boost to learning performance that has been observed in human participants using WM in a comparable task. Adolescent and adult mice also showed comparable performance, with no change in learning rate or softmax beta parameters with adolescent development and task experience. However, reliance on a one-back perseverative, win-stay strategy increased with development in males in both odor set sizes, but was not dependent on gonadal hormones. Our findings advance a simple conditional associative learning task and new models to enable the isolation and quantification of reinforcement learning alongside other strategies mice use while learning to associate stimuli with rewards within a single behavioral session. These data and methods can inform and aid comparative study of reinforcement learning across species.

Symbolic regression with polynomial neural networks and polynomial neural ordinary differential equations (ODEs) are two recent and powerful approaches for equation recovery of many science and engineering problems. However, these methods provide point estimates for the model parameters and are currently unable to accommodate noisy data. We address this challenge by developing and validating the following Bayesian inference methods: the Laplace approximation, Markov Chain Monte Carlo (MCMC) sampling methods, and variational inference. We have found the Laplace approximation to be the best method for this class of problems. Our work can be easily extended to the broader class of symbolic neural networks to which the polynomial neural network belongs.

Monosynaptic connectivity mapping is crucial for building circuit-level models of neural computation. Two-photon optogenetic stimulation, when combined with whole-cell recording, enables large-scale mapping of physiological circuit parameters. In this experimental setup, recorded postsynaptic currents are used to infer the presence and strength of connections. For many cell types, nearby connections are those we expect to be strongest. However, when the postsynaptic cell expresses opsin, optical excitation of nearby cells can induce direct photocurrents in the postsynaptic cell. These photocurrent artifacts contaminate synaptic currents, making it difficult or impossible to probe connectivity for nearby cells. To overcome this problem, we developed a computational tool, Photocurrent Removal with Constraints (PhoRC). Our method is based on a constrained matrix factorization model which leverages the fact that photocurrent kinetics are less variable than those of synaptic currents. We demonstrate on real and simulated data that PhoRC consistently removes photocurrents while preserving synaptic currents, despite variations in photocurrent kinetics across datasets. Our method allows the discovery of synaptic connections which would have been otherwise obscured by photocurrent artifacts, and may thus reveal a more complete picture of synaptic connectivity. PhoRC runs faster than real time and is available as open source software.

The self-organization of cells relies on the profound complexity of protein-protein interactions. Challenges in directly observing these events have hindered progress toward understanding their diverse behaviors. One notable example is the interaction between molecular motors and cytoskeletal systems that combine to perform a variety of cellular functions. In this work, we leverage theory and experiments to identify and quantify the rate-limiting mechanism of the initial association between a cargo-bound kinesin motor and a microtubule track. Recent advances in optical tweezers provide binding times for several lengths of kinesin motors trapped at varying distances from a microtubule, empowering the investigation of competing models. We first explore a diffusion-limited model of binding. Through Brownian dynamics simulations and simulation-based inference, we find this simple diffusion model fails to explain the experimental binding times, but an extended model that accounts for the ADP state of the molecular motor agrees closely with the data, even under the scrutiny of penalizing for additional model complexity. We provide quantification of both kinetic rates and biophysical parameters underlying the proposed binding process. Our model suggests that a typical binding event is limited by ADP state rather than physical search. Lastly, we predict how these association rates can be modulated in distinct ways through variation of environmental concentrations and physical properties.

Implicit adaptation has been regarded as a rigid process that automatically operates in response to movement errors to keep the sensorimotor system precisely calibrated. This hypothesis has been challenged by recent evidence suggesting flexibility in this learning process. One compelling line of evidence comes from work suggesting that this form of learning is context-dependent, with the rate of learning modulated by error history. Specifically, learning was attenuated in the presence of perturbations exhibiting high variance compared to when the perturbation is fixed. However, these findings are confounded by the fact that the adaptation system corrects for errors of different magnitudes in a non-linear manner, with the adaptive response increasing in a proportional manner to small errors and saturating to large errors. Through simulations, we show that this non-linear motor correction function is sufficient to explain the effect of perturbation variance without referring to an experience-dependent change in error sensitivity. Moreover, by controlling the distribution of errors experienced during training, we provide empirical evidence showing that there is no measurable effect of perturbation variance on implicit adaptation. As such, we argue that the evidence to date remains consistent with the rigidity assumption.

Behavioral speech tasks have been widely used to understand the mechanisms of speech motor control in typical speakers as well as in various clinical populations. However, determining which neural functions differ between typical speakers and clinical populations based on behavioral data alone is difficult because multiple mechanisms may lead to the same behavioral differences. For example, individuals with cerebellar ataxia (CA) produce atypically large compensatory responses to pitch perturbations in their auditory feedback, compared to typical speakers, but this pattern could have many explanations. Here, computational modeling techniques were used to address this challenge. Bayesian inference was used to fit a state feedback control (SFC) model of voice fundamental frequency (fo) control to the behavioral pitch perturbation responses of speakers with CA and typical speakers. This fitting process resulted in estimates of posterior likelihood distributions for five model parameters (sensory feedback delays, absolute and relative levels of auditory and somatosensory feedback noise, and controller gain), which were compared between the two groups. Results suggest that the speakers with CA may proportionally weight auditory and somatosensory feedback differently from typical speakers. Specifically, the CA group showed a greater relative sensitivity to auditory feedback than the control group. There were also large group differences in the controller gain parameter, suggesting increased motor output responses to target errors in the CA group. These modeling results generate hypotheses about how CA may affect the speech motor system, which could help guide future empirical investigations in CA. This study also demonstrates the overall proof-of-principle of using this Bayesian inference approach to understand behavioral speech data in terms of interpretable parameters of speech motor control models.

Understanding the mechanisms of the cellular aging processes is crucial for attempting to extend organismal lifespan and for studying age-related degenerative diseases. Yeast cells divide through budding, providing a classical biological model for studying cellular aging. With their powerful genetics, relatively short cell cycle, and well-established signaling pathways also found in animals, yeast cells offer valuable insights into the aging process. Recent experiments suggested the existence of two aging modes in yeast characterized by nucleolar and mitochondrial declines, respectively. By analyzing experimental data, this study shows that cells evolving into those two aging modes behave differently when they are young. While buds grow linearly in both modes, cells that consistently generate spherical buds throughout their lifespan demonstrate greater efficacy in controlling bud size and growth rate at young ages. A three-dimensional multiscale chemical-mechanical model was developed and used to suggest and test hypothesized impacts of aging on bud morphogenesis. Experimentally calibrated model simulations showed that during the early stage of budding, tubular bud shape in one aging mode could be generated by locally inserting new materials at the bud tip, a process guided by the polarized Cdc42 signal. Furthermore, the aspect ratio of the tubular bud could be stabilized during the late stage as observed in experiments in this work. The model simulation results suggest that the localization of new cell surface material insertion, regulated by chemical signal polarization, could be weakened due to cellular aging in yeast and other cell types, leading to the change and stabilization of the bud aspect ratio.

Computational cognitive models have been used extensively to formalize cognitive processes. Model parameters offer a simple way to quantify individual differences in how humans process information. Similarly, model comparison allows researchers to identify which theories, embedded in different models, provide the best accounts of the data. Cognitive modeling uses statistical tools to quantitatively relate models to data that often rely on computing/estimating the likelihood of the data under the model. However, this likelihood is computationally intractable for a substantial number of models. These relevant models may embody reasonable theories of cognition, but are often under-explored due to the limited range of tools available to relate them to data. We contribute to filling this gap in a simple way using artificial neural networks (ANNs) to map data directly onto model identity and parameters, bypassing the likelihood estimation. We test our instantiation of an ANN as a cognitive model fitting tool on classes of cognitive models with strong inter-trial dependencies (such as reinforcement learning models), which offer unique challenges to most methods. We show that we can adequately perform both parameter estimation and model identification using our ANN approach, including for models that cannot be fit using traditional likelihood-based methods. We further discuss our work in the context of the ongoing research leveraging simulation-based approaches to parameter estimation and model identification, and how these approaches broaden the class of cognitive models researchers can quantitatively investigate.

Contrast transport models are widely used to quantify blood flow and transport in dynamic contrast-enhanced magnetic resonance imaging. These models analyze the time course of the contrast agent concentration, providing diagnostic and prognostic value for many biological systems. Thus, ensuring accuracy and repeatability of the model parameter estimation is a fundamental concern. In this work, we analyze the structural and practical identifiability of a class of nested compartment models pervasively used in analysis of MRI data. We combine artificial and real data to study the role of noise in model parameter estimation. We observe that although all the models are structurally identifiable, practical identifiability strongly depends on the data characteristics. We analyze the impact of increasing data noise on parameter identifiability and show how the latter can be recovered with increased data quality. To complete the analysis, we show that the results do not depend on specific tissue characteristics or the type of enhancement patterns of contrast agent signal.

The brain produces diverse functions, from perceiving sounds to producing arm reaches, through the collective activity of populations of many neurons. Determining if and how the features of these exogenous variables (e.g., sound frequency, reach angle) are reflected in population neural activity is important for understanding how the brain operates. Often, high-dimensional neural population activity is confined to low-dimensional latent spaces. However, many current methods fail to extract latent spaces that are clearly structured by exogenous variables. This has contributed to a debate about whether or not brains should be thought of as dynamical systems or representational systems. Here, we developed a new latent process Bayesian regression framework, the orthogonal stochastic linear mixing model (OSLMM) which introduces an orthogonality constraint amongst time-varying mixture coefficients, and provide Markov chain Monte Carlo inference procedures. We demonstrate superior performance of OSLMM on latent trajectory recovery in synthetic experiments and show superior computational efficiency and prediction performance on several real-world benchmark data sets. We primarily focus on demonstrating the utility of OSLMM in two neural data sets: μECoG recordings from rat auditory cortex during presentation of pure tones and multi-single unit recordings form monkey motor cortex during complex arm reaching. We show that OSLMM achieves superior or comparable predictive accuracy of neural data and decoding of external variables (e.g., reach velocity). Most importantly, in both experimental contexts, we demonstrate that OSLMM latent trajectories directly reflect features of the sounds and reaches, demonstrating that neural dynamics are structured by neural representations. Together, these results demonstrate that OSLMM will be useful for the analysis of diverse, large-scale biological time-series datasets.

We combine an off-lattice agent-based mathematical model and experimentation to explore filamentous growth of a yeast colony. Under environmental stress, Saccharomyces cerevisiae yeast cells can transition from a bipolar (sated) to unipolar (pseudohyphal) budding mechanism, where cells elongate and bud end-to-end. This budding asymmetry yields spatially non-uniform growth, where filaments extend away from the colony centre, foraging for food. We use approximate Bayesian computation to quantify how individual cell budding mechanisms give rise to spatial patterns observed in experiments. We apply this method of parameter inference to experimental images of colonies of two strains of S. cerevisiae, in low and high nutrient environments. The colony size at the transition from sated to pseudohyphal growth, and a forking mechanism for pseudohyphal cell proliferation are the key features driving colony morphology. Simulations run with the most likely inferred parameters produce colony morphologies that closely resemble experimental results. Author summary: Yeasts are one of the most-studied organisms in biology due to their widespread use in food and beverage production and their role as a model organism in biomedical research. In this work, we combine mathematical modelling with experimentation to better understand the growth mechanisms of a yeast colony. Typically, yeast cells are circular, but under environmental stress new daughter cells can become elongated. This change in the cell shape creates non-uniform growth in the colony, where filaments consisting of chains of elongated cells extend away from the colony centre, foraging for food. We use computer simulations to understand how individual cell shapes and reproduction mechanisms produce different colony patterns observed in experiments. The colony size at the sated-to-pseudohyphal transition and a forking in the cell proliferation mechanism were the key features that determine the shape and structure of a colony. We showed that by using the most likely set of parameters, our computer simulation is able to generate colonies resembling those seen in experiments. [ABSTRACT FROM AUTHOR]

Reducing spillover of zoonotic pathogens is an appealing approach to preventing human disease and minimizing the risk of future epidemics and pandemics. Although the immediate human health benefit of reducing spillover is clear, over time, spillover reduction could lead to counterintuitive negative consequences for human health. Here, we use mathematical models and computer simulations to explore the conditions under which unanticipated consequences of spillover reduction can occur in systems where the severity of disease increases with age at infection. Our results demonstrate that, because the average age at infection increases as spillover is reduced, programs that reduce spillover can actually increase population-level disease burden if the clinical severity of infection increases sufficiently rapidly with age. If, however, immunity wanes over time and reinfection is possible, our results reveal that negative health impacts of spillover reduction become substantially less likely. When our model is parameterized using published data on Lassa virus in West Africa, it predicts that negative health outcomes are possible, but likely to be restricted to a small subset of populations where spillover is unusually intense. Together, our results suggest that adverse consequences of spillover reduction programs are unlikely but that the public health gains observed immediately after spillover reduction may fade over time as the age structure of immunity gradually re-equilibrates to a reduced force of infection. Author summary: Many pathogens, such as rabies, coronaviruses, and hantaviruses primarily circulate within wild animals but can infect humans when the opportunity arises. This pervasive challenge to public health has motivated the development of new methods designed to reduce the frequency of these spillover events. Although reducing spillover infection of humans appears to be an obvious win for public health, it is conceivable that altering historical patterns of spillover could change the age structure of human immunity in a way that undermines human health. Using mathematical and computational models we evaluate the conditions required for these counterintuitive impacts to occur. Our analyses demonstrate that reducing spillover from wild animals will generally improve public health and that negative outcomes can occur in only rare and unusual circumstances. Although negative impacts of spillover reduction are likely to be rare, our results show that the public health benefits of spillover reduction may fade over time unless a near total elimination of spillover can be achieved. [ABSTRACT FROM AUTHOR]

Mitochondria form a network in the cell that rapidly changes through fission, fusion, and motility. Dysregulation of this four-dimensional (4D: x,y,z,time) network is implicated in numerous diseases ranging from cancer to neurodegeneration. While lattice light-sheet microscopy has recently made it possible to image mitochondria in 4D, quantitative analysis methods for the resulting datasets have been lacking. Here we present MitoTNT, the first-in-class software for Mitochondrial Temporal Network Tracking in 4D live-cell fluorescence microscopy data. MitoTNT uses spatial proximity and network topology to compute an optimal tracking assignment. To validate the accuracy of tracking, we created a reaction-diffusion simulation to model mitochondrial network motion and remodeling events. We found that our tracking is >90% accurate for ground-truth simulations and agrees well with published motility results for experimental data. We used MitoTNT to quantify 4D mitochondrial networks from human induced pluripotent stem cells. First, we characterized sub-fragment motility and analyzed network branch motion patterns. We revealed that the skeleton node motion is correlated along branch nodes and is uncorrelated in time. Second, we identified fission and fusion events with high spatiotemporal resolution. We found that mitochondrial skeleton nodes near the fission/fusion sites move nearly twice as fast as random skeleton nodes and that microtubules play a role in mediating selective fission/fusion. Finally, we developed graph-based transport simulations that model how material would distribute on experimentally measured mitochondrial temporal networks. We showed that pharmacological perturbations increase network reachability but decrease network resilience through a combination of altered mitochondrial fission/fusion dynamics and motility. MitoTNT's easy-to-use tracking module, interactive 4D visualization capability, and powerful post-tracking analyses aim at making temporal network tracking accessible to the wider mitochondria research community.

Sequence simulators are fundamental tools in bioinformatics, as they allow us to test data processing and inference tools, and are an essential component of some inference methods. The ongoing surge in available sequence data is however testing the limits of our bioinformatics software. One example is the large number of SARS-CoV-2 genomes available, which are beyond the processing power of many methods, and simulating such large datasets is also proving difficult. Here, we present a new algorithm and software for efficiently simulating sequence evolution along extremely large trees (e.g. > 100, 000 tips) when the branches of the tree are short, as is typical in genomic epidemiology. Our algorithm is based on the Gillespie approach, and it implements an efficient multi-layered search tree structure that provides high computational efficiency by taking advantage of the fact that only a small proportion of the genome is likely to mutate at each branch of the considered phylogeny. Our open source software allows easy integration with other Python packages as well as a variety of evolutionary models, including indel models and new hypermutability models that we developed to more realistically represent SARS-CoV-2 genome evolution.

Identifying host factors that influence infectious disease transmission is an important step toward developing interventions to reduce disease incidence. Recent advances in methods for reconstructing infectious disease transmission events using pathogen genomic and epidemiological data open the door for investigation of host factors that affect onward transmission. While most transmission reconstruction methods are designed to work with densely sampled outbreaks, these methods are making their way into surveillance studies, where the fraction of sampled cases with sequenced pathogens could be relatively low. Surveillance studies that use transmission event reconstruction then use the reconstructed events as response variables (i.e., infection source status of each sampled case) and use host characteristics as predictors (e.g., presence of HIV infection) in regression models. We use simulations to study estimation of the effect of a host factor on probability of being an infection source via this multi-step inferential procedure. Using TransPhylo-a widely-used method for Bayesian estimation of infectious disease transmission events-and logistic regression, we find that low sensitivity of identifying infection sources leads to dilution of the signal, biasing logistic regression coefficients toward zero. We show that increasing the proportion of sampled cases improves sensitivity and some, but not all properties of the logistic regression inference. Application of these approaches to real world data from a population-based TB study in Botswana fails to detect an association between HIV infection and probability of being a TB infection source. We conclude that application of a pipeline, where one first uses TransPhylo and sparsely sampled surveillance data to infer transmission events and then estimates effects of host characteristics on probabilities of these events, should be accompanied by a realistic simulation study to better understand biases stemming from imprecise transmission event inference.

Stem cell maintenance in multilayered shoot apical meristems (SAMs) of plants requires strict regulation of cell growth and division. Exactly how the complex milieu of chemical and mechanical signals interact in the central region of the SAM to regulate cell division plane orientation is not well understood. In this paper, simulations using a newly developed multiscale computational model are combined with experimental studies to suggest and test three hypothesized mechanisms for the regulation of cell division plane orientation and the direction of anisotropic cell expansion in the corpus. Simulations predict that in the Apical corpus, WUSCHEL and cytokinin regulate the direction of anisotropic cell expansion, and cells divide according to tensile stress on the cell wall. In the Basal corpus, model simulations suggest dual roles for WUSCHEL and cytokinin in regulating both the direction of anisotropic cell expansion and cell division plane orientation. Simulation results are followed by a detailed analysis of changes in cell characteristics upon manipulation of WUSCHEL and cytokinin in experiments that support model predictions. Moreover, simulations predict that this layer-specific mechanism maintains both the experimentally observed shape and structure of the SAM as well as the distribution of WUSCHEL in the tissue. This provides an additional link between the roles of WUSCHEL, cytokinin, and mechanical stress in regulating SAM growth and proper stem cell maintenance in the SAM.

The sinoatrial node (SAN), the primary pacemaker of the heart, is responsible for the initiation and robust regulation of sinus rhythm. 3D mapping studies of the ex-vivo human heart suggested that the robust regulation of sinus rhythm relies on specialized fibrotically-insulated pacemaker compartments (head, center and tail) with heterogeneous expressions of key ion channels and receptors. They also revealed up to five sinoatrial conduction pathways (SACPs), which electrically connect the SAN with neighboring right atrium (RA). To elucidate the role of these structural-molecular factors in the functional robustness of human SAN, we developed comprehensive biophysical computer models of the SAN based on 3D structural, functional and molecular mapping of ex-vivo human hearts. Our key finding is that the electrical insulation of the SAN except SACPs, the heterogeneous expression of If, INa currents and adenosine A1 receptors (A1R) across SAN pacemaker-conduction compartments are required to experimentally reproduce observed SAN activation patterns and important phenomena such as shifts of the leading pacemaker and preferential SACP. In particular, we found that the insulating border between the SAN and RA, is required for robust SAN function and protection from SAN arrest during adenosine challenge. The heterogeneity in the expression of A1R within the human SAN compartments underlies the direction of pacemaker shift and preferential SACPs in the presence of adenosine. Alterations of INa current and fibrotic remodelling in SACPs can significantly modulate SAN conduction and shift the preferential SACP/exit from SAN. Finally, we show that disease-induced fibrotic remodeling, INa suppression or increased adenosine make the human SAN vulnerable to pacing-induced exit blocks and reentrant arrhythmia. In summary, our computer model recapitulates the structural and functional features of the human SAN and can be a valuable tool for investigating mechanisms of SAN automaticity and conduction as well as SAN arrhythmia mechanisms under different pathophysiological conditions. Author summary: The human heart is driven and modulated by the sinoatrial node (SAN), our body's natural pacemaker. Recent studies using explanted human hearts discovered that to regulate heart rhythm robustly, the SAN has three pacemaker compartments–SAN head, center and tail characterized by heterogeneous expression of key ion channels and receptors. In addition, the fibrotically-insulated SAN electrically connects with the right atrium (RA) through up to five sinoatrial conduction pathways (SACPs). Due to the complexities of the human 3D structure and limited functional data on SAN conduction, the specific role of the SAN insulation/border, distinct SACPs and intranodal pacemaker molecular heterogeneity in regulating sinus rhythm in health and diseased hearts remain debatable. The goals of this study were to define the key factors influencing human SAN pacemaking function and SAN dysfunction by developing and utilizing computer models of the human SAN. This study presents the first comprehensive biophysical computer model of the human SAN complex based on direct molecular, structural and functional studies in the ex-vivo human heart. Our data show that the computer models can closely replicate pacemaking, SAN activation patterns and exit sites/earliest atrial activation through preferential SACPs, as well as physiological changes including the shift of the leading pacemaker in the presence of adenosine reported in the human heart ex-vivo. More importantly, the novel computer modeling simulation results illustrate the crucial role of the structural and electrical heterogeneity of the human SAN in pacemaking and conduction. Our human-specific SAN computer model represents a valuable tool for investigating SAN automaticity, conduction and arrhythmia mechanisms under normal and various disease conditions. [ABSTRACT FROM AUTHOR]

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

Defects in axonal transport may partly underpin the differences between the observed pathophysiology of Alzheimer's disease (AD) and that of other non-amyloidogenic tauopathies. Particularly, pathological tau variants may have molecular properties that dysregulate motor proteins responsible for the anterograde-directed transport of tau in a disease-specific fashion. Here we develop the first computational model of tau-modified axonal transport that produces directional biases in the spread of tau pathology. We simulated the spatiotemporal profiles of soluble and insoluble tau species in a multicompartment, two-neuron system using biologically plausible parameters and time scales. Changes in the balance of tau transport feedback parameters can elicit anterograde and retrograde biases in the distributions of soluble and insoluble tau between compartments in the system. Aggregation and fragmentation parameters can also perturb this balance, suggesting a complex interplay between these distinct molecular processes. Critically, we show that the model faithfully recreates the characteristic network spread biases in both AD-like and non-AD-like mouse tauopathy models. Tau transport feedback may therefore help link microscopic differences in tau conformational states and the resulting variety in clinical presentations.

The effectiveness of immune responses depends on the precision of stimulus-responsive gene expression programs. Cells specify which genes to express by activating stimulus-specific combinations of stimulus-induced transcription factors (TFs). Their activities are decoded by a gene regulatory strategy (GRS) associated with each response gene. Here, we examined whether the GRSs of target genes may be inferred from stimulus-response (input-output) datasets, which remains an unresolved model-identifiability challenge. We developed a mechanistic modeling framework and computational workflow to determine the identifiability of all possible combinations of synergistic (AND) or non-synergistic (OR) GRSs involving three transcription factors. Considering different sets of perturbations for stimulus-response studies, we found that two thirds of GRSs are easily distinguishable but that substantially more quantitative data is required to distinguish the remaining third. To enhance the accuracy of the inference with timecourse experimental data, we developed an advanced error model that avoids error overestimates by distinguishing between value and temporal error. Incorporating this error model into a Bayesian framework, we show that GRS models can be identified for individual genes by considering multiple datasets. Our analysis rationalizes the allocation of experimental resources by identifying most informative TF stimulation conditions. Applying this computational workflow to experimental data of immune response genes in macrophages, we found that a much greater fraction of genes are combinatorially controlled than previously reported by considering compensation among transcription factors. Specifically, we revealed that a group of known NFκB target genes may also be regulated by IRF3, which is supported by chromatin immuno-precipitation analysis. Our study provides a computational workflow for designing and interpreting stimulus-response gene expression studies to identify underlying gene regulatory strategies and further a mechanistic understanding.

Genetic predisposition for complex traits often acts through multiple tissues at different time points during development. As a simple example, the genetic predisposition for obesity could be manifested either through inherited variants that control metabolism through regulation of genes expressed in the brain, or that control fat storage through dysregulation of genes expressed in adipose tissue, or both. Here we describe a statistical approach that leverages tissue-specific expression quantitative trait loci (eQTLs) corresponding to tissue-specific genes to prioritize a relevant tissue underlying the genetic predisposition of a given individual for a complex trait. Unlike existing approaches that prioritize relevant tissues for the trait in the population, our approach probabilistically quantifies the tissue-wise genetic contribution to the trait for a given individual. We hypothesize that for a subgroup of individuals the genetic contribution to the trait can be mediated primarily through a specific tissue. Through simulations using the UK Biobank, we show that our approach can predict the relevant tissue accurately and can cluster individuals according to their tissue-specific genetic architecture. We analyze body mass index (BMI) and waist to hip ratio adjusted for BMI (WHRadjBMI) in the UK Biobank to identify subgroups of individuals whose genetic predisposition act primarily through brain versus adipose tissue, and adipose versus muscle tissue, respectively. Notably, we find that these individuals have specific phenotypic features beyond BMI and WHRadjBMI that distinguish them from random individuals in the data, suggesting biological effects of tissue-specific genetic contribution for these traits.

State-dependent sodium channel blockers are often prescribed to treat cardiac arrhythmias, but many sodium channel blockers are known to have pro-arrhythmic side effects. While the anti and proarrhythmic potential of a sodium channel blocker is thought to depend on the characteristics of its rate-dependent block, the mechanisms linking these two attributes are unclear. Furthermore, how specific properties of rate-dependent block arise from the binding kinetics of a particular drug is poorly understood. Here, we examine the rate-dependent effects of the sodium channel blocker lidocaine by constructing and analyzing a novel drug-channel interaction model. First, we identify the predominant mode of lidocaine binding in a 24 variable Markov model for lidocaine-sodium channel interaction by Moreno et al. Specifically, we find that (1) the vast majority of lidocaine bound to sodium channels is in the neutral form, i.e., the binding of charged lidocaine to sodium channels is negligible, and (2) neutral lidocaine binds almost exclusively to inactivated channels and, upon binding, immobilizes channels in the inactivated state. We then develop a novel 3-variable lidocaine-sodium channel interaction model that incorporates only the predominant mode of drug binding. Our low-dimensional model replicates an extensive amount of the voltage-clamp data used to parameterize the Moreno et al. model. Furthermore, the effects of lidocaine on action potential upstroke velocity and conduction velocity in our model are similar to those predicted by the Moreno et al. model. By exploiting the low-dimensionality of our model, we derive an algebraic expression for level of rate-dependent block as a function of pacing frequency, restitution properties, diastolic and plateau potentials, and drug binding rate constants. Our model predicts that the level of rate-dependent block is sensitive to alterations in restitution properties and increases in diastolic potential, but it is insensitive to variations in the shape of the action potential waveform and lidocaine binding rates.

Many cellular processes require cell polarization to be maintained as the cell changes shape, grows or moves. Without feedback mechanisms relaying information about cell shape to the polarity molecular machinery, the coordination between cell polarization and morphogenesis, movement or growth would not be possible. Here we theoretically and computationally study the role of a genetically-encoded mechanical feedback (in the Cell Wall Integrity pathway) as a potential coordination mechanism between cell morphogenesis and polarity during budding yeast mating projection growth. We developed a coarse-grained continuum description of the coupled dynamics of cell polarization and morphogenesis as well as 3D stochastic simulations of the molecular polarization machinery in the evolving cell shape. Both theoretical approaches show that in the absence of mechanical feedback (or in the presence of weak feedback), cell polarity cannot be maintained at the projection tip during growth, with the polarization cap wandering off the projection tip, arresting morphogenesis. In contrast, for mechanical feedback strengths above a threshold, cells can robustly maintain cell polarization at the tip and simultaneously sustain mating projection growth. These results indicate that the mechanical feedback encoded in the Cell Wall Integrity pathway can provide important positional information to the molecular machinery in the cell, thereby enabling the coordination of cell polarization and morphogenesis.

Mathematical models of metabolic networks utilize simulation to study system-level mechanisms and functions. Various approaches have been used to model the steady state behavior of metabolic networks using genome-scale reconstructions, but formulating dynamic models from such reconstructions continues to be a key challenge. Here, we present the Mass Action Stoichiometric Simulation Python (MASSpy) package, an open-source computational framework for dynamic modeling of metabolism. MASSpy utilizes mass action kinetics and detailed chemical mechanisms to build dynamic models of complex biological processes. MASSpy adds dynamic modeling tools to the COnstraint-Based Reconstruction and Analysis Python (COBRApy) package to provide an unified framework for constraint-based and kinetic modeling of metabolic networks. MASSpy supports high-performance dynamic simulation through its implementation of libRoadRunner: the Systems Biology Markup Language (SBML) simulation engine. Three examples are provided to demonstrate how to use MASSpy: (1) a validation of the MASSpy modeling tool through dynamic simulation of detailed mechanisms of enzyme regulation; (2) a feature demonstration using a workflow for generating ensemble of kinetic models using Monte Carlo sampling to approximate missing numerical values of parameters and to quantify biological uncertainty, and (3) a case study in which MASSpy is utilized to overcome issues that arise when integrating experimental data with the computation of functional states of detailed biological mechanisms. MASSpy represents a powerful tool to address challenges that arise in dynamic modeling of metabolic networks, both at small and large scales.

Coronary artery thrombosis is the major risk associated with Kawasaki disease (KD). Long-term management of KD patients with persistent aneurysms requires a thrombotic risk assessment and clinical decisions regarding the administration of anticoagulation therapy. Computational fluid dynamics has demonstrated that abnormal KD coronary artery hemodynamics can be associated with thrombosis. However, the underlying mechanisms of clot formation are not yet fully understood. Here we present a new model incorporating data from patient-specific simulated velocity fields to track platelet activation and accumulation. We use a system of Reaction-Advection-Diffusion equations solved with a stabilized finite element method to describe the evolution of non-activated platelets and activated platelet concentrations [AP], local concentrations of adenosine diphosphate (ADP) and poly-phosphate (PolyP). The activation of platelets is modeled as a function of shear-rate exposure and local concentration of agonists. We compared the distribution of activated platelets in a healthy coronary case and six cases with coronary artery aneurysms caused by KD, including three with confirmed thrombosis. Results show spatial correlation between regions of higher concentration of activated platelets and the reported location of the clot, suggesting predictive capabilities of this model towards identifying regions at high risk for thrombosis. Also, the concentration levels of ADP and PolyP in cases with confirmed thrombosis are higher than the reported critical values associated with platelet aggregation (ADP) and activation of the intrinsic coagulation pathway (PolyP). These findings suggest the potential initiation of a coagulation pathway even in the absence of an extrinsic factor. Finally, computational simulations show that in regions of flow stagnation, biochemical activation, as a result of local agonist concentration, is dominant. Identifying the leading factors to a pro-coagulant environment in each case-mechanical or biochemical-could help define improved strategies for thrombosis prevention tailored for each patient.

A key question in SARS-CoV-2 infection is why viral loads and patient outcomes vary dramatically across individuals. Because spatial-temporal dynamics of viral spread and immune response are challenging to study in vivo, we developed Spatial Immune Model of Coronavirus (SIMCoV), a scalable computational model that simulates hundreds of millions of lung cells, including respiratory epithelial cells and T cells. SIMCoV replicates viral growth dynamics observed in patients and shows how spatially dispersed infections can lead to increased viral loads. The model also shows how the timing and strength of the T cell response can affect viral persistence, oscillations, and control. By incorporating spatial interactions, SIMCoV provides a parsimonious explanation for the dramatically different viral load trajectories among patients by varying only the number of initial sites of infection and the magnitude and timing of the T cell immune response. When the branching airway structure of the lung is explicitly represented, we find that virus spreads faster than in a 2D layer of epithelial cells, but much more slowly than in an undifferentiated 3D grid or in a well-mixed differential equation model. These results illustrate how realistic, spatially explicit computational models can improve understanding of within-host dynamics of SARS-CoV-2 infection.

The broad concept of emergence is instrumental in various of the most challenging open scientific questions-yet, few quantitative theories of what constitutes emergent phenomena have been proposed. This article introduces a formal theory of causal emergence in multivariate systems, which studies the relationship between the dynamics of parts of a system and macroscopic features of interest. Our theory provides a quantitative definition of downward causation, and introduces a complementary modality of emergent behaviour-which we refer to as causal decoupling. Moreover, the theory allows practical criteria that can be efficiently calculated in large systems, making our framework applicable in a range of scenarios of practical interest. We illustrate our findings in a number of case studies, including Conway's Game of Life, Reynolds' flocking model, and neural activity as measured by electrocorticography.

Mosquitoes are important vectors for pathogens that infect humans and other vertebrate animals. Some aspects of adult mosquito behavior and mosquito ecology play an important role in determining the capacity of vector populations to transmit pathogens. Here, we re-examine factors affecting the transmission of pathogens by mosquitoes using a new approach. Unlike most previous models, this framework considers the behavioral states and state transitions of adult mosquitoes through a sequence of activity bouts. We developed a new framework for individual-based simulation models called MBITES (Mosquito Bout-based and Individual-based Transmission Ecology Simulator). In MBITES, it is possible to build models that simulate the behavior and ecology of adult mosquitoes in exquisite detail on complex resource landscapes generated by spatial point processes. We also developed an ordinary differential equation model which is the Kolmogorov forward equations for models developed in MBITES under a specific set of simplifying assumptions. While mosquito infection and pathogen development are one possible part of a mosquito's state, that is not our main focus. Using extensive simulation using some models developed in MBITES, we show that vectorial capacity can be understood as an emergent property of simple behavioral algorithms interacting with complex resource landscapes, and that relative density or sparsity of resources and the need to search can have profound consequences for mosquito populations' capacity to transmit pathogens.

Recent advances in electron microscopy have enabled the imaging of single cells in 3D at nanometer length scale resolutions. An uncharted frontier for in silico biology is the ability to simulate cellular processes using these observed geometries. Enabling such simulations requires watertight meshing of electron micrograph images into 3D volume meshes, which can then form the basis of computer simulations of such processes using numerical techniques such as the finite element method. In this paper, we describe the use of our recently rewritten mesh processing software, GAMer 2, to bridge the gap between poorly conditioned meshes generated from segmented micrographs and boundary marked tetrahedral meshes which are compatible with simulation. We demonstrate the application of a workflow using GAMer 2 to a series of electron micrographs of neuronal dendrite morphology explored at three different length scales and show that the resulting meshes are suitable for finite element simulations. This work is an important step towards making physical simulations of biological processes in realistic geometries routine. Innovations in algorithms to reconstruct and simulate cellular length scale phenomena based on emerging structural data will enable realistic physical models and advance discovery at the interface of geometry and cellular processes. We posit that a new frontier at the intersection of computational technologies and single cell biology is now open.

Glioblastomas are aggressive primary brain tumors known for their inter- and intratumor heterogeneity. This disease is uniformly fatal, with intratumor heterogeneity the major reason for treatment failure and recurrence. Just like the nature vs nurture debate, heterogeneity can arise from intrinsic or environmental influences. Whilst it is impossible to clinically separate observed behavior of cells from their environmental context, using a mathematical framework combined with multiscale data gives us insight into the relative roles of variation from different sources. To better understand the implications of intratumor heterogeneity on therapeutic outcomes, we created a hybrid agent-based mathematical model that captures both the overall tumor kinetics and the individual cellular behavior. We track single cells as agents, cell density on a coarser scale, and growth factor diffusion and dynamics on a finer scale over time and space. Our model parameters were fit utilizing serial MRI imaging and cell tracking data from ex vivo tissue slices acquired from a growth-factor driven glioblastoma murine model. When fitting our model to serial imaging only, there was a spectrum of equally-good parameter fits corresponding to a wide range of phenotypic behaviors. When fitting our model using imaging and cell scale data, we determined that environmental heterogeneity alone is insufficient to match the single cell data, and intrinsic heterogeneity is required to fully capture the migration behavior. The wide spectrum of in silico tumors also had a wide variety of responses to an application of an anti-proliferative treatment. Recurrent tumors were generally less proliferative than pre-treatment tumors as measured via the model simulations and validated from human GBM patient histology. Further, we found that all tumors continued to grow with an anti-migratory treatment alone, but the anti-proliferative/anti-migratory combination generally showed improvement over an anti-proliferative treatment alone. Together our results emphasize the need to better understand the underlying phenotypes and tumor heterogeneity present in a tumor when designing therapeutic regimens.

Sequential sampling models such as the drift diffusion model (DDM) have a long tradition in research on perceptual decision-making, but mounting evidence suggests that these models can account for response time (RT) distributions that arise during reinforcement learning and value-based decision-making. Building on this previous work, we implemented the DDM as the choice rule in inter-temporal choice (temporal discounting) and risky choice (probability discounting) using hierarchical Bayesian parameter estimation. We validated our approach in data from nine patients with focal lesions to the ventromedial prefrontal cortex / medial orbitofrontal cortex (vmPFC/mOFC) and nineteen age- and education-matched controls. Model comparison revealed that, for both tasks, the data were best accounted for by a variant of the drift diffusion model including a non-linear mapping from value-differences to trial-wise drift rates. Posterior predictive checks confirmed that this model provided a superior account of the relationship between value and RT. We then applied this modeling framework and 1) reproduced our previous results regarding temporal discounting in vmPFC/mOFC patients and 2) showed in a previously unpublished data set on risky choice that vmPFC/mOFC patients exhibit increased risk-taking relative to controls. Analyses of DDM parameters revealed that patients showed substantially increased non-decision times and reduced response caution during risky choice. In contrast, vmPFC/mOFC damage abolished neither scaling nor asymptote of the drift rate. Relatively intact value processing was also confirmed using DDM mixture models, which revealed that in both groups >98% of trials were better accounted for by a DDM with value modulation than by a null model without value modulation. Our results highlight that novel insights can be gained from applying sequential sampling models in studies of inter-temporal and risky decision-making in cognitive neuroscience.

Calmodulin-dependent kinase II (CaMKII) has long been known to play an important role in learning and memory as well as long term potentiation (LTP). More recently it has been suggested that it might be involved in the time averaging of synaptic signals, which can then lead to the high precision of information stored at a single synapse. However, the role of the scaffolding molecule, neurogranin (Ng), in governing the dynamics of CaMKII is not yet fully understood. In this work, we adopt a rule-based modeling approach through the Monte Carlo method to study the effect of Ca2+ signals on the dynamics of CaMKII phosphorylation in the postsynaptic density (PSD). Calcium surges are observed in synaptic spines during an EPSP and back-propagating action potential due to the opening of NMDA receptors and voltage dependent calcium channels. Using agent-based models, we computationally investigate the dynamics of phosphorylation of CaMKII monomers and dodecameric holoenzymes. The scaffolding molecule, Ng, when present in significant concentration, limits the availability of free calmodulin (CaM), the protein which activates CaMKII in the presence of calcium. We show that Ng plays an important modulatory role in CaMKII phosphorylation following a surge of high calcium concentration. We find a non-intuitive dependence of this effect on CaM concentration that results from the different affinities of CaM for CaMKII depending on the number of calcium ions bound to the former. It has been shown previously that in the absence of phosphatase, CaMKII monomers integrate over Ca2+ signals of certain frequencies through autophosphorylation (Pepke et al, Plos Comp. Bio., 2010). We also study the effect of multiple calcium spikes on CaMKII holoenzyme autophosphorylation, and show that in the presence of phosphatase, CaMKII behaves as a leaky integrator of calcium signals, a result that has been recently observed in vivo. Our models predict that the parameters of this leaky integrator are finely tuned through the interactions of Ng, CaM, CaMKII, and PP1, providing a mechanism to precisely control the sensitivity of synapses to calcium signals. Author Summary not valid for PLOS ONE submissions.

Interpretations of elevated blood levels of alanine aminotransferase (ALT) for drug-induced liver injury often assume that the biomarker is released passively from dying cells. However, the mechanisms driving that release have not been explored experimentally. The usefulness of ALT and related biomarkers will improve by developing mechanism-based explanations of elevated levels that can be expanded and elaborated incrementally. We provide the means to challenge the ability of closely related model mechanisms to generate patterns of simulated hepatic injury and ALT release that scale (or not) to be quantitatively similar to the wet-lab validation targets, which are elevated plasma ALT values following acetaminophen (APAP) exposure in mice. We build on a published model mechanism that helps explain the generation of characteristic spatiotemporal features of APAP hepatotoxicity within hepatic lobules. Discrete event and agent-oriented software methods are most prominent. We instantiate and leverage a small constellation of concrete model mechanisms. Their details during execution help bring into focus ways in which particular sources of uncertainty become entangled with cause-effect details within and across several levels. We scale ALT amounts in virtual mice directly to target plasma ALT values in individual mice. A virtual experiment comprises a set of Monte Carlo simulations. We challenge the sufficiency of four potentially explanatory theories for ALT release. The first of the tested model theories failed to achieve the initial validation target, but each of the three others succeeded. Results for one of the three model mechanisms matched all target ALT values quantitatively. It explains how ALT externalization is the combined consequence of lobular-location-dependent drug-induced cellular damage and hepatocyte death. Falsification of one (or more) of the model mechanisms provides new knowledge and incrementally shrinks the constellation of model mechanisms. The modularity and biomimicry of our explanatory models enable seamless transition from mice to humans.

Motivated by growing evidence for pathway heterogeneity and alternative functions of molecular machines, we demonstrate a computational approach for investigating two questions: (1) Are there multiple mechanisms (state-space pathways) by which a machine can perform a given function, such as cotransport across a membrane? (2) How can additional functionality, such as proofreading/error-correction, be built into machine function using standard biochemical processes? Answers to these questions will aid both the understanding of molecular-scale cell biology and the design of synthetic machines. Focusing on transport in this initial study, we sample a variety of mechanisms by employing Metropolis Markov chain Monte Carlo. Trial moves adjust transition rates among an automatically generated set of conformational and binding states while maintaining fidelity to thermodynamic principles and a user-supplied fitness/functionality goal. Each accepted move generates a new model. The simulations yield both single and mixed reaction pathways for cotransport in a simple environment with a single substrate along with a driving ion. In a "competitive" environment including an additional decoy substrate, several qualitatively distinct reaction pathways are found which are capable of extremely high discrimination coupled to a leak of the driving ion, akin to proofreading. The array of functional models would be difficult to find by intuition alone in the complex state-spaces of interest.

Recent years have seen rising incidence of dengue and large outbreaks of Zika and chikungunya, which are all caused by viruses transmitted by Aedes aegypti mosquitoes. In most settings, the primary intervention against Aedes-transmitted viruses is vector control, such as indoor, ultra-low volume (ULV) spraying. Targeted indoor residual spraying (TIRS) has the potential to more effectively impact Aedes-borne diseases, but its implementation requires careful planning and evaluation. The optimal time to deploy these interventions and their relative epidemiological effects are, however, not well understood. We used an agent-based model of dengue virus transmission calibrated to data from Iquitos, Peru to assess the epidemiological effects of these interventions under differing strategies for deploying them. Specifically, we compared strategies where spray application was initiated when incidence rose above a threshold based on incidence in recent years to strategies where spraying occurred at the same time(s) each year. In the absence of spraying, the model predicted 361,000 infections [inter-quartile range (IQR): 347,000-383,000] in the period 2000-2010. The ULV strategy with the fewest median infections was spraying twice yearly, in March and October, which led to a median of 172,000 infections [IQR: 158,000-183,000], a 52% reduction from baseline. Compared to spraying once yearly in September, the best threshold-based strategy utilizing ULV had fewer median infections (254,000 vs. 261,000), but required more spraying (351 vs. 274 days). For TIRS, the best strategy was threshold-based, which led to the fewest infections of all strategies tested (9,900; [IQR: 8,720-11,400], a 94% reduction), and required fewer days spraying than the equivalent ULV strategy (280). Although spraying twice each year is likely to avert the most infections, our results indicate that a threshold-based strategy can become an alternative to better balance the translation of spraying effort into impact, particularly if used with a residual insecticide.

Infectious disease surveillance systems provide vital data for guiding disease prevention and control policies, yet the formalization of methods to optimize surveillance networks has largely been overlooked. Decisions surrounding surveillance design parameters-such as the number and placement of surveillance sites, target populations, and case definitions-are often determined by expert opinion or deference to operational considerations, without formal analysis of the influence of design parameters on surveillance objectives. Here we propose a simulation framework to guide evidence-based surveillance network design to better achieve specific surveillance goals with limited resources. We define evidence-based surveillance design as an optimization problem, acknowledging the many operational constraints under which surveillance systems operate, the many dimensions of surveillance system design, the multiple and competing goals of surveillance, and the complex and dynamic nature of disease systems. We describe an analytical framework-the Disease Surveillance Informatics Optimization and Simulation (DIOS) framework-for the identification of optimal surveillance designs through mathematical representations of disease and surveillance processes, definition of objective functions, and numerical optimization. We then apply the framework to the problem of selecting candidate sites to expand an existing surveillance network under alternative objectives of: (1) improving spatial prediction of disease prevalence at unmonitored sites; or (2) estimating the observed effect of a risk factor on disease. Results of this demonstration illustrate how optimal designs are sensitive to both surveillance goals and the underlying spatial pattern of the target disease. The findings affirm the value of designing surveillance systems through quantitative and adaptive analysis of network characteristics and performance. The framework can be applied to the design of surveillance systems tailored to setting-specific disease transmission dynamics and surveillance needs, and can yield improved understanding of tradeoffs between network architectures.

Top-down proteomics has enabled the elucidation of heterogeneous protein complexes with different cofactors, post-translational modifications, and protein membership. This heterogeneity is believed to play a previously unknown role in cellular processes. The different molecular forms of a protein complex have come to be called "complex isoform" or "complexoform". Despite the elucidation of the complexoform, it remains unclear how and whether cellular circuits control the distribution of a complexoform. To help address this issue, we first simulate a generic three-protein complexoform to reveal the control of its distribution by the timing of gene transcription, mRNA translation, and protein transport. Overall, we ran 265 computational experiments: each averaged over 1,000 stochastic simulations. Based on the experiments, we show that genes arranged in a single operon, a cascade, or as two operons all give rise to the different protein composition of complexoform because of timing differences in protein-synthesis order. We also show that changes in the kinetics of expression, protein transport, or protein binding dramatically alter the distribution of the complexoform. Furthermore, both stochastic and transient kinetics control the assembly of the complexoform when the expression and assembly occur concurrently. We test our model against the biological cellulosome system. With biologically relevant rates, we find that the genetic circuitry controls the average final complexoform assembly and the variation in the assembly structure. Our results highlight the importance of both the genetic circuit architecture and kinetics in determining the distribution of a complexoform. Our work has a broad impact on our understanding of non-equilibrium processes in both living and synthetic biological systems.

For a chemical signal to propagate across a cell, it must navigate a tortuous environment involving a variety of organelle barriers. In this work we study mathematical models for a basic chemical signal, the arrival times at the nuclear membrane of proteins that are activated at the cell membrane and diffuse throughout the cytosol. Organelle surfaces within human B cells are reconstructed from soft X-ray tomographic images, and modeled as reflecting barriers to the molecules' diffusion. We show that signal inactivation sharpens signals, reducing variability in the arrival time at the nuclear membrane. Inactivation can also compensate for an observed slowdown in signal propagation induced by the presence of organelle barriers, leading to arrival times at the nuclear membrane that are comparable to models in which the cytosol is treated as an open, empty region. In the limit of strong signal inactivation this is achieved by filtering out molecules that traverse non-geodesic paths.

Birth-death processes have given biologists a model-based framework to answer questions about changes in the birth and death rates of lineages in a phylogenetic tree. Therefore birth-death models are central to macroevolutionary as well as phylodynamic analyses. Early approaches to studying temporal variation in birth and death rates using birth-death models faced difficulties due to the restrictive choices of birth and death rate curves through time. Sufficiently flexible time-varying birth-death models are still lacking. We use a piecewise-constant birth-death model, combined with both Gaussian Markov random field (GMRF) and horseshoe Markov random field (HSMRF) prior distributions, to approximate arbitrary changes in birth rate through time. We implement these models in the widely used statistical phylogenetic software platform RevBayes, allowing us to jointly estimate birth-death process parameters, phylogeny, and nuisance parameters in a Bayesian framework. We test both GMRF-based and HSMRF-based models on a variety of simulated diversification scenarios, and then apply them to both a macroevolutionary and an epidemiological dataset. We find that both models are capable of inferring variable birth rates and correctly rejecting variable models in favor of effectively constant models. In general the HSMRF-based model has higher precision than its GMRF counterpart, with little to no loss of accuracy. Applied to a macroevolutionary dataset of the Australian gecko family Pygopodidae (where birth rates are interpretable as speciation rates), the GMRF-based model detects a slow decrease whereas the HSMRF-based model detects a rapid speciation-rate decrease in the last 12 million years. Applied to an infectious disease phylodynamic dataset of sequences from HIV subtype A in Russia and Ukraine (where birth rates are interpretable as the rate of accumulation of new infections), our models detect a strongly elevated rate of infection in the 1990s.

The lumenal pH of an organelle is one of its defining characteristics and central to its biological function. Experiments have elucidated many of the key pH regulatory elements and how they vary from compartment-to-compartment, and continuum mathematical models have played an important role in understanding how these elements (proton pumps, counter-ion fluxes, membrane potential, buffering capacity, etc.) work together to achieve specific pH setpoints. While continuum models have proven successful in describing ion regulation at the cellular length scale, it is unknown if they are valid at the subcellular level where volumes are small, ion numbers may fluctuate wildly, and biochemical heterogeneity is large. Here, we create a discrete, stochastic (DS) model of vesicular acidification to answer this question. We used this simplified model to analyze pH measurements of isolated vesicles containing single proton pumps and compared these results to solutions from a continuum, ordinary differential equations (ODE)-based model. Both models predict similar parameter estimates for the mean proton pumping rate, membrane permeability, etc., but, as expected, the ODE model fails to report on the fluctuations in the system. The stochastic model predicts that pH fluctuations decrease during acidification, but noise analysis of single-vesicle data confirms our finding that the experimental noise is dominated by the fluorescent dye, and it reveals no insight into the true noise in the proton fluctuations. Finally, we again use the reduced DS model explore the acidification of large, lysosome-like vesicles to determine how stochastic elements, such as variations in proton-pump copy number and cycling between on and off states, impact the pH setpoint and fluctuations around this setpoint.

Response to acid stress is critical for Escherichia coli to successfully complete its life-cycle by passing through the stomach to colonize the digestive tract. To develop a fundamental understanding of this response, we established a molecular mechanistic description of acid stress mitigation responses in E. coli and integrated them with a genome-scale model of its metabolism and macromolecular expression (ME-model). We considered three known mechanisms of acid stress mitigation: 1) change in membrane lipid fatty acid composition, 2) change in periplasmic protein stability over external pH and periplasmic chaperone protection mechanisms, and 3) change in the activities of membrane proteins. After integrating these mechanisms into an established ME-model, we could simulate their responses in the context of other cellular processes. We validated these simulations using RNA sequencing data obtained from five E. coli strains grown under external pH ranging from 5.5 to 7.0. We found: i) that for the differentially expressed genes accounted for in the ME-model, 80% of the upregulated genes were correctly predicted by the ME-model, and ii) that these genes are mainly involved in translation processes (45% of genes), membrane proteins and related processes (18% of genes), amino acid metabolism (12% of genes), and cofactor and prosthetic group biosynthesis (8% of genes). We also demonstrated several intervention strategies on acid tolerance that can be simulated by the ME-model. We thus established a quantitative framework that describes, on a genome-scale, the acid stress mitigation response of E. coli that has both scientific and practical uses.

Bacterial populations that colonize a host can play important roles in host health, including serving as a reservoir that transmits to other hosts and from which invasive strains emerge, thus emphasizing the importance of understanding rates of acquisition and clearance of colonizing populations. Studies of colonization dynamics have been based on assessment of whether serial samples represent a single population or distinct colonization events. With the use of whole genome sequencing to determine genetic distance between isolates, a common solution to estimate acquisition and clearance rates has been to assume a fixed genetic distance threshold below which isolates are considered to represent the same strain. However, this approach is often inadequate to account for the diversity of the underlying within-host evolving population, the time intervals between consecutive measurements, and the uncertainty in the estimated acquisition and clearance rates. Here, we present a fully Bayesian model that provides probabilities of whether two strains should be considered the same, allowing us to determine bacterial clearance and acquisition from genomes sampled over time. Our method explicitly models the within-host variation using population genetic simulation, and the inference is done using a combination of Approximate Bayesian Computation (ABC) and Markov Chain Monte Carlo (MCMC). We validate the method with multiple carefully conducted simulations and demonstrate its use in practice by analyzing a collection of methicillin resistant Staphylococcus aureus (MRSA) isolates from a large recently completed longitudinal clinical study. An R-code implementation of the method is freely available at: https://github.com/mjarvenpaa/bacterial-colonization-model.

Quantitative viral outgrowth assays (QVOA) use limiting dilutions of CD4+ T cells to measure the size of the latent HIV-1 reservoir, a major obstacle to curing HIV-1. Efforts to reduce the reservoir require assays that can reliably quantify its size in blood and tissues. Although QVOA is regarded as a "gold standard" for reservoir measurement, little is known about its accuracy and precision or about how cell storage conditions or laboratory-specific practices affect results. Owing to this lack of knowledge, confidence intervals around reservoir size estimates-as well as judgments of the ability of therapeutic interventions to alter the size of the replication-competent but transcriptionally inactive latent reservoir-rely on theoretical statistical assumptions about dilution assays. To address this gap, we have carried out a Bayesian statistical analysis of QVOA reliability on 75 split samples of peripheral blood mononuclear cells (PBMC) from 5 antiretroviral therapy (ART)-suppressed participants, measured using four different QVOAs at separate labs, estimating assay precision and the effect of frozen cell storage on estimated reservoir size. We found that typical assay results are expected to differ from the true value by a factor of 1.6 to 1.9 up or down. Systematic assay differences comprised a 24-fold range between the assays with highest and lowest scales, likely reflecting differences in viral outgrowth readout and input cell stimulation protocols. We also found that controlled-rate freezing and storage of samples did not cause substantial differences in QVOA compared to use of fresh cells (95% probability of < 2-fold change), supporting continued use of frozen storage to allow transport and batched analysis of samples. Finally, we simulated an early-phase clinical trial to demonstrate that batched analysis of pre- and post-therapy samples may increase power to detect a three-fold reservoir reduction by 15 to 24 percentage points.

Elaboration of Bayesian phylogenetic inference methods has continued at pace in recent years with major new advances in nearly all aspects of the joint modelling of evolutionary data. It is increasingly appreciated that some evolutionary questions can only be adequately answered by combining evidence from multiple independent sources of data, including genome sequences, sampling dates, phenotypic data, radiocarbon dates, fossil occurrences, and biogeographic range information among others. Including all relevant data into a single joint model is very challenging both conceptually and computationally. Advanced computational software packages that allow robust development of compatible (sub-)models which can be composed into a full model hierarchy have played a key role in these developments. Developing such software frameworks is increasingly a major scientific activity in its own right, and comes with specific challenges, from practical software design, development and engineering challenges to statistical and conceptual modelling challenges. BEAST 2 is one such computational software platform, and was first announced over 4 years ago. Here we describe a series of major new developments in the BEAST 2 core platform and model hierarchy that have occurred since the first release of the software, culminating in the recent 2.5 release.

Prophylactic vaccination is a powerful tool for reducing the burden of infectious diseases, due to a combination of direct protection of vaccinees and indirect protection of others via herd immunity. Computational models play an important role in devising strategies for vaccination by making projections of its impacts on public health. Such projections are subject to uncertainty about numerous factors, however. For example, many vaccine efficacy trials focus on measuring protection against disease rather than protection against infection, leaving the extent of breakthrough infections (i.e., disease ameliorated but infection unimpeded) among vaccinees unknown. Our goal in this study was to quantify the extent to which uncertainty about breakthrough infections results in uncertainty about vaccination impact, with a focus on vaccines for dengue. To realistically account for the many forms of heterogeneity in dengue virus (DENV) transmission, which could have implications for the dynamics of indirect protection, we used a stochastic, agent-based model for DENV transmission informed by more than a decade of empirical studies in the city of Iquitos, Peru. Following 20 years of routine vaccination of nine-year-old children at 80% coverage, projections of the proportion of disease episodes averted varied by a factor of 1.76 (95% CI: 1.54-2.06) across the range of uncertainty about breakthrough infections. This was equivalent to the range of vaccination impact projected across a range of uncertainty about vaccine efficacy of 0.268 (95% CI: 0.210-0.329). Until uncertainty about breakthrough infections can be addressed empirically, our results demonstrate the importance of accounting for it in models of vaccination impact.

Multi-scale computational modeling is a major branch of computational biology as evidenced by the US federal interagency Multi-Scale Modeling Consortium and major international projects. It invariably involves specific and detailed sequences of data analysis and simulation, often with multiple tools and datasets, and the community recognizes improved modularity, reuse, reproducibility, portability and scalability as critical unmet needs in this area. Scientific workflows are a well-recognized strategy for addressing these needs in scientific computing. While there are good examples if the use of scientific workflows in bioinformatics, medical informatics, biomedical imaging and data analysis, there are fewer examples in multi-scale computational modeling in general and cardiac electrophysiology in particular. Cardiac electrophysiology simulation is a mature area of multi-scale computational biology that serves as an excellent use case for developing and testing new scientific workflows. In this article, we develop, describe and test a computational workflow that serves as a proof of concept of a platform for the robust integration and implementation of a reusable and reproducible multi-scale cardiac cell and tissue model that is expandable, modular and portable. The workflow described leverages Python and Kepler-Python actor for plotting and pre/post-processing. During all stages of the workflow design, we rely on freely available open-source tools, to make our workflow freely usable by scientists.

Toll-like receptors (TLRs) play a central role in both the innate and adaptive immune systems by recognizing pathogen-associated molecular patterns and inducing the release of the effector molecules of the immune system. The dysregulation of the TLR system may cause various autoimmune diseases and septic shock. A series of molecular dynamics simulations and free energy calculations were performed to investigate the ligand-free, lipopolysaccharide (LPS)-bound, and neoseptin3-bound (TLR4-MD2)2 tetramers. Compared to earlier simulations done by others, our simulations showed that TLR4 structure was well maintained with stable interfaces. Free energy decomposition by molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method suggests critical roles that two hydrophobic clusters I85-L87-P88 and I124-L125-P127 of MD2, together with LPS and neoseptin3, may play in TLR4 activation. We propose that 1) direct contacts between TLR4 convex surface and LPS and neoseptin3 at the region around L442 significantly increase the binding and 2) binding of LPS and neoseptin3 in the central hydrophobic cavity of MD2 triggers burial of F126 and exposure of I85-L87-P88 that facilitate formation of (TLR4-MD2)2 tetramer and activation of TLR4 system.

An expanded chemical space is essential for improved identification of small molecules for emerging therapeutic targets. However, the identification of targets for novel compounds is biased towards the synthesis of known scaffolds that bind familiar protein families, limiting the exploration of chemical space. To change this paradigm, we validated a new pipeline that identifies small molecule-protein interactions and works even for compounds lacking similarity to known drugs. Based on differential mRNA profiles in multiple cell types exposed to drugs and in which gene knockdowns (KD) were conducted, we showed that drugs induce gene regulatory networks that correlate with those produced after silencing protein-coding genes. Next, we applied supervised machine learning to exploit drug-KD signature correlations and enriched our predictions using an orthogonal structure-based screen. As a proof-of-principle for this regimen, top-10/top-100 target prediction accuracies of 26% and 41%, respectively, were achieved on a validation of set 152 FDA-approved drugs and 3104 potential targets. We then predicted targets for 1680 compounds and validated chemical interactors with four targets that have proven difficult to chemically modulate, including non-covalent inhibitors of HRAS and KRAS. Importantly, drug-target interactions manifest as gene expression correlations between drug treatment and both target gene KD and KD of genes that act up- or down-stream of the target, even for relatively weak binders. These correlations provide new insights on the cellular response of disrupting protein interactions and highlight the complex genetic phenotypes of drug treatment. With further refinement, our pipeline may accelerate the identification and development of novel chemical classes by screening compound-target interactions.

How does pattern formation occur accurately when confronted with tissue growth and stochastic fluctuations (noise) in gene expression? Dorso-ventral (D-V) patterning of the mandibular arch specifies upper versus lower jaw skeletal elements through a combination of Bone morphogenetic protein (Bmp), Endothelin-1 (Edn1), and Notch signaling, and this system is highly robust. We combine NanoString experiments of early D-V gene expression with live imaging of arch development in zebrafish to construct a computational model of the D-V mandibular patterning network. The model recapitulates published genetic perturbations in arch development. Patterning is most sensitive to changes in Bmp signaling, and the temporal order of gene expression modulates the response of the patterning network to noise. Thus, our integrated systems biology approach reveals non-intuitive features of the complex signaling system crucial for craniofacial development, including novel insights into roles of gene expression timing and stochasticity in signaling and gene regulation.

Early afterdepolarizations (EADs) are spontaneous depolarizations during the repolarization phase of an action potential in cardiac myocytes. It is widely known that EADs are promoted by increasing inward currents and/or decreasing outward currents, a condition called reduced repolarization reserve. Recent studies based on bifurcation theories show that EADs are caused by a dual Hopf-homoclinic bifurcation, bringing in further mechanistic insights into the genesis and dynamics of EADs. In this study, we investigated the EAD properties, such as the EAD amplitude, the inter-EAD interval, and the latency of the first EAD, and their major determinants. We first made predictions based on the bifurcation theory and then validated them in physiologically more detailed action potential models. These properties were investigated by varying one parameter at a time or using parameter sets randomly drawn from assigned intervals. The theoretical and simulation results were compared with experimental data from the literature. Our major findings are that the EAD amplitude and takeoff potential exhibit a negative linear correlation; the inter-EAD interval is insensitive to the maximum ionic current conductance but mainly determined by the kinetics of ICa,L and the dual Hopf-homoclinic bifurcation; and both inter-EAD interval and latency vary largely from model to model. Most of the model results generally agree with experimental observations in isolated ventricular myocytes. However, a major discrepancy between modeling results and experimental observations is that the inter-EAD intervals observed in experiments are mainly between 200 and 500 ms, irrespective of species, while those of the mathematical models exhibit a much wider range with some models exhibiting inter-EAD intervals less than 100 ms. Our simulations show that the cause of this discrepancy is likely due to the difference in ICa,L recovery properties in different mathematical models, which needs to be addressed in future action potential model development.