Publications

One approach to analyzing the dynamics of a physical system is to search for long-lived patterns in its motions. This approach has been particularly successful for molecular dynamics data, where slowly decorrelating patterns can indicate large-scale conformational changes. Detecting such patterns is the central objective of the variational approach to conformational dynamics (VAC), as well as the related methods of time-lagged independent component analysis and Markov state modeling. In VAC, the search for slowly decorrelating patterns is formalized as a variational problem solved by the eigenfunctions of the system’s transition operator. VAC computes solutions to this variational problem by optimizing a linear or nonlinear model of the eigenfunctions using time series data. Here, we build on VAC’s success by addressing two practical limitations. First, VAC can give poor eigenfunction estimates when the lag time parameter is chosen poorly. Second, VAC can overfit when using flexible parametrizations such as artificial neural networks with insufficient regularization. To address these issues, we propose an extension that we call integrated VAC (IVAC). IVAC integrates over multiple lag times before solving the variational problem, making its results more robust and reproducible than VAC’s.

Multitaper estimators have enjoyed significant success in estimating spectral densities from finite samples using as tapers Slepian functions defined on the acquisition domain. Unfortunately, the numerical calculation of these Slepian tapers is only tractable for certain symmetric domains, such as rectangles or disks. In addition, no performance bounds are currently available for the mean squared error of the spectral density estimate. This situation is inadequate for applications such as cryo-electron microscopy, where noise models must be estimated from irregular domains with small sample sizes. We show that the multitaper estimator only depends on the linear space spanned by the tapers. As a result, Slepian tapers may be replaced by proxy tapers spanning the same subspace (validating the common practice of using partially converged solutions to the Slepian eigenproblem as tapers). These proxies may consequently be calculated using standard numerical algorithms for block diagonalization. We also prove a set of performance bounds for multitaper estimators on arbitrary domains. The method is demonstrated on synthetic and experimental datasets from cryo-electron microscopy, where it reduces mean squared error by a factor of two or more compared to traditional methods.

One of the most powerful approaches to imaging at the nanometer or subnanometer length scale is coherent diffraction imaging using X-ray sources. For amorphous (non-crystalline) samples, the raw data can be interpreted as the modulus of the continuous Fourier transform of the unknown object. Making use of prior information about the sample (such as its support), a natural goal is to recover the phase through computational means, after which the unknown object can be visualized at high resolution. While many algorithms have been proposed for this phase retrieval problem, careful analysis of its well-posedness has received relatively little attention. In this paper, we show that the problem is, in general, not well-posed and describe some of the underlying issues that are responsible for the ill-posedness. We then show how this analysis can be used to develop experimental protocols that lead to better conditioned inverse problems.

In the present paper we describe a class of algorithms for the solution of Laplace's equation on polygonal domains with Neumann boundary conditions. It is well known that in such cases the solutions have singularities near the corners which poses a challenge for many existing methods. If the boundary data is smooth on each edge of the polygon, then in the vicinity of each corner the solution to the corresponding boundary integral equation has an expansion in terms of certain (analytically available) singular powers. Using the known behavior of the solution, universal discretizations have been constructed for the solution of the Dirichlet problem. However, the leading order behavior of solutions to the Neumann problem is $O(t^{\mu})$ for $\mu \in (-1/2,0)$ depending on the angle at the corner (compared to $O(C+t^{\mu})$ with $\mu>1/2$ for the Dirichlet problem); this presents a significant challenge in the design of universal discretizations. Our approach is based on using the discretization for the Dirichlet problem in order to compute a solution in the "weak sense" by solving an adjoint linear system; namely, it can be used to compute inner products with smooth functions accurately, but it cannot be interpolated. Furthermore we present a procedure to obtain accurate solutions arbitrarily close to the corner, by solving a sequence of small local subproblems in the vicinity of that corner. The results are illustrated with several numerical examples.

When testing for a rare disease, prevalence estimates can be highly sensitive to uncertainty in the specificity and sensitivity of the test. Bayesian inference is a natural way to propagate these uncertainties, with hierarchical modeling capturing variation in these parameters across experiments. Another concern is the people in the sample not being representative of the general population. Statistical adjustment cannot without strong assumptions correct for selection bias in an opt-in sample, but multilevel regression and poststratification can at least adjust for known differences between sample and population. We demonstrate these models with code in R and Stan and discuss their application to a controversial recent study of COVID-19 antibodies in a sample of people from the Stanford University area. Wide posterior intervals make it impossible to evaluate the quantitative claims of that study regarding the number of unreported infections. For future studies, the methods described here should facilitate more accurate estimates of disease prevalence from imperfect tests performed on non-representative samples

The Quijote simulations are a set of 44,100 full N-body simulations spanning more than 7000 cosmological models in the $\{{{\rm{\Omega }}}_{{\rm{m}}},{{\rm{\Omega }}}_{{\rm{b}}},h,{n}_{s},{\sigma }_{8},{M}_{\nu },w\}$ hyperplane. At a single redshift, the simulations contain more than 8.5 trillion particles over a combined volume of 44,100 ${\left({h}^{-1}\mathrm{Gpc}\right)}^{3};$ each simulation follows the evolution of \(256^{3}, 512^{3}, or 1024^{3}\) particles in a box of 1 h −1 Gpc length. Billions of dark matter halos and cosmic voids have been identified in the simulations, whose runs required more than 35 million core hours. The Quijote simulations have been designed for two main purposes: (1) to quantify the information content on cosmological observables and (2) to provide enough data to train machine-learning algorithms. In this paper, we describe the simulations and show a few of their applications. We also release the petabyte of data generated, comprising hundreds of thousands of simulation snapshots at multiple redshifts; halo and void catalogs; and millions of summary statistics, such as power spectra, bispectra, correlation functions, marked power spectra, and estimated probability density functions

Dorsal root ganglia (DRG) are promising sites for recording sensory activity. Current technologies for DRG recording are stiff and typically do not have sufficient site density for high-density neural data techniques. We demonstrate neural recordings in feline sacral DRG using a flexible polyimide microelectrode array with 30-40 μm site spacing. We delivered arrays into DRG with ultrananocrystalline diamond shuttles designed for high stiffness with small footprint. We recorded neural activity during sensory activation, including cutaneous brushing and bladder filling. We successfully delivered arrays in 5/6 experiments and recorded sensory activity in 4. Median signal amplitude was 55 μV and the maximum unique units recorded at one array position was 260, with 157 driven by sensory or electrical stimulation. We used specialized high-density neural signal analysis software to sort neural signals and, in one experiment, track 8 signals as the array was retracted 500 μm. This study is the first demonstration of ultrathin, flexible, high-density electronics delivered into DRG, with capabilities for recording and tracking sensory information that are a significant improvement over conventional DRG interfaces.Competing Interest StatementT.M.B. is a named inventors on a granted patent (US9622671B2; assigned to University of Pittsburgh) which is on the monitoring of physiological states via microelectrodes at DRG. The authors declare no other personal or institutional interest with regards to the authorship and/or publication of this manuscript.

The method of Helsing and co-workers evaluates Laplace and related layer potentials generated by a panel (composite) quadrature on a curve, efficiently and with high-order accuracy for arbitrarily close targets. Since it exploits complex analysis, its use has been restricted to two dimensions (2D). We first explain its loss of accuracy as panels become curved, using a classical complex approximation result of Walsh that can be interpreted as "electrostatic shielding" of a Schwarz singularity. We then introduce a variant that swaps the target singularity for one at its complexified parameter preimage; in the latter space the panel is flat, hence the convergence rate can be much higher. The preimage is found robustly by Newton iteration. This idea also enables, for the first time, a near-singular quadrature for potentials generated by smooth curves in 3D, building on recurrences of Tornberg-Gustavsson. We apply this to accurate evaluation of the Stokes flow near to a curved filament in the slender body approximation. Our 3D method is several times more efficient (both in terms of kernel evaluations, and in speed in a C implementation) than the only existing alternative, namely, adaptive integration.

A novel optimization formulation is provided for holographic phase retrieval with data subjected to Poisson shot noise in the low-photon regime. This formulation, based on regularized maximum likelihood estimation, consistently provides improved image reconstruction.

We present a neural network architecture that is fully equivariant with respect to transformations under the Lorentz group, a fundamental symmetry of space and time in physics. The architecture is based on the theory of the finite-dimensional representations of the Lorentz group and the equivariant nonlinearity involves the tensor product. For classification tasks in particle physics, we show that such an equivariant architecture leads to drastically simpler models that have relatively few learnable parameters and are much more physically interpretable than leading approaches that use CNNs and point cloud approaches. The performance of the network is tested on a public classification dataset [https://zenodo.org/record/2603256] for tagging top quark decays given energy-momenta of jet constituents produced in proton-proton collisions.