Stefan Engblom

Stefan Engblom, Lars Ferm, Andreas Hellander et al.

Stochastic chemical systems with diffusion are modeled with a reaction-diffusion master equation. On a macroscopic level, the governing equation is a reaction-diffusion equation for the averages of the chemical species. On a mesoscopic level, the master equation for a well stirred chemical system is combined with Brownian motion in space to obtain the reaction-diffusion master equation. The space is covered by an unstructured mesh and the diffusion coefficients on the mesoscale are obtained from a finite element discretization of the Laplace operator on the macroscale. The resulting method is a flexible hybrid algorithm in that the diffusion can be handled either on the meso- or on the macroscale level. The accuracy and the efficiency of the method are illustrated in three numerical examples inspired by molecular biology.

Stefan Engblom, Minh Do-Quang, Gustav Amberg et al.

An existing phase-field model of two immiscible fluids with a single soluble surfactant present is discussed in detail. We analyze the well-posedness of the model and provide strong evidence that it is mathematically ill-posed for a large set of physically relevant parameters. As a consequence, critical modifications to the model are suggested that substantially increase the domain of validity. Carefully designed numerical simulations offer informative demonstrations as to the sharpness of our theoretical results and the qualities of the physical model. A fully coupled hydrodynamic test-case demonstrates the potential to capture also non-trivial effects on the overall flow.

Stefan Engblom

A version of the time-parallel algorithm parareal is analyzed and applied to stochastic models in chemical kinetics. A fast predictor at the macroscopic scale (evaluated in serial) is available in the form of the usual reaction rate equations. A stochastic simulation algorithm is used to obtain an exact realization of the process at the mesoscopic scale (in parallel). The underlying stochastic description is a jump process driven by the Poisson measure. A convergence result in this arguably difficult setting is established suggesting that a homogenization of the solution is advantageous. We devise a simple but highly general such technique. Three numerical experiments on models representative to the field of computational systems biology illustrate the method. For non-stiff problems, it is shown that the method is able to quickly converge even when stochastic effects are present. For stiff problems we are instead able to obtain fast convergence to a homogenized solution. Overall, the method builds an attractive bridge between on the one hand, macroscopic deterministic scales and, on the other hand, mesoscopic stochastic ones. This construction is clearly possible to apply also to stochastic models within other fields.

Lina Meinecke, Stefan Engblom, Andreas Hellander et al.

In computational system biology, the mesoscopic model of reaction-diffusion kinetics is described by a continuous time, discrete space Markov process. To simulate diffusion stochastically, the jump coefficients are obtained by a discretization of the diffusion equation. Using unstructured meshes to represent complicated geometries may lead to negative coefficients when using piecewise linear finite elements. Several methods have been proposed to modify the coefficients to enforce the non-negativity needed in the stochastic setting. In this paper, we present a method to quantify the error introduced by that change. We interpret the modified discretization matrix as the exact finite element discretization of a perturbed equation. The forward error, the error between the analytical solutions to the original and the perturbed equations, is bounded by the backward error, the error between the diffusion of the two equations. We present a backward analysis algorithm to compute the diffusion coefficient from a given discretization matrix. The analysis suggests a new way of deriving non-negative jump coefficients that minimizes the backward error. The theory is tested in numerical experiments indicating that the new method is superior and minimizes also the forward error.

Stefan Engblom

Motivated by the lack of a suitable constructive framework for analyzing popular stochastic models of Systems Biology, we devise conditions for existence and uniqueness of solutions to certain jump stochastic differential equations (SDEs). Working from simple examples we find reasonable and explicit assumptions on the driving coefficients for the SDE representation to make sense. By `reasonable' we mean that stronger assumptions generally do not hold for systems of practical interest. In particular, we argue against the traditional use of global Lipschitz conditions and certain common growth restrictions. By `explicit', finally, we like to highlight the fact that the various constants occurring among our assumptions all can be determined once the model is fixed. We show how basic long time estimates and some limit results for perturbations can be derived in this setting such that these can be contrasted with the corresponding estimates from deterministic dynamics. The main complication is that the natural path-wise representation is generated by a counting measure with an intensity that depends nonlinearly on the state.

Doghonay Arjmand, Stefan Engblom, Gunilla Kreiss

We consider a multiscale strategy addressing the disparate scales in the Landau-Lifschitz equations in micro-magnetism. At the microscopic scale, the dynamics of magnetic moments are driven by a high frequency field. On the macroscopic scale we are interested in simulating the dynamics of the magnetisation without fully resolving the microscopic scales. The method follows the framework of heterogeneous multiscale methods and it has two main ingredients: a micro- and a macroscale model. The microscopic model is assumed to be known exactly whereas the macro model is incomplete as it lacks effective quantities. The two models use different temporal and spatial scales and effective parameter values for the macro model are computed on the fly, allowing for improved efficiency over traditional one-scale schemes. For the analysis, we consider a single spin under a high frequency field and show that effective quantities can be obtained accurately with step-sizes much larger than the size of the microscopic scales required to resolve the microscopic features. Numerical results both for a single magnetic particle as well as a chain of interacting magnetic particles are given to validate the theory.

Stefan Engblom, Per Lötstedt, Lina Meinecke

We develop a mesoscopic modeling framework for diffusion in a crowded environment, particularly targeting applications in the modeling of living cells. Through homogenization techniques we effectively coarse-grain a detailed microscopic description into a previously developed internal state diffusive framework. The observables in the mesoscopic model correspond to solutions of macroscopic partial differential equations driven by stochastically varying diffusion fields in space and time. Analytical solutions and numerical experiments illustrate the framework.

Anders Goude, Stefan Engblom

We develop an efficient and high order panel method with applications in airfoil design. Through the use of analytic work and careful considerations near singularities our approach is quadrature-free. The resulting method is examined with respect to accuracy and efficiency and we discuss the different trade-offs in approximation order and computational complexity. A reference implementation within a package for a two-dimensional fast multipole method is distributed freely.

Fredrik Wrede, Robin Eriksson, Richard Jiang et al.

State-of-the-art neural network-based methods for learning summary statistics have delivered promising results for simulation-based likelihood-free parameter inference. Existing approaches require density estimation as a post-processing step building upon deterministic neural networks, and do not take network prediction uncertainty into account. This work proposes a robust integrated approach that learns summary statistics using Bayesian neural networks, and directly estimates the posterior density using categorical distributions. An adaptive sampling scheme selects simulation locations to efficiently and iteratively refine the predictive posterior of the network conditioned on observations. This allows for more efficient and robust convergence on comparatively large prior spaces. We demonstrate our approach on benchmark examples and compare against related methods.

Jing Liu, Stefan Engblom, Carl Nettelblad

Modern Flash X-ray diffraction Imaging (FXI) acquires diffraction signals from single biomolecules at a high repetition rate from X-ray Free Electron Lasers (XFELs), easily obtaining millions of 2D diffraction patterns from a single experiment. Due to the stochastic nature of FXI experiments and the massive volumes of data, retrieving 3D electron densities from raw 2D diffraction patterns is a challenging and time-consuming task. We propose a semi-automatic data analysis pipeline for FXI experiments, which includes four steps: hit finding and preliminary filtering, pattern classification, 3D Fourier reconstruction, and post analysis. We also include a recently developed bootstrap methodology in the post-analysis step for uncertainty analysis and quality control. To achieve the best possible resolution, we further suggest using background subtraction, signal windowing, and convex optimization techniques when retrieving the Fourier phases in the post-analysis step. As an application example, we quantified the 3D electron structure of the PR772 virus using the proposed data-analysis pipeline. The retrieved structure was above the detector-edge resolution and clearly showed the pseudo-icosahedral capsid of the PR772.

Jing Liu, Gijs van der Schot, Stefan Engblom

Current Flash X-ray single-particle diffraction Imaging (FXI) experiments, which operate on modern X-ray Free Electron Lasers (XFELs), can record millions of interpretable diffraction patterns from individual biomolecules per day. Due to the stochastic nature of the XFELs, those patterns will to a varying degree include scatterings from contaminated samples. Also, the heterogeneity of the sample biomolecules is unavoidable and complicates data processing. Reducing the data volumes and selecting high-quality single-molecule patterns are therefore critical steps in the experimental set-up. In this paper, we present two supervised template-based learning methods for classifying FXI patterns. Our Eigen-Image and Log-Likelihood classifier can find the best-matched template for a single-molecule pattern within a few milliseconds. It is also straightforward to parallelize them so as to fully match the XFEL repetition rate, thereby enabling processing at site.

Augustin Chevallier, Stefan Engblom

Stochastic computational models in the form of pure jump processes occur frequently in the description of chemical reactive processes, of ion channel dynamics, and of the spread of infections in populations. For spatially extended models, the computational complexity can be rather high such that approximate multiscale models are attractive alternatives. Within this framework some variables are described stochastically, while others are approximated with a macroscopic point value. We devise theoretical tools for analyzing the pathwise multiscale convergence of this type of variable splitting methods, aiming specifically at spatially extended models. Notably, the conditions we develop guarantee well-posedness of the approximations without requiring explicit assumptions of \textit{a priori} bounded solutions. We are also able to quantify the effect of the different sources of errors, namely the \emph{multiscale error} and the \emph{splitting error}, respectively, by developing suitable error bounds. Computational experiments on selected problems serve to illustrate our findings.

Stefan Engblom, Carl Nettelblad, Jing Liu

Modern technology for producing extremely bright and coherent X-ray laser pulses provides the possibility to acquire a large number of diffraction patterns from individual biological nanoparticles, including proteins, viruses, and DNA. These two-dimensional diffraction patterns can be practically reconstructed and retrieved down to a resolution of a few \angstrom. In principle, a sufficiently large collection of diffraction patterns will contain the required information for a full three-dimensional reconstruction of the biomolecule. The computational methodology for this reconstruction task is still under development and highly resolved reconstructions have not yet been produced. We analyze the Expansion-Maximization-Compression scheme, the current state of the art approach for this very challenging application, by isolating different sources of uncertainty. Through numerical experiments on synthetic data we evaluate their respective impact. We reach conclusions of relevance for handling actual experimental data, as well as pointing out certain improvements to the underlying estimation algorithm. We also introduce a practically applicable computational methodology in the form of bootstrap procedures for assessing reconstruction uncertainty in the real data case. We evaluate the sharpness of this approach and argue that this type of procedure will be critical in the near future when handling the increasing amount of data.

Stefan Engblom, Dimitar Lukarski

We develop and implement in this paper a fast sparse assembly algorithm, the fundamental operation which creates a compressed matrix from raw index data. Since it is often a quite demanding and sometimes critical operation, it is of interest to design a highly efficient implementation. We show how to do this, and moreover, we show how our implementation can be parallelized to utilize the power of modern multicore computers. Our freely available code, fully Matlab compatible, achieves about a factor of 5 times in speedup on a typical 6-core machine and 10 times on a dual-socket 16 core machine compared to the built-in serial implementation.

Stefan Engblom

Mesoscopic models in the reaction-diffusion framework have gained recognition as a viable approach to describing chemical processes in cell biology. The resulting computational problem is a continuous-time Markov chain on a discrete and typically very large state space. Due to the many temporal and spatial scales involved many different types of computationally more effective multiscale models have been proposed, typically coupling different types of descriptions within the Markov chain framework. In this work we look at the strong convergence properties of the basic first order Strang, or Lie-Trotter, split-step method, which is formed by decoupling the dynamics in finite time-steps. Thanks to its simplicity and flexibility, this approach has been tried in many different combinations. We develop explicit sufficient conditions for path-wise well-posedness and convergence of the method, including error estimates, and we illustrate our findings with numerical examples. In doing so, we also suggest a certain partition of unity representation for the split-step method, which in turn implies a concrete simulation algorithm under which trajectories may be compared in a path-wise sense.

Emilie Blanc, Stefan Engblom, Andreas Hellander et al.

Subdiffusion has been proposed as an explanation of various kinetic phenomena inside living cells. In order to fascilitate large-scale computational studies of subdiffusive chemical processes, we extend a recently suggested mesoscopic model of subdiffusion into an accurate and consistent reaction-subdiffusion computational framework. Two different possible models of chemical reaction are revealed and some basic dynamic properties are derived. In certain cases those mesoscopic models have a direct interpretation at the macroscopic level as fractional partial differential equations in a bounded time interval. Through analysis and numerical experiments we estimate the macroscopic effects of reactions under subdiffusive mixing. The models display properties observed also in experiments: for a short time interval the behavior of the diffusion and the reaction is ordinary, in an intermediate interval the behavior is anomalous, and at long times the behavior is ordinary again.

Tomas Ekeberg, Stefan Engblom, Jing Liu

The classical method of determining the atomic structure of complex molecules by analyzing diffraction patterns is currently undergoing drastic developments. Modern techniques for producing extremely bright and coherent X-ray lasers allow a beam of streaming particles to be intercepted and hit by an ultrashort high energy X-ray beam. Through machine learning methods the data thus collected can be transformed into a three-dimensional volumetric intensity map of the particle itself. The computational complexity associated with this problem is very high such that clusters of data parallel accelerators are required. We have implemented a distributed and highly efficient algorithm for inversion of large collections of diffraction patterns targeting clusters of hundreds of GPUs. With the expected enormous amount of diffraction data to be produced in the foreseeable future, this is the required scale to approach real time processing of data at the beam site. Using both real and synthetic data we look at the scaling properties of the application and discuss the overall computational viability of this exciting and novel imaging technique.