Two-Locus Likelihoods under Variable Population Size and Fine-Scale Recombination Rate Estimation

Two-locus sampling probabilities have played a central role in devising an efficient composite likelihood method for estimating fine-scale recombination rates. Due to mathematical and computational challenges, these sampling probabilities are typically computed under the unrealistic assumption of a constant population size, and simulation studies have shown that resulting recombination rate estimates can be severely biased in certain cases of historical population size changes. To alleviate this problem, we develop here new methods to compute the sampling probability for variable population size functions that are piecewise constant. Our main theoretical result, implemented in a new software package called LDpop, is a novel formula for the sampling probability that can be evaluated by numerically exponentiating a large but sparse matrix. This formula can handle moderate sample sizes ($n \leq 50$) and demographic size histories with a large number of epochs ($\mathcal{D} \geq 64$). In addition, LDpop implements an approximate formula for the sampling probability that is reasonably accurate and scales to hundreds in sample size ($n \geq 256$). Finally, LDpop includes an importance sampler for the posterior distribution of two-locus genealogies, based on a new result for the optimal proposal distribution in the variable-size setting. Using our methods, we study how a sharp population bottleneck followed by rapid growth affects the correlation between partially linked sites. Then, through an extensive simulation study, we show that accounting for population size changes under such a demographic model leads to substantial improvements in fine-scale recombination rate estimation. LDpop is freely available for download at https://github.com/popgenmethods/ldpopComment: 32 pages, 13 figure

The Diversity and Abundance of the Benthic Macroinvertebrates in an Oligo-Mesotrophic Central Florida Lake

Benthic macroinvertebrates and physicochemical parameters were sampled monthly in Spring Lake, Florida from July, 1978, to June, 1979. Spring Lake is a slightly acidic, sand bottom lake located in the sandhill region of the Central Highlands. While submersed vegetation is not dense, the lake does contain an abundance of the endemic submersed plant Mayaca aubletii. The littoral zone is dominated by plants belonging to the genera Panicum, Nuphar, Hydrocotyle, Nymphaea, Satittaria, and Typha. The benthic macroinvertebrates collected consisted of 51 species; approximately 50 percent were in the family Chironomidae. The annual mean number of individuals was 947/m2. The mayfly Hexagenia munda Orlando was the most numerous species (18.4 percent of the annual mean); the Chironomidae was the most numerous family (31.6 percent of the annual mean). The annual mean value for the Simpson\u27s Index was 0.25 while the annual mean value for the Shannon Index was 2.60

A Likelihood-Free Inference Framework for Population Genetic Data using Exchangeable Neural Networks

An explosion of high-throughput DNA sequencing in the past decade has led to a surge of interest in population-scale inference with whole-genome data. Recent work in population genetics has centered on designing inference methods for relatively simple model classes, and few scalable general-purpose inference techniques exist for more realistic, complex models. To achieve this, two inferential challenges need to be addressed: (1) population data are exchangeable, calling for methods that efficiently exploit the symmetries of the data, and (2) computing likelihoods is intractable as it requires integrating over a set of correlated, extremely high-dimensional latent variables. These challenges are traditionally tackled by likelihood-free methods that use scientific simulators to generate datasets and reduce them to hand-designed, permutation-invariant summary statistics, often leading to inaccurate inference. In this work, we develop an exchangeable neural network that performs summary statistic-free, likelihood-free inference. Our framework can be applied in a black-box fashion across a variety of simulation-based tasks, both within and outside biology. We demonstrate the power of our approach on the recombination hotspot testing problem, outperforming the state-of-the-art.Comment: 9 pages, 8 figure

Inference of Population History using Coalescent HMMs: Review and Outlook

Studying how diverse human populations are related is of historical and anthropological interest, in addition to providing a realistic null model for testing for signatures of natural selection or disease associations. Furthermore, understanding the demographic histories of other species is playing an increasingly important role in conservation genetics. A number of statistical methods have been developed to infer population demographic histories using whole-genome sequence data, with recent advances focusing on allowing for more flexible modeling choices, scaling to larger data sets, and increasing statistical power. Here we review coalescent hidden Markov models, a powerful class of population genetic inference methods that can effectively utilize linkage disequilibrium information. We highlight recent advances, give advice for practitioners, point out potential pitfalls, and present possible future research directions.Comment: 12 pages, 2 figure

Sharp preasymptotic error bounds for the Helmholtz hh-FEM

In the analysis of the $h$-version of the finite-element method (FEM), with fixed polynomial degree $p$, applied to the Helmholtz equation with wavenumber $k\gg 1$, the $\textit{asymptotic regime}$ is when $(hk)^p C_{\rm sol}$ is sufficiently small and the sequence of Galerkin solutions are quasioptimal; here $C_{\rm sol}$ is the norm of the Helmholtz solution operator, normalised so that $C_{\rm sol} \sim k$ for nontrapping problems. In the $\textit{preasymptotic regime}$, one expects that if $(hk)^{2p}C_{\rm sol}$ is sufficiently small, then (for physical data) the relative error of the Galerkin solution is controllably small. In this paper, we prove the natural error bounds in the preasymptotic regime for the variable-coefficient Helmholtz equation in the exterior of a Dirichlet, or Neumann, or penetrable obstacle (or combinations of these) and with the radiation condition $\textit{either}$ realised exactly using the Dirichlet-to-Neumann map on the boundary of a ball $\textit{or}$ approximated either by a radial perfectly-matched layer (PML) or an impedance boundary condition. Previously, such bounds for $p>1$ were only available for Dirichlet obstacles with the radiation condition approximated by an impedance boundary condition. Our result is obtained via a novel generalisation of the "elliptic-projection" argument (the argument used to obtain the result for $p=1$) which can be applied to a wide variety of abstract Helmholtz-type problems

Modeling radiation belt radial diffusion in ULF wave fields: 1. Quantifying ULF wave power at geosynchronous orbit in observations and in global MHD model

[1] To provide critical ULF wave field information for radial diffusion studies in the radiation belts, we quantify ULF wave power (f = 0.5–8.3 mHz) in GOES observations and magnetic field predictions from a global magnetospheric model. A statistical study of 9 years of GOES data reveals the wave local time distribution and power at geosynchronous orbit in field-aligned coordinates as functions of wave frequency, solar wind conditions (Vx, ΔPd and IMF Bz) and geomagnetic activity levels (Kp, Dst and AE). ULF wave power grows monotonically with increasing solar wind Vx, dynamic pressure variations ΔPd and geomagnetic indices in a highly correlated way. During intervals of northward and southward IMF Bz, wave activity concentrates on the dayside and nightside sectors, respectively, due to different wave generation mechanisms in primarily open and closed magnetospheric configurations. Since global magnetospheric models have recently been used to trace particles in radiation belt studies, it is important to quantify the wave predictions of these models at frequencies relevant to electron dynamics (mHz range). Using 27 days of real interplanetary conditions as model inputs, we examine the ULF wave predictions modeled by the Lyon-Fedder-Mobarry magnetohydrodynamic code. The LFM code does well at reproducing, in a statistical sense, the ULF waves observed by GOES. This suggests that the LFM code is capable of modeling variability in the magnetosphere on ULF time scales during typical conditions. The code provides a long-missing wave field model needed to quantify the interaction of radiation belt electrons with realistic, global ULF waves throughout the inner magnetosphere

Eigenvalues of the truncated Helmholtz solution operator under strong trapping

For the Helmholtz equation posed in the exterior of a Dirichlet obstacle, we prove that if there exists a family of quasimodes (as is the case when the exterior of the obstacle has stable trapped rays), then there exist near-zero eigenvalues of the standard variational formulation of the exterior Dirichlet problem (recall that this formulation involves truncating the exterior domain and applying the exterior Dirichlet-to-Neumann map on the truncation boundary). Our motivation for proving this result is that a) the finite-element method for computing approximations to solutions of the Helmholtz equation is based on the standard variational formulation, and b) the location of eigenvalues, and especially near-zero ones, plays a key role in understanding how iterative solvers such as the generalised minimum residual method (GMRES) behave when used to solve linear systems, in particular those arising from the finite-element method. The result proved in this paper is thus the first step towards rigorously understanding how GMRES behaves when applied to discretisations of high-frequency Helmholtz problems under strong trapping (the subject of the companion paper [Marchand, Galkowski, Spence, Spence, 2021])

High-frequency estimates on boundary integral operators for the Helmholtz exterior Neumann problem

We study a commonly-used second-kind boundary-integral equation for solving the Helmholtz exterior Neumann problem at high frequency, where, writing $\Gamma$ for the boundary of the obstacle, the relevant integral operators map $L^2(\Gamma)$ to itself. We prove new frequency-explicit bounds on the norms of both the integral operator and its inverse. The bounds on the norm are valid for piecewise-smooth $\Gamma$ and are sharp up to factors of $\log k$ (where $k$ is the wavenumber), and the bounds on the norm of the inverse are valid for smooth $\Gamma$ and are observed to be sharp at least when $\Gamma$ is smooth with strictly-positive curvature. Together, these results give bounds on the condition number of the operator on $L^2(\Gamma)$; this is the first time $L^2(\Gamma)$ condition-number bounds have been proved for this operator for obstacles other than balls