bp: Blood pressure analysis in R

Despite the world-wide prevalence of hypertension, there is a lack in open-source software for analyzing blood pressure data. The R package bp fills this gap by providing functionality for blood pressure data processing, visualization, and feature extraction. In addition to the comprehensive functionality, the package includes six sample data sets covering continuous arterial pressure data (AP), home blood pressure monitoring data (HBPM) and ambulatory blood pressure monitoring data (ABPM), making it easier for researchers to get started. The R package bp is publicly available on CRAN and at https://github.com/johnschwenck/bp.

Highly scalable maximum likelihood and conjugate Bayesian inference for ERGMs on graph sets with equivalent vertices

The exponential family random graph modeling (ERGM) framework provides a highly flexible approach for the statistical analysis of networks (i.e., graphs). As ERGMs with dyadic dependence involve normalizing factors that are extremely costly to compute, practical strategies for ERGMs inference generally employ a variety of approximations or other workarounds. Markov Chain Monte Carlo maximum likelihood (MCMC MLE) provides a powerful tool to approximate the maximum likelihood estimator (MLE) of ERGM parameters, and is generally feasible for typical models on single networks with as many as a few thousand nodes. MCMC-based algorithms for Bayesian analysis are more expensive, and high-quality answers are challenging to obtain on large graphs. For both strategies, extension to the pooled case—in which we observe multiple networks from a common generative process—adds further computational cost, with both time and memory scaling linearly in the number of graphs. This becomes prohibitive for large networks, or cases in which large numbers of graph observations are available. Here, we exploit some basic properties of the discrete exponential families to develop an approach for ERGM inference in the pooled case that (where applicable) allows an arbitrarily large number of graph observations to be fit at no additional computational cost beyond preprocessing the data itself. Moreover, a variant of our approach can also be used to perform Bayesian inference under conjugate priors, again with no additional computational cost in the estimation phase. The latter can be employed either for single graph observations, or for observations from graph sets. As we show, the conjugate prior is easily specified, and is well-suited to applications such as regularization. Simulation studies show that the pooled method leads to estimates with good frequentist properties, and posterior estimates under the conjugate prior are well-behaved. We demonstrate the usefulness of our approach with applications to pooled analysis of brain functional connectivity networks and to replicated x-ray crystal structures of hen egg-white lysozyme.

Accelerated Optimization on Riemannian Manifolds via Discrete Constrained Variational Integrators

A variational formulation for accelerated optimization on normed vector spaces was recently introduced in Wibisono et al. (PNAS 113:E7351-E7358, 2016), and later generalized to the Riemannian manifold setting in Duruisseaux and Leok (SJMDS, 2022a). This variational framework was exploited on normed vector spaces in Duruisseaux et al. (SJSC 43:A2949-A2980, 2021) using time-adaptive geometric integrators to design efficient explicit algorithms for symplectic accelerated optimization, and it was observed that geometric discretizations which respect the time-rescaling invariance and symplecticity of the Lagrangian and Hamiltonian flows were substantially less prone to stability issues, and were therefore more robust, reliable, and computationally efficient. As such, it is natural to develop time-adaptive Hamiltonian variational integrators for accelerated optimization on Riemannian manifolds. In this paper, we consider the case of Riemannian manifolds embedded in a Euclidean space that can be characterized as the level set of a submersion. We will explore how holonomic constraints can be incorporated in discrete variational integrators to constrain the numerical discretization of the Riemannian Hamiltonian system to the Riemannian manifold, and we will test the performance of the resulting algorithms by solving eigenvalue and Procrustes problems formulated as optimization problems on the unit sphere and Stiefel manifold.

Data-Driven Modeling of Linear Dynamical Systems with Quadratic Output in the AAA Framework

We extend the Adaptive Antoulas-Anderson (AAA) algorithm to develop a data-driven modeling framework for linear systems with quadratic output (LQO). Such systems are characterized by two transfer functions: one corresponding to the linear part of the output and another one to the quadratic part. We first establish the joint barycentric representations and the interpolation theory for the two transfer functions of LQO systems. This analysis leads to the proposed AAA-LQO algorithm. We show that by interpolating the transfer function values on a subset of samples together with imposing a least-squares minimization on the rest, we construct reliable data-driven LQO models. Two numerical test cases illustrate the efficiency of the proposed method.

Individual heterogeneity influences the effects of translocation on urban dispersal of an invasive reptile

BACKGROUND

Invasive reptiles pose a serious threat to global biodiversity, but early detection of individuals in an incipient population is often hindered by their cryptic nature, sporadic movements, and variation among individuals. Little is known about the mechanisms that affect the movement of these species, which limits our understanding of their dispersal. Our aim was to determine whether translocation or small-scale landscape features affect movement patterns of brown treesnakes (Boiga irregularis), a destructive invasive predator on the island of Guam.

METHODS

We conducted a field experiment to compare the movements of resident (control) snakes to those of snakes translocated from forests and urban areas into new urban habitats. We developed a Bayesian hierarchical model to analyze snake movement mechanisms and account for attributes unique to invasive reptiles by incorporating multiple behavioral states and individual heterogeneity in movement parameters.

RESULTS

We did not observe strong differences in mechanistic movement parameters (turning angle or step length) among experimental treatment groups. We found some evidence that translocated snakes from both forests and urban areas made longer movements than resident snakes, but variation among individuals within treatment groups weakened this effect. Snakes translocated from forests moved more frequently from pavement than those translocated from urban areas. Snakes translocated from urban areas moved less frequently from buildings than resident snakes. Resident snakes had high individual heterogeneity in movement probability.

CONCLUSIONS

Our approach to modeling movement improved our understanding of invasive reptile dispersal by allowing us to examine the mechanisms that influence their movement. We also demonstrated the importance of accounting for individual heterogeneity in population-level analyses, especially when management goals involve eradication of an invasive species.

Mast seeding promotes evolution of scatter-hoarding

Many plant species worldwide are dispersed by scatter-hoarding granivores: animals that hide seeds in numerous, small caches for future consumption. Yet, the evolution of scatter-hoarding is difficult to explain because undefended caches are at high risk of pilferage. Previous models have attempted to solve this problem by giving cache owners large advantages in cache recovery, by kin selection, or by introducing reciprocal pilferage of 'shared' seed resources. However, the role of environmental variability has been so far overlooked in this context. One important form of such variability is masting, which is displayed by many plant species dispersed by scatterhoarders. We use a mathematical model to investigate the influence of masting on the evolution of scatter-hoarding. The model accounts for periodically varying annual seed fall, caching and pilfering behaviour, and the demography of scatterhoarders. The parameter values are based mostly on research on European beech (Fagus sylvatica) and yellow-necked mice (Apodemus flavicollis). Starvation of scatterhoarders between mast years decreases the population density that enters masting events, which leads to reduced seed pilferage. Satiation of scatterhoarders during mast events lowers the reproductive cost of caching (i.e. the cost of caching for the future rather than using seeds for current reproduction). These reductions promote the evolution of scatter-hoarding behaviour especially when interannual variation in seed fall and the period between masting events are large. This article is part of the theme issue 'The ecology and evolution of synchronized seed production in plants'.

Supracellular organization confers directionality and mechanical potency to migrating pairs of cardiopharyngeal progenitor cells

Physiological and pathological morphogenetic events involve a wide array of collective movements, suggesting that multicellular arrangements confer biochemical and biomechanical properties contributing to tissue-scale organization. The Ciona cardiopharyngeal progenitors provide the simplest model of collective cell migration, with cohesive bilateral cell pairs polarized along the leader-trailer migration path while moving between the ventral epidermis and trunk endoderm. We use the Cellular Potts Model to computationally probe the distributions of forces consistent with shapes and collective polarity of migrating cell pairs. Combining computational modeling, confocal microscopy, and molecular perturbations, we identify cardiopharyngeal progenitors as the simplest cell collective maintaining supracellular polarity with differential distributions of protrusive forces, cell-matrix adhesion, and myosin-based retraction forces along the leader-trailer axis. 4D simulations and experimental observations suggest that cell-cell communication helps establish a hierarchy to align collective polarity with the direction of migration, as observed with three or more cells in silico and in vivo. Our approach reveals emerging properties of the migrating collective: cell pairs are more persistent, migrating longer distances, and presumably with higher accuracy. Simulations suggest that cell pairs can overcome mechanical resistance of the trunk endoderm more effectively when they are polarized collectively. We propose that polarized supracellular organization of cardiopharyngeal progenitors confers emergent physical properties that determine mechanical interactions with their environment during morphogenesis.

mbImpute: an accurate and robust imputation method for microbiome data

A critical challenge in microbiome data analysis is the existence of many non-biological zeros, which distort taxon abundance distributions, complicate data analysis, and jeopardize the reliability of scientific discoveries. To address this issue, we propose the first imputation method for microbiome data-mbImpute-to identify and recover likely non-biological zeros by borrowing information jointly from similar samples, similar taxa, and optional metadata including sample covariates and taxon phylogeny. We demonstrate that mbImpute improves the power of identifying disease-related taxa from microbiome data of type 2 diabetes and colorectal cancer, and mbImpute preserves non-zero distributions of taxa abundances.

Scheduling fixed length quarantines to minimize the total number of fatalities during an epidemic

We consider a susceptible, infected, removed (SIR) system where the transmission rate may be temporarily reduced for a fixed amount of time. We show that in order to minimize the total number of fatalities, the transmission rate should be reduced on a single contiguous time interval, and we characterize this interval via an integral condition. We conclude with a few numerical simulations showing the actual reduction obtained.

Novel Autosegmentation Spatial Similarity Metrics Capture the Time Required to Correct Segmentations Better Than Traditional Metrics in a Thoracic Cavity Segmentation Workflow

Automated segmentation templates can save clinicians time compared to de novo segmentation but may still take substantial time to review and correct. It has not been thoroughly investigated which automated segmentation-corrected segmentation similarity metrics best predict clinician correction time. Bilateral thoracic cavity volumes in 329 CT scans were segmented by a UNet-inspired deep learning segmentation tool and subsequently corrected by a fourth-year medical student. Eight spatial similarity metrics were calculated between the automated and corrected segmentations and associated with correction times using Spearman’s rank correlation coefficients. Nine clinical variables were also associated with metrics and correction times using Spearman’s rank correlation coefficients or Mann–Whitney U tests. The added path length, false negative path length, and surface Dice similarity coefficient correlated better with correction time than traditional metrics, including the popular volumetric Dice similarity coefficient (respectively ρ = 0.69, ρ = 0.65, ρ =  − 0.48 versus ρ =  − 0.25; correlation p values < 0.001). Clinical variables poorly represented in the autosegmentation tool’s training data were often associated with decreased accuracy but not necessarily with prolonged correction time. Metrics used to develop and evaluate autosegmentation tools should correlate with clinical time saved. To our knowledge, this is only the second investigation of which metrics correlate with time saved. Validation of our findings is indicated in other anatomic sites and clinical workflows. Novel spatial similarity metrics may be preferable to traditional metrics for developing and evaluating autosegmentation tools that are intended to save clinicians time.

Three pre-vaccine responses to Covid-like epidemics

This paper contains a theoretical study of epidemic control. It is inspired by current events but not intended to be an accurate depiction of the SARS-CoV-2 pandemic. We consider the emergence of a highly transmissible pathogen, focusing on metropolitan areas. To ensure some degree of realism, we present a conceptual model of the outbreak and early attempts to stave off the onslaught, including the use of lockdowns. Model outputs show strong qualitative—in some respects even quantitative—resemblance to the events of Spring 2020 in many cities worldwide. We then use this model to project forward in time to examine different paths in epidemic control after the initial surge is tamed and before the arrival of vaccines. Three very different control strategies are analyzed, leading to vastly different outcomes in terms of economic recovery and total infected population (or progress toward herd immunity). Our model, which is a version of the SEIQR model, is a time-dependent dynamical system with feedback-control. One of the main conclusions of this analysis is that the course of the epidemic is not entirely dictated by the virus: how the population responds to it can play an equally important role in determining the eventual outcome.

Incorporating false negative tests in epidemiological models for SARS-CoV-2 transmission and reconciling with seroprevalence estimates

Susceptible-Exposed-Infected-Removed (SEIR)-type epidemiologic models, modeling unascertained infections latently, can predict unreported cases and deaths assuming perfect testing. We apply a method we developed to account for the high false negative rates of diagnostic RT-PCR tests for detecting an active SARS-CoV-2 infection in a classic SEIR model. The number of unascertained cases and false negatives being unobservable in a real study, population-based serosurveys can help validate model projections. Applying our method to training data from Delhi, India, during March 15-June 30, 2020, we estimate the underreporting factor for cases at 34-53 (deaths: 8-13) on July 10, 2020, largely consistent with the findings of the first round of serosurveys for Delhi (done during June 27-July 10, 2020) with an estimated 22.86% IgG antibody prevalence, yielding estimated underreporting factors of 30-42 for cases. Together, these imply approximately 96-98% cases in Delhi remained unreported (July 10, 2020). Updated calculations using training data during March 15-December 31, 2020 yield estimated underreporting factor for cases at 13-22 (deaths: 3-7) on January 23, 2021, which are again consistent with the latest (fifth) round of serosurveys for Delhi (done during January 15-23, 2021) with an estimated 56.13% IgG antibody prevalence, yielding an estimated range for the underreporting factor for cases at 17-21. Together, these updated estimates imply approximately 92-96% cases in Delhi remained unreported (January 23, 2021). Such model-based estimates, updated with latest data, provide a viable alternative to repeated resource-intensive serosurveys for tracking unreported cases and deaths and gauging the true extent of the pandemic.

Using excess deaths and testing statistics to determine COVID-19 mortalities

Factors such as varied definitions of mortality, uncertainty in disease prevalence, and biased sampling complicate the quantification of fatality during an epidemic. Regardless of the employed fatality measure, the infected population and the number of infection-caused deaths need to be consistently estimated for comparing mortality across regions. We combine historical and current mortality data, a statistical testing model, and an SIR epidemic model, to improve estimation of mortality. We find that the average excess death across the entire US from January 2020 until February 2021 is 9[Formula: see text] higher than the number of reported COVID-19 deaths. In some areas, such as New York City, the number of weekly deaths is about eight times higher than in previous years. Other countries such as Peru, Ecuador, Mexico, and Spain exhibit excess deaths significantly higher than their reported COVID-19 deaths. Conversely, we find statistically insignificant or even negative excess deaths for at least most of 2020 in places such as Germany, Denmark, and Norway.

Unified model of short- and long-term HIV viral rebound for clinical trial planning

Antiretroviral therapy (ART) effectively controls HIV infection, suppressing HIV viral loads. Typically suspension of therapy is rapidly followed by rebound of viral loads to high, pre-therapy levels. Indeed, a recent study showed that approximately 90% of treatment interruption study participants show viral rebound within at most a few months of therapy suspension, but the remaining 10%, showed viral rebound some months, or years, after ART suspension. Some may even never rebound. We investigate and compare branching process models aimed at gaining insight into these viral dynamics. Specifically, we provide a theory that explains both short- and long-term viral rebounds, and post-treatment control, via a multitype branching process with time-inhomogeneous rates, validated with data from Li et al. (Li et al. 2016 AIDS30, 343-353. (doi:10.1097/QAD.0000000000000953)). We discuss the associated biological interpretation and implications of our best-fit model. To test the effectiveness of an experimental intervention in delaying or preventing rebound, the standard practice is to suspend therapy and monitor the study participants for rebound. We close with a discussion of an important application of our modelling in the design of such clinical trials.

Powerline right-of-way management and flower-visiting insects: How vegetation management can promote pollinator diversity

Loss in the availability of early successional habitat is a threat to pollinator populations. Given that powerline rights-of-way (ROW) must be managed to maintain early successional habitat, preventing vegetation from interfering with electrical lines, they have the potential to provide conservation benefits for wild pollinators. Moreover, it is possible to provide conservation benefits with no additional cost to land managers. We surveyed flower-visiting insects over two years in different vegetation management treatments in a long-term research ROW to determine which best promoted pollinator abundance and species richness. We found that the ROW had stabilized in an early successional state soon after its establishment and that this early successional state could be maintained with low levels of periodic maintenance. We collected a high diversity of flower-visiting insects (126 bee species and 179 non-bee morphospecies) in six ROW plots. Higher levels of herbicide application had a negative effect on bee species richness, but low levels of herbicide application were compatible with a high abundance and species richness of flower-visiting insects, including several rare species. Moreover, this effect was seen only in the bee community, and not in non-bee flower-visiting insects. Our results suggest further research into the conservation value of ROW for pollinators is warranted. We demonstrate that there is substantial potential for pollinator conservation in ROW, compatible with low-cost vegetation management.

Decoding SARS-CoV-2 Transmission and Evolution and Ramifications for COVID-19 Diagnosis, Vaccine, and Medicine

Tremendous effort has been given to the development of diagnostic tests, preventive vaccines, and therapeutic medicines for coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Much of this development has been based on the reference genome collected on January 5, 2020. Based on the genotyping of 15 140 genome samples collected up to June 1, 2020, we report that SARS-CoV-2 has undergone 8309 single mutations which can be clustered into six subtypes. We introduce mutation ratio and mutation h-index to characterize the protein conservativeness and unveil that SARS-CoV-2 envelope protein, main protease, and endoribonuclease protein are relatively conservative, while SARS-CoV-2 nucleocapsid protein, spike protein, and papain-like protease are relatively nonconservative. In particular, we have identified mutations on 40% of nucleotides in the nucleocapsid gene in the population level, signaling potential impacts on the ongoing development of COVID-19 diagnosis, vaccines, and antibody and small-molecular drugs.

A duality between scattering poles and transmission eigenvalues in scattering theory

In this paper, we develop a conceptually unified approach for characterizing and determining scattering poles and interior eigenvalues for a given scattering problem. Our approach explores a duality stemming from interchanging the roles of incident and scattered fields in our analysis. Both sets are related to the kernel of the relative scattering operator mapping incident fields to scattered fields, corresponding to the exterior scattering problem for the interior eigenvalues and the interior scattering problem for scattering poles. Our discussion includes the scattering problem for a Dirichlet obstacle where duality is between scattering poles and Dirichlet eigenvalues, and the inhomogeneous scattering problem where the duality is between scattering poles and transmission eigenvalues. Our new characterization of the scattering poles suggests a numerical method for their computation in terms of scattering data for the corresponding interior scattering problem.

Strong intracellular signal inactivation produces sharper and more robust signaling from cell membrane to nucleus

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.

The inside of cells is a complex spatial environment, filled with organelles, filaments and proteins. It is an open question how cell signaling pathways function robustly in the presence of such spatial heterogeneity. In this work we study how organelle barriers influence the most basic of chemical signals; the diffusive propagation of an activated protein from the cell membrane to nucleus. Three-dimensional B cell organelle and membrane geometries reconstructed from soft X-ray tomographic images are used in building mathematical models of the signal propagation process. Our models demonstrate that organelle barriers significantly increase the time required for a diffusing protein to traverse from the cell membrane to nucleus when compared to a cell with an empty cytosolic space. We also show that signal inactivation, a fundamental component of all signaling pathways, can provide robustness in the signal arrival time in two ways. Increasing rates of signal inactivation reduce variability in the arrival time, while also dramatically reducing the degree to which organelle barriers increase the arrival time (in comparison to a cell with an empty cytosol).

PhANNs, a fast and accurate tool and web server to classify phage structural proteins

For any given bacteriophage genome or phage-derived sequences in metagenomic data sets, we are unable to assign a function to 50–90% of genes, or more. Structural protein-encoding genes constitute a large fraction of the average phage genome and are among the most divergent and difficult-to-identify genes using homology-based methods. To understand the functions encoded by phages, their contributions to their environments, and to help gauge their utility as potential phage therapy agents, we have developed a new approach to classify phage ORFs into ten major classes of structural proteins or into an “other” category. The resulting tool is named PhANNs (Phage Artificial Neural Networks). We built a database of 538,213 manually curated phage protein sequences that we split into eleven subsets (10 for cross-validation, one for testing) using a novel clustering method that ensures there are no homologous proteins between sets yet maintains the maximum sequence diversity for training. An Artificial Neural Network ensemble trained on features extracted from those sets reached a test F₁-score of 0.875 and test accuracy of 86.2%. PhANNs can rapidly classify proteins into one of the ten structural classes or, if not predicted to fall in one of the ten classes, as “other,” providing a new approach for functional annotation of phage proteins. PhANNs is open source and can be run from our web server or installed locally.

Bacteriophages (phages, viruses that infect bacteria) are the most abundant biological entity on Earth. They outnumber bacteria by a factor of ten. As phages are very different from each other and from bacteria, and we have relatively few phage genes in our database compared to bacterial genes, we are unable to assign function to 50–90% of phage genes. In this work, we developed PhANNs, a machine learning tool that can classify a phage gene as one of ten structural roles, or “other”. This approach does not require a similar gene to be known.

Estimating effective population size changes from preferentially sampled genetic sequences

Coalescent theory combined with statistical modeling allows us to estimate effective population size fluctuations from molecular sequences of individuals sampled from a population of interest. When sequences are sampled serially through time and the distribution of the sampling times depends on the effective population size, explicit statistical modeling of sampling times improves population size estimation. Previous work assumed that the genealogy relating sampled sequences is known and modeled sampling times as an inhomogeneous Poisson process with log-intensity equal to a linear function of the log-transformed effective population size. We improve this approach in two ways. First, we extend the method to allow for joint Bayesian estimation of the genealogy, effective population size trajectory, and other model parameters. Next, we improve the sampling time model by incorporating additional sources of information in the form of time-varying covariates. We validate our new modeling framework using a simulation study and apply our new methodology to analyses of population dynamics of seasonal influenza and to the recent Ebola virus outbreak in West Africa.

Estimating changes in the number of individuals in a given population is a challenging problem in some settings. For example, estimating population size trajectories of the number of people infected by a pathogen (e.g., Influenza virus) is a difficult problem, because many infections in a large population remain unobserved/hidden. One indirect way of assessing population size changes is to take a sample of individuals from the population of interest and analyze genetic sequences from these individuals (e.g., Influenza virus genomes). Intuitively, genetic data is informative about population size changes, because genetic diversity increases/decreases together with the population size. However, if we sample more individuals when the population size increases and less when it decreases, this strategy produces biased results. To avoid this bias, we propose a method that explicitly and flexibly models potential dependency of genetic sequence sampling on the population size. An added bonus of this new modeling framework is more precise estimation of population size changes. We demonstrate strengths of our new methodology on simulated data and on genetic sequences of Influenza and Ebola viruses.

Rapid prediction of crucial hotspot interactions for icosahedral viral capsid self-assembly by energy landscape atlasing validated by mutagenesis

Icosahedral viruses are under a micrometer in diameter, their infectious genome encapsulated by a shell assembled by a multiscale process, starting from an integer multiple of 60 viral capsid or coat protein (VP) monomers. We predict and validate inter-atomic hotspot interactions between VP monomers that are important for the assembly of 3 types of icosahedral viral capsids: Adeno Associated Virus serotype 2 (AAV2) and Minute Virus of Mice (MVM), both T = 1 single stranded DNA viruses, and Bromo Mosaic Virus (BMV), a T = 3 single stranded RNA virus. Experimental validation is by in-vitro, site-directed mutagenesis data found in literature. We combine ab-initio predictions at two scales: at the interface-scale, we predict the importance (cruciality) of an interaction for successful subassembly across each interface between symmetry-related VP monomers; and at the capsid-scale, we predict the cruciality of an interface for successful capsid assembly. At the interface-scale, we measure cruciality by changes in the capsid free-energy landscape partition function when an interaction is removed. The partition function computation uses atlases of interface subassembly landscapes, rapidly generated by a novel geometric method and curated opensource software EASAL (efficient atlasing and search of assembly landscapes). At the capsid-scale, cruciality of an interface for successful assembly of the capsid is based on combinatorial entropy. Our study goes all the way from resource-light, multiscale computational predictions of crucial hotspot inter-atomic interactions to validation using data on site-directed mutagenesis' effect on capsid assembly. By reliably and rapidly narrowing down target interactions, (no more than 1.5 hours per interface on a laptop with Intel Core i5-2500K @ 3.2 Ghz CPU and 8GB of RAM) our predictions can inform and reduce time-consuming in-vitro and in-vivo experiments, or more computationally intensive in-silico analyses.

Renewal Reward Perspective on Linear Switching Diffusion Systems in Models of Intracellular Transport

In many biological systems, the movement of individual agents is characterized having multiple qualitatively distinct behaviors that arise from a variety of biophysical states. For example, in cells the movement of vesicles, organelles, and other intracellular cargo is affected by their binding to and unbinding from cytoskeletal filaments such as microtubules through molecular motor proteins. A typical goal of theoretical or numerical analysis of models of such systems is to investigate effective transport properties and their dependence on model parameters. While the effective velocity of particles undergoing switching diffusion dynamics is often easily characterized in terms of the long-time fraction of time that particles spend in each state, the calculation of the effective diffusivity is more complicated because it cannot be expressed simply in terms of a statistical average of the particle transport state at one moment of time. However, it is common that these systems are regenerative, in the sense that they can be decomposed into independent cycles marked by returns to a base state. Using decompositions of this kind, we calculate effective transport properties by computing the moments of the dynamics within each cycle and then applying renewal reward theory. This method provides a useful alternative large-time analysis to direct homogenization for linear advection-reaction-diffusion partial differential equation models. Moreover, it applies to a general class of semi-Markov processes and certain stochastic differential equations that arise in models of intracellular transport. Applications of the proposed renewal reward framework are illustrated for several case studies such as mRNA transport in developing oocytes and processive cargo movement by teams of molecular motor proteins.

The de Rham-Hodge Analysis and Modeling of Biomolecules

Biological macromolecules have intricate structures that underpin their biological functions. Understanding their structure-function relationships remains a challenge due to their structural complexity and functional variability. Although de Rham-Hodge theory, a landmark of twentieth-century mathematics, has had a tremendous impact on mathematics and physics, it has not been devised for macromolecular modeling and analysis. In this work, we introduce de Rham-Hodge theory as a unified paradigm for analyzing the geometry, topology, flexibility, and Hodge mode analysis of biological macromolecules. Geometric characteristics and topological invariants are obtained either from the Helmholtz-Hodge decomposition of the scalar, vector, and/or tensor fields of a macromolecule or from the spectral analysis of various Laplace-de Rham operators defined on the molecular manifolds. We propose Laplace-de Rham spectral-based models for predicting macromolecular flexibility. We further construct a Laplace-de Rham-Helfrich operator for revealing cryo-EM natural frequencies. Extensive experiments are carried out to demonstrate that the proposed de Rham-Hodge paradigm is one of the most versatile tools for the multiscale modeling and analysis of biological macromolecules and subcellular organelles. Accurate, reliable, and topological structure-preserving algorithms for implementing discrete exterior calculus (DEC) have been developed to facilitate the aforementioned modeling and analysis of biological macromolecules. The proposed de Rham-Hodge paradigm has potential applications to subcellular organelles and the structure construction from medium- or low-resolution cryo-EM maps, and functional predictions from massive biomolecular datasets.

Collective vibrations of a hydrodynamic active lattice

Recent experiments show that quasi-one-dimensional lattices of self-propelled droplets exhibit collective instabilities in the form of out-of-phase oscillations and solitary-like waves. This hydrodynamic lattice is driven by the external forcing of a vertically vibrating fluid bath, which invokes a field of subcritical Faraday waves on the bath surface, mediating the spatio-temporal droplet coupling. By modelling the droplet lattice as a memory-endowed system with spatially non-local coupling, we herein rationalize the form and onset of instability in this new class of dynamical oscillator. We identify the memory-driven instability of the lattice as a function of the number of droplets, and determine equispaced lattice configurations precluded by geometrical constraints. Each memory-driven instability is then classified as either a super- or subcritical Hopf bifurcation via a systematic weakly nonlinear analysis, rationalizing experimental observations. We further discover a previously unreported symmetry-breaking instability, manifest as an oscillatory-rotary motion of the lattice. Numerical simulations support our findings and prompt further investigations of this nonlinear dynamical system.

A Note on Exploratory Item Factor Analysis by Singular Value Decomposition

We revisit a singular value decomposition (SVD) algorithm given in Chen et al. (Psychometrika 84:124-146, 2019b) for exploratory item factor analysis (IFA). This algorithm estimates a multidimensional IFA model by SVD and was used to obtain a starting point for joint maximum likelihood estimation in Chen et al. (2019b). Thanks to the analytic and computational properties of SVD, this algorithm guarantees a unique solution and has computational advantage over other exploratory IFA methods. Its computational advantage becomes significant when the numbers of respondents, items, and factors are all large. This algorithm can be viewed as a generalization of principal component analysis to binary data. In this note, we provide the statistical underpinning of the algorithm. In particular, we show its statistical consistency under the same double asymptotic setting as in Chen et al. (2019b). We also demonstrate how this algorithm provides a scree plot for investigating the number of factors and provide its asymptotic theory. Further extensions of the algorithm are discussed. Finally, simulation studies suggest that the algorithm has good finite sample performance.

Simultaneous confidence corridors for mean functions in functional data analysis of imaging data

Motivated by recent work involving the analysis of biomedical imaging data, we present a novel procedure for constructing simultaneous confidence corridors for the mean of imaging data. We propose to use flexible bivariate splines over triangulations to handle an irregular domain of the images that is common in brain imaging studies and in other biomedical imaging applications. The proposed spline estimators of the mean functions are shown to be consistent and asymptotically normal under some regularity conditions. We also provide a computationally efficient estimator of the covariance function and derive its uniform consistency. The procedure is also extended to the two-sample case in which we focus on comparing the mean functions from two populations of imaging data. Through Monte Carlo simulation studies, we examine the finite sample performance of the proposed method. Finally, the proposed method is applied to analyze brain positron emission tomography data in two different studies. One data set used in preparation of this article was obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database.

Modeling Stripe Formation on Growing Zebrafish Tailfins

As zebrafish develop, black and gold stripes form across their skin due to the interactions of brightly colored pigment cells. These characteristic patterns emerge on the growing fish body, as well as on the anal and caudal fins. While wild-type stripes form parallel to a horizontal marker on the body, patterns on the tailfin gradually extend distally outward. Interestingly, several mutations lead to altered body patterns without affecting fin stripes. Through an exploratory modeling approach, our goal is to help better understand these differences between body and fin patterns. By adapting a prior agent-based model of cell interactions on the fish body, we present an in silico study of stripe development on tailfins. Our main result is a demonstration that two cell types can produce stripes on the caudal fin. We highlight several ways that bone rays, growth, and the body-fin interface may be involved in patterning, and we raise questions for future work related to pattern robustness.

Pointwise Bounds for Joint Eigenfunctions of Quantum Completely Integrable Systems

Let (M, g) be a compact Riemannian manifold of dimension n and P 1 : = - h 2 Δ g + V ( x ) - E 1 so that d p 1 ≠ 0 on p 1 = 0 . We assume that P 1 is quantum completely integrable (ACI) in the sense that there exist functionally independent pseuodifferential operators P 2 , ⋯ P n with [ P i , P j ] = 0 , i , j = 1 , ⋯ n . We study the pointwise bounds for the joint eigenfunctions, u h of the system { P i } i = 1 n with P 1 u h = E 1 u h + o ( 1 ) . In Theorem 1, we first give polynomial improvements over the standard Hörmander bounds for typical points in M. In two and three dimensions, these estimates agree with the Hardy exponent h - 1 - n 4 and in higher dimensions we obtain a gain of h 1 2 over the Hörmander bound. In our second main result (Theorem 3), under a real-analyticity assumption on the QCI system, we give exponential decay estimates for joint eigenfunctions at points outside the projection of invariant Lagrangian tori; that is at points x ∈ M in the "microlocally forbidden" region p 1 - 1 ( E 1 ) ∩ ⋯ ∩ p n - 1 ( E n ) ∩ T x ∗ M = ∅ . These bounds are sharp locally near the projection of the invariant tori.

3D mesh processing using GAMer 2 to enable reaction-diffusion simulations in realistic cellular geometries

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.

Striated myocyte structural integrity: Automated analysis of sarcomeric z-discs

As sarcomeres produce the force necessary for contraction, assessment of sarcomere order is paramount in evaluation of cardiac and skeletal myocytes. The uniaxial force produced by sarcomeres is ideally perpendicular to their z-lines, which couple parallel myofibrils and give cardiac and skeletal myocytes their distinct striated appearance. Accordingly, sarcomere structure is often evaluated by staining for z-line proteins such as α-actinin. However, due to limitations of current analysis methods, which require manual or semi-manual handling of images, the mechanism by which sarcomere and by extension z-line architecture can impact contraction and which characteristics of z-line architecture should be used to assess striated myocytes has not been fully explored. Challenges such as isolating z-lines from regions of off-target staining that occur along immature stress fibers and cell boundaries and choosing metrics to summarize overall z-line architecture have gone largely unaddressed in previous work. While an expert can qualitatively appraise tissues, these challenges leave researchers without robust, repeatable tools to assess z-line architecture across different labs and experiments. Additionally, the criteria used by experts to evaluate sarcomeric architecture have not been well-defined. We address these challenges by providing metrics that summarize different aspects of z-line architecture that correspond to expert tissue quality assessment and demonstrate their efficacy through an examination of engineered tissues and single cells. In doing so, we have elucidated a mechanism by which highly elongated cardiomyocytes become inefficient at producing force. Unlike previous manual or semi-manual methods, characterization of z-line architecture using the metrics discussed and implemented in this work can quantitatively evaluate engineered tissues and contribute to a robust understanding of the development and mechanics of striated muscles.

Multistability in the epithelial-mesenchymal transition network

BACKGROUND

The transitions between epithelial (E) and mesenchymal (M) cell phenotypes are essential in many biological processes like tissue development and cancer metastasis. Previous studies, both modeling and experimental, suggested that in addition to E and M states, the network responsible for these phenotypes exhibits intermediate phenotypes between E and M states. The number and importance of such states is subject to intense discussion in the epithelial-mesenchymal transition (EMT) community.

RESULTS

Previous modeling efforts used traditional bifurcation analysis to explore the number of the steady states that correspond to E, M and intermediate states by varying one or two parameters at a time. Since the system has dozens of parameters that are largely unknown, it remains a challenging problem to fully describe the potential set of states and their relationship across all parameters. We use the computational tool DSGRN (Dynamic Signatures Generated by Regulatory Networks) to explore the intermediate states of an EMT model network by computing summaries of the dynamics across all of parameter space. We find that the only attractors in the system are equilibria, that E and M states dominate across parameter space, but that bistability and multistability are common. Even at extreme levels of some of the known inducers of the transition, there is a certain proportion of the parameter space at which an E or an M state co-exists with other stable steady states.

CONCLUSIONS

Our results suggest that the multistability is broadly present in the EMT network across parameters and thus response of cells to signals may strongly depend on the particular cell line and genetic background.

Transient amplifiers of selection and reducers of fixation for death-Birth updating on graphs

The spatial structure of an evolving population affects the balance of natural selection versus genetic drift. Some structures amplify selection, increasing the role that fitness differences play in determining which mutations become fixed. Other structures suppress selection, reducing the effect of fitness differences and increasing the role of random chance. This phenomenon can be modeled by representing spatial structure as a graph, with individuals occupying vertices. Births and deaths occur stochastically, according to a specified update rule. We study death-Birth updating: An individual is chosen to die and then its neighbors compete to reproduce into the vacant spot. Previous numerical experiments suggested that amplifiers of selection for this process are either rare or nonexistent. We introduce a perturbative method for this problem for weak selection regime, meaning that mutations have small fitness effects. We show that fixation probability under weak selection can be calculated in terms of the coalescence times of random walks. This result leads naturally to a new definition of effective population size. Using this and other methods, we uncover the first known examples of transient amplifiers of selection (graphs that amplify selection for a particular range of fitness values) for the death-Birth process. We also exhibit new families of “reducers of fixation”, which decrease the fixation probability of all mutations, whether beneficial or deleterious.

Natural selection is often thought of as “survival of the fittest”, but random chance plays a significant role in which mutations persist and which are eliminated. The balance of selection versus randomness is affected by spatial structure—how individuals are arranged within their habitat. Some structures amplify the effects of selection, so that only the fittest mutations are likely to persist. Others suppress the effects of selection, making the survival of genes primarily a matter of random chance. We study this question using a mathematical model called the “death-Birth process”. Previous studies have found that spatial structure rarely, if ever, amplifies selection for this process. Here we report that spatial structure can indeed amplify selection, at least for mutations with small fitness effects. We also identify structures that reduce the spread of any new mutation, whether beneficial or deleterious. Our work introduces new mathematical techniques for assessing how population structure affects natural selection.

Identifiability and numerical algebraic geometry

A common problem when analyzing models, such as mathematical modeling of a biological process, is to determine if the unknown parameters of the model can be determined from given input-output data. Identifiable models are models such that the unknown parameters can be determined to have a finite number of values given input-output data. The total number of such values over the complex numbers is called the identifiability degree of the model. Unidentifiable models are models such that the unknown parameters can have an infinite number of values given input-output data. For unidentifiable models, a set of identifiable functions of the parameters are sought so that the model can be reparametrized in terms of these functions yielding an identifiable model. In this work, we use numerical algebraic geometry to determine if a model given by polynomial or rational ordinary differential equations is identifiable or unidentifiable. For identifiable models, we present a novel approach to compute the identifiability degree. For unidentifiable models, we present a novel numerical differential algebra technique aimed at computing a set of algebraically independent identifiable functions. Several examples are used to demonstrate the new techniques.

A computational study on the role of glutamate and NMDA receptors on cortical spreading depression using a multidomain electrodiffusion model

Cortical spreading depression (SD) is a spreading disruption of ionic homeostasis in the brain during which neurons experience complete and prolonged depolarizations. SD is the basis of migraine aura and is increasingly associated with many other brain pathologies. Here, we study the role of glutamate and NMDA receptor dynamics in the context of an ionic electrodiffusion model. We perform simulations in one (1D) and two (2D) spatial dimension. Our 1D simulations reproduce the "inverted saddle" shape of the extracellular voltage signal for the first time. Our simulations suggest that SD propagation depends on two overlapping mechanisms; one dependent on extracellular glutamate diffusion and NMDA receptors and the other dependent on extracellular potassium diffusion and persistent sodium channel conductance. In 2D simulations, we study the dynamics of spiral waves. We study the properties of the spiral waves in relation to the planar 1D wave, and also compute the energy expenditure associated with the recurrent SD spirals.

Social inhibition maintains adaptivity and consensus of honeybees foraging in dynamic environments

To effectively forage in natural environments, organisms must adapt to changes in the quality and yield of food sources across multiple timescales. Individuals foraging in groups act based on both their private observations and the opinions of their neighbours. How do these information sources interact in changing environments? We address this problem in the context of honeybee colonies whose inhibitory social interactions promote adaptivity and consensus needed for effective foraging. Individual and social interactions within a mathematical model of collective decisions shape the nutrition yield of a group foraging from feeders with temporally switching quality. Social interactions improve foraging from a single feeder if temporal switching is fast or feeder quality is low. When the colony chooses from multiple feeders, the most beneficial form of social interaction is direct switching, whereby bees flip the opinion of nest-mates foraging at lower-yielding feeders. Model linearization shows that effective social interactions increase the fraction of the colony at the correct feeder (consensus) and the rate at which bees reach that feeder (adaptivity). Our mathematical framework allows us to compare a suite of social inhibition mechanisms, suggesting experimental protocols for revealing effective colony foraging strategies in dynamic environments.

Isolating phyllotactic patterns embedded in the secondary growth of sweet cherry (Prunus avium L.) using magnetic resonance imaging

BACKGROUND

Epicormic branches arise from dormant buds patterned during the growth of previous years. Dormant epicormic buds remain just below the surface of trees, pushed outward from the pith during secondary growth, but maintain vascular connections. Epicormic buds can be activated to elongate into a new shoot, either through natural processes or horticultural intervention, to potentially rejuvenate orchards and restructure tree architecture. Because epicormic structures are embedded within secondary growth, tomographic approaches are a useful method to study them and understand their development.

RESULTS

We apply techniques from image processing to determine the locations of epicormic vascular traces embedded within secondary growth of sweet cherry (Prunus avium L.), revealing the juvenile phyllotactic pattern in the trunk of an adult tree. Techniques include the flood fill algorithm to find the pith of the tree, edge detection to approximate the radius, and a conversion to polar coordinates to threshold and segment phyllotactic features. Intensity values from magnetic resonance imaging (MRI) of the trunk are projected onto the surface of a perfect cylinder to find the locations of traces in the "boundary image". Mathematical phyllotaxy provides a means to capture the patterns in the boundary image by modeling phyllotactic parameters. Our cherry tree specimen has the conspicuous parastichy pair (2,3), phyllotactic fraction 2/5, and divergence angle of approximately 143°.

CONCLUSIONS

The methods described provide a framework not only for studying phyllotaxy, but also for processing of volumetric image data in plants. Our results have practical implications for orchard rejuvenation and directed approaches to influence tree architecture. The study of epicormic structures, which are hidden within secondary growth, using tomographic methods also opens the possibility of studying genetic and environmental influences such structures.

Detection of arterial wall abnormalities via Bayesian model selection

Patient-specific modelling of haemodynamics in arterial networks has so far relied on parameter estimation for inexpensive or small-scale models. We describe here a Bayesian uncertainty quantification framework which makes two major advances: an efficient parallel implementation, allowing parameter estimation for more complex forward models, and a system for practical model selection, allowing evidence-based comparison between distinct physical models. We demonstrate the proposed methodology by generating simulated noisy flow velocity data from a branching arterial tree model in which a structural defect is introduced at an unknown location; our approach is shown to accurately locate the abnormality and estimate its physical properties even in the presence of significant observational and systemic error. As the method readily admits real data, it shows great potential in patient-specific parameter fitting for haemodynamical flow models.

Effects of persistent sodium current blockade in respiratory circuits depend on the pharmacological mechanism of action and network dynamics

The mechanism(s) of action of most commonly used pharmacological blockers of voltage-gated ion channels are well understood; however, this knowledge is rarely considered when interpreting experimental data. Effects of blockade are often assumed to be equivalent, regardless of the mechanism of the blocker involved. Using computer simulations, we demonstrate that this assumption may not always be correct. We simulate the blockade of a persistent sodium current (I_NaP), proposed to underlie rhythm generation in pre-Bötzinger complex (pre-BötC) respiratory neurons, via two distinct pharmacological mechanisms: (1) pore obstruction mediated by tetrodotoxin and (2) altered inactivation dynamics mediated by riluzole. The reported effects of experimental application of tetrodotoxin and riluzole in respiratory circuits are diverse and seemingly contradictory and have led to considerable debate within the field as to the specific role of I_NaP in respiratory circuits. The results of our simulations match a wide array of experimental data spanning from the level of isolated pre-BötC neurons to the level of the intact respiratory network and also generate a series of experimentally testable predictions. Specifically, in this study we: (1) provide a mechanistic explanation for seemingly contradictory experimental results from in vitro studies of I_NaP block, (2) show that the effects of I_NaP block in in vitro preparations are not necessarily equivalent to those in more intact preparations, (3) demonstrate and explain why riluzole application may fail to effectively block I_NaP in the intact respiratory network, and (4) derive the prediction that effective block of I_NaP by low concentration tetrodotoxin will stop respiratory rhythm generation in the intact respiratory network. These simulations support a critical role for I_NaP in respiratory rhythmogenesis in vivo and illustrate the importance of considering mechanism when interpreting and simulating data relating to pharmacological blockade.

The application of pharmacological agents that affect transmembrane ionic currents in neurons is a commonly used experimental technique. A simplistic interpretation of experiments involving these agents suggests that antagonist application removes the impacted current and that subsequently observed changes in activity are attributable to the loss of that current’s effects. The more complex reality, however, is that different drugs may have distinct mechanisms of action, some corresponding not to a removal of a current but rather to a changing of its properties. We use computational modeling to explore the implications of the distinct mechanisms associated with two drugs, riluzole and tetrodotoxin, that are often characterized as sodium channel blockers. Through this approach, we offer potential explanations for disparate findings observed in experiments on neural respiratory circuits and show that the experimental results are consistent with a key role for the persistent sodium current in respiratory rhythm generation.

Invariant Measures for the Stochastic One-Dimensional Compressible Navier-Stokes Equations

We investigate the long-time behavior of solutions to a stochastically forced one-dimensional Navier-Stokes system, describing the motion of a compressible viscous fluid, in the case of linear pressure law. We prove existence of an invariant measure for the Markov process generated by strong solutions. We overcome the difficulties of working with non-Feller Markov semigroups on non-complete metric spaces by generalizing the classical Krylov-Bogoliubov method, and by providing suitable polynomial and exponential moment bounds on the solution, together with pathwise estimates.

Predictions of time to HIV viral rebound following ART suspension that incorporate personal biomarkers

Antiretroviral therapy (ART) effectively controls HIV infection, suppressing HIV viral loads. Suspension of therapy is followed by rebound of viral loads to high, pre-therapy levels. However, there is significant heterogeneity in speed of rebound, with some rebounds occurring within days, weeks, or sometimes years. We present a stochastic mathematical model to gain insight into these post-treatment dynamics, specifically characterizing the dynamics of short term viral rebounds (≤ 60 days). Li et al. (2016) report that the size of the expressed HIV reservoir, i.e., cell-associated HIV RNA levels, and drug regimen correlate with the time between ART suspension and viral rebound to detectable levels. We incorporate this information and viral rebound times to parametrize our model. We then investigate insights offered by our model into the underlying dynamics of the latent reservoir. In particular, we refine previous estimates of viral recrudescence after ART interruption by accounting for heterogeneity in infection rebound dynamics, and determine a recrudescence rate of once every 2-4 days. Our parametrized model can be used to aid in design of clinical trials to study viral dynamics following analytic treatment interruption. We show how to derive informative personalized testing frequencies from our model and offer a proof-of-concept example. Our results represent first steps towards a model that can make predictions on a person living with HIV (PLWH)'s rebound time distribution based on biomarkers, and help identify PLWH with long viral rebound delays.

Hydrodynamics of transient cell-cell contact: The role of membrane permeability and active protrusion length

In many biological settings, two or more cells come into physical contact to form a cell-cell interface. In some cases, the cell-cell contact must be transient, forming on timescales of seconds. One example is offered by the T cell, an immune cell which must attach to the surface of other cells in order to decipher information about disease. The aspect ratio of these interfaces (tens of nanometers thick and tens of micrometers in diameter) puts them into the thin-layer limit, or "lubrication limit", of fluid dynamics. A key question is how the receptors and ligands on opposing cells come into contact. What are the relative roles of thermal undulations of the plasma membrane and deterministic forces from active filopodia? We use a computational fluid dynamics algorithm capable of simulating 10-nanometer-scale fluid-structure interactions with thermal fluctuations up to seconds- and microns-scales. We use this to simulate two opposing membranes, variously including thermal fluctuations, active forces, and membrane permeability. In some regimes dominated by thermal fluctuations, proximity is a rare event, which we capture by computing mean first-passage times using a Weighted Ensemble rare-event computational method. Our results demonstrate a parameter regime in which the time it takes for an active force to drive local contact actually increases if the cells are being held closer together (e.g., by nonspecific adhesion), a phenomenon we attribute to the thin-layer effect. This leads to an optimal initial cell-cell separation for fastest receptor-ligand binding, which could have relevance for the role of cellular protrusions like microvilli. We reproduce a previous experimental observation that fluctuation spatial scales are largely unaffected, but timescales are dramatically slowed, by the thin-layer effect. We also find that membrane permeability would need to be above physiological levels to abrogate the thin-layer effect.

Evolution of cooperation on large networks with community structure

Cooperation is a major factor in the evolution of human societies. The structure of social networks, which affects the dynamics of cooperation and other interpersonal phenomena, have common structural signatures. One of these signatures is the tendency to organize as groups. This tendency gives rise to networks with community structure, which are composed of distinct modules. In this paper, we study analytically the evolutionary game dynamics on large modular networks in the limit of weak selection. We obtain novel analytical conditions such that natural selection favours cooperation over defection. We calculate the transition point for each community to favour cooperation. We find that a critical inter-community link creation probability exists for given group density, such that the overall network supports cooperation even if individual communities inhibit it. As a byproduct, we present solutions for the critical benefit-to-cost ratio which perform with remarkable accuracy for diverse generative network models, including those with community structure and heavy-tailed degree distributions. We also demonstrate the generalizability of the results to arbitrary two-player games.

Dissipative particle dynamics simulation of critical pore size in a lipid bilayer membrane

We investigate with computer simulations the critical radius of pores in a lipid bilayer membrane. Ilton et al. (Ilton et al. 2016 Phys. Rev. Lett. 117, 257801 (doi:10.1103/PhysRevLett.117.257801)) recently showed that nucleated pores in a homopolymer film can increase or decrease in size, depending on whether they are larger or smaller than a critical size which scales linearly with film thickness. Using dissipative particle dynamics, a particle-based simulation method, we investigate the same scenario for a lipid bilayer membrane whose structure is determined by lipid-water interactions. We simulate a perforated membrane in which holes larger than a critical radius grow, while holes smaller than the critical radius close, as in the experiment of Ilton et al. (Ilton et al. 2016 Phys. Rev. Lett. 117, 257801 (doi:10.1103/PhysRevLett.117.257801)). By altering key system parameters such as the number of particles per lipid and the periodicity, we also describe scenarios in which pores of any initial size can seal or even remain stable, showing a fundamental difference in the behaviour of lipid membranes from polymer films.

Drug sensitivity prediction with high-dimensional mixture regression

This paper proposes a mixture regression model-based method for drug sensitivity prediction. The proposed method explicitly addresses two fundamental issues in drug sensitivity prediction, namely, population heterogeneity and feature selection pertaining to each of the subpopulations. The mixture regression model is estimated using the imputation-conditional consistency algorithm, and the resulting estimator is consistent. This paper also proposes an average-BIC criterion for determining the number of components for the mixture regression model. The proposed method is applied to the CCLE dataset, and the numerical results indicate that the proposed method can make a drastic improvement over the existing ones, such as random forest, support vector regression, and regularized linear regression, in both drug sensitivity prediction and feature selection. The p-values for the comparisons in drug sensitivity prediction can reach the order O(10-8) or lower for the drugs with heterogeneous populations.

Boundedness and Decay for the Teukolsky Equation on Kerr Spacetimes I: The Case | a | ≪ M

We prove boundedness and polynomial decay statements for solutions of the spin ± 2 Teukolsky equation on a Kerr exterior background with parameters satisfying | a | ≪ M . The bounds are obtained by introducing generalisations of the higher order quantities P and P _ used in our previous work on the linear stability of Schwarzschild. The existence of these quantities in the Schwarzschild case is related to the transformation theory of Chandrasekhar. In a followup paper, we shall extend this result to the general sub-extremal range of parameters | a | < M . As in the Schwarzschild case, these bounds provide the first step in proving the full linear stability of the Kerr metric to gravitational perturbations.

Top-down inputs drive neuronal network rewiring and context-enhanced sensory processing in olfaction

Much of the computational power of the mammalian brain arises from its extensive top-down projections. To enable neuron-specific information processing these projections have to be precisely targeted. How such a specific connectivity emerges and what functions it supports is still poorly understood. We addressed these questions in silico in the context of the profound structural plasticity of the olfactory system. At the core of this plasticity are the granule cells of the olfactory bulb, which integrate bottom-up sensory inputs and top-down inputs delivered by vast top-down projections from cortical and other brain areas. We developed a biophysically supported computational model for the rewiring of the top-down projections and the intra-bulbar network via adult neurogenesis. The model captures various previous physiological and behavioral observations and makes specific predictions for the cortico-bulbar network connectivity that is learned by odor exposure and environmental contexts. Specifically, it predicts that-after learning-the granule-cell receptive fields with respect to sensory and with respect to cortical inputs are highly correlated. This enables cortical cells that respond to a learned odor to enact disynaptic inhibitory control specifically of bulbar principal cells that respond to that odor. For this the reciprocal nature of the granule cell synapses with the principal cells is essential. Functionally, the model predicts context-enhanced stimulus discrimination in cluttered environments ('olfactory cocktail parties') and the ability of the system to adapt to its tasks by rapidly switching between different odor-processing modes. These predictions are experimentally testable. At the same time they provide guidance for future experiments aimed at unraveling the cortico-bulbar connectivity.

Invariance properties for the error function used for multilinear regression

The connections between the error function used in multilinear regression and the expected, or assumed, properties of the data are investigated. It is shown that two of the most basic properties often required in data analysis, scale and rotational invariance, are incompatible. With this, it is established that multilinear regression using an error function derived from a geometric mean is both scale and reflectively invariant. The resulting error function is also shown to have the property that its minimizer, under certain conditions, is well approximated using the centroid of the error simplex. It is then applied to several multidimensional real world data sets, and compared to other regression methods.

Probing the ecological and evolutionary history of a thermophilic cyanobacterial population via statistical properties of its microdiversity

Despite extensive DNA sequencing data derived from natural microbial communities, it remains a major challenge to identify the key evolutionary and ecological forces that shape microbial populations. We have focused on the extensive microdiversity of the cyanobacterium Synechococcus sp., which is a dominant member of the dense phototrophic biofilms in the hot springs of Yellowstone National Park. From deep amplicon sequencing of many loci and statistical analyses of these data, we showed previously that the population has undergone an unexpectedly high degree of homologous recombination, unlinking synonymous SNP-pair correlations even on intragenic length scales. Here, we analyze the genic amino acid diversity, which provides new evidence of selection and insights into the evolutionary history of the population. Surprisingly, some features of the data, including the spectrum of distances between genic-alleles, appear consistent with primarily asexual neutral drift. Yet the non-synonymous site frequency spectrum has too large an excess of low-frequency polymorphisms to result from negative selection on deleterious mutations given the distribution of coalescent times that we infer. And our previous analyses showed that the population is not asexual. Taken together, these apparently contradictory data suggest that selection, epistasis, and hitchhiking all play essential roles in generating and stabilizing the diversity. We discuss these as well as potential roles of ecological niches at genomic and genic levels. From quantitative properties of the diversity and comparative genomic data, we infer aspects of the history and inter-spring dispersal of the meta-population since it was established in the Yellowstone Caldera. Our investigations illustrate the need for combining multiple types of sequencing data and quantitative statistical analyses to develop an understanding of microdiversity in natural microbial populations.

VIPER: variability-preserving imputation for accurate gene expression recovery in single-cell RNA sequencing studies

We develop a method, VIPER, to impute the zero values in single-cell RNA sequencing studies to facilitate accurate transcriptome quantification at the single-cell level. VIPER is based on nonnegative sparse regression models and is capable of progressively inferring a sparse set of local neighborhood cells that are most predictive of the expression levels of the cell of interest for imputation. A key feature of our method is its ability to preserve gene expression variability across cells after imputation. We illustrate the advantages of our method through several well-designed real data-based analytical experiments.

Time-varying predatory behavior is primary predictor of fine-scale movement of wildland-urban cougars

BACKGROUND

While many species have suffered from the detrimental impacts of increasing human population growth, some species, such as cougars (Puma concolor), have been observed using human-modified landscapes. However, human-modified habitat can be a source of both increased risk and increased food availability, particularly for large carnivores. Assessing preferential use of the landscape is important for managing wildlife and can be particularly useful in transitional habitats, such as at the wildland-urban interface. Preferential use is often evaluated using resource selection functions (RSFs), which are focused on quantifying habitat preference using either a temporally static framework or researcher-defined temporal delineations. Many applications of RSFs do not incorporate time-varying landscape availability or temporally-varying behavior, which may mask conflict and avoidance behavior.

METHODS

Contemporary approaches to incorporate landscape availability into the assessment of habitat selection include spatio-temporal point process models, step selection functions, and continuous-time Markov chain (CTMC) models; in contrast with the other methods, the CTMC model allows for explicit inference on animal movement in continuous-time. We used a hierarchical version of the CTMC framework to model speed and directionality of fine-scale movement by a population of cougars inhabiting the Front Range of Colorado, U.S.A., an area exhibiting rapid population growth and increased recreational use, as a function of individual variation and time-varying responses to landscape covariates.

RESULTS

We found evidence for individual- and daily temporal-variability in cougar response to landscape characteristics. Distance to nearest kill site emerged as the most important driver of movement at a population-level. We also detected seasonal differences in average response to elevation, heat loading, and distance to roads. Motility was also a function of amount of development, with cougars moving faster in developed areas than in undeveloped areas.

CONCLUSIONS

The time-varying framework allowed us to detect temporal variability that would be masked in a generalized linear model, and improved the within-sample predictive ability of the model. The high degree of individual variation suggests that, if agencies want to minimize human-wildlife conflict management options should be varied and flexible. However, due to the effect of recursive behavior on cougar movement, likely related to the location and timing of potential kill-sites, kill-site identification tools may be useful for identifying areas of potential conflict.