Search

Reproducibility is a foundational standard for validating scientific claims in computational research. Stochastic computational models are employed across diverse fields such as systems biology, financial modelling and environmental sciences. Existing infrastructure and software tools support various aspects of reproducible model development, application, and dissemination, but do not adequately address independently reproducing simulation results that form the basis of scientific conclusions. To bridge this gap, we introduce the Empirical Characteristic Function Equality Convergence Test (EFECT), a data-driven method to quantify the reproducibility of stochastic simulation results. EFECT employs empirical characteristic functions to compare reported results with those independently generated by assessing distributional inequality, termed EFECT error, a metric to quantify the likelihood of equality. Additionally, we establish the EFECT convergence point, a metric for determining the required number of simulation runs to achieve an EFECT error value of a priori statistical significance, setting a reproducibility benchmark. EFECT supports all real-valued and bounded results irrespective of the model or method that produced them, and accommodates stochasticity from intrinsic model variability and random sampling of model inputs. We tested EFECT with stochastic differential equations, agent-based models, and Boolean networks, demonstrating its broad applicability and effectiveness. EFECT standardizes stochastic simulation reproducibility, establishing a workflow that guarantees reliable results, supporting a wide range of stakeholders, and thereby enhancing validation of stochastic simulation studies, across a model's lifecycle. To promote future standardization efforts, we are developing open source software library libSSR in diverse programming languages for easy integration of EFECT., Comment: 25 pages, 4 figures

Medical digital twins are computational models of human biology relevant to a given medical condition, which can be tailored to an individual patient, thereby predicting the course of disease and individualized treatments, an important goal of personalized medicine. The immune system, which has a central role in many diseases, is highly heterogeneous between individuals, and thus poses a major challenge for this technology. If medical digital twins are to faithfully capture the characteristics of a patient's immune system, we need to answer many questions, such as: What do we need to know about the immune system to build mathematical models that reflect features of an individual? What data do we need to collect across the different scales of immune system action? What are the right modeling paradigms to properly capture immune system complexity? In February 2023, an international group of experts convened in Lake Nona, FL for two days to discuss these and other questions related to digital twins of the immune system. The group consisted of clinicians, immunologists, biologists, and mathematical modelers, representative of the interdisciplinary nature of medical digital twin development. A video recording of the entire event is available. This paper presents a synopsis of the discussions, brief descriptions of ongoing digital twin projects at different stages of progress. It also proposes a 5-year action plan for further developing this technology. The main recommendations are to identify and pursue a small number of promising use cases, to develop stimulation-specific assays of immune function in a clinical setting, and to develop a database of existing computational immune models, as well as advanced modeling technology and infrastructure.

Neural networks (NNs) have proven to be a viable alternative to traditional direct numerical algorithms, with the potential to accelerate computational time by several orders of magnitude. In the present paper we study the use of encoder-decoder convolutional neural network (CNN) as surrogates for steady-state diffusion solvers. The construction of such surrogates requires the selection of an appropriate task, network architecture, training set structure and size, loss function, and training algorithm hyperparameters. It is well known that each of these factors can have a significant impact on the performance of the resultant model. Our approach employs an encoder-decoder CNN architecture, which we posit is particularly well-suited for this task due to its ability to effectively transform data, as opposed to merely compressing it. We systematically evaluate a range of loss functions, hyperparameters, and training set sizes. Our results indicate that increasing the size of the training set has a substantial effect on reducing performance fluctuations and overall error. Additionally, we observe that the performance of the model exhibits a logarithmic dependence on the training set size. Furthermore, we investigate the effect on model performance by using different subsets of data with varying features. Our results highlight the importance of sampling the configurational space in an optimal manner, as this can have a significant impact on the performance of the model and the required training time. In conclusion, our results suggest that training a model with a pre-determined error performance bound is not a viable approach, as it does not guarantee that edge cases with errors larger than the bound do not exist. Furthermore, as most surrogate tasks involve a high dimensional landscape, an ever increasing training set size is, in principle, needed, however it is not a practical solution., Comment: 21 ps, 17 figs, 8 ts

Modeling and simulation (M&S) has revolutionized the way modern engineered products are designed, tested and evaluated. Yet modeling and simulation is much less frequently applied in the study of natural biological and active matter systems. Two of the greatest challenges in applying M\&S approaches to natural biological systems are (1) difficulty in specifying a model, and developing a simulation from it, and (2) tuning and optimizing parameters. Here we address the first challenge directly, and develop a software library that can lead to progress in the second.

In many mechanistic medical, biological, physical and engineered spatiotemporal dynamic models the numerical solution of partial differential equations (PDEs) can make simulations impractically slow. Biological models require the simultaneous calculation of the spatial variation of concentration of dozens of diffusing chemical species. Machine learning surrogates, neural networks trained to provide approximate solutions to such complicated numerical problems, can often provide speed-ups of several orders of magnitude compared to direct calculation. PDE surrogates enable use of larger models than are possible with direct calculation and can make including such simulations in real-time or near-real time workflows practical. Creating a surrogate requires running the direct calculation tens of thousands of times to generate training data and then training the neural network, both of which are computationally expensive. We use a Convolutional Neural Network to approximate the stationary solution to the diffusion equation in the case of two equal-diameter, circular, constant-value sources located at random positions in a two-dimensional square domain with absorbing boundary conditions. To improve convergence during training, we apply a training approach that uses roll-back to reject stochastic changes to the network that increase the loss function. The trained neural network approximation is about 1e3 times faster than the direct calculation for individual replicas. Because different applications will have different criteria for acceptable approximation accuracy, we discuss a variety of loss functions and accuracy estimators that can help select the best network for a particular application., Comment: 17 pp, 2 tables, 11 figs, 1sm, 12sm-figs, code available at GitHub

Generative models rely on the key idea that data can be represented in terms of latent variables which are uncorrelated by definition. Lack of correlation is important because it suggests that the latent space manifold is simpler to understand and manipulate. Generative models are widely used in deep learning, e.g., variational autoencoders (VAEs) and generative adversarial networks (GANs). Here we propose a method to build a set of linearly independent vectors in the latent space of a GANs, which we call quasi-eigenvectors. These quasi-eigenvectors have two key properties: i) They span all the latent space, ii) A set of these quasi-eigenvectors map to each of the labeled features one-on-one. We show that in the case of the MNIST, while the number of dimensions in latent space is large by construction, 98% of the data in real space map to a sub-domain of latent space of dimensionality equal to the number of labels. We then show how the quasi-eigenvalues can be used for Latent Spectral Decomposition (LSD), which has applications in denoising images and for performing matrix operations in latent space that map to feature transformations in real space. We show how this method provides insight into the latent space topology. The key point is that the set of quasi-eigenvectors form a basis set in latent space and each direction corresponds to a feature in real space., Comment: 9pp, 8 figs, 1 pseudocode, code available

Active Matter models commonly consider particles with overdamped dynamics subject to a force (speed) with constant modulus and random direction. Some models include also random noise in particle displacement (Wiener process) resulting in a diffusive motion at short time scales. On the other hand, Ornstein-Uhlenbeck processes consider Langevin dynamics for the particle velocity and predict a motion that is not diffusive at short time scales. However, experiments show that migrating cells may present a varying speed as well as a short-time diffusive behavior. While Ornstein-Uhlenbeck processes can describe the varying speed, Active Mater models can explain the short-time diffusive behavior. Isotropic models cannot explain both: short-time diffusion renders instantaneous velocity ill-defined, hence impeding dynamical equations that consider velocity time-derivatives. On the other hand, both models apply for migrating biological cells and must, in some limit, yield the same observable predictions. Here we propose and analytically solve an Anisotropic Ornstein-Uhlenbeck process that considers polarized particles, with a Langevin dynamics for the particle movement in the polarization direction while following a Wiener process for displacement in the orthogonal direction. Our characterization provides a theoretically robust way to compare movement in dimensionless simulations to movement in dimensionful experiments, besides proposing a procedure to deal with inevitable finite precision effects in experiments or simulations.

The convergence of HPC and data-intensive methodologies provide a promising approach to major performance improvements. This paper provides a general description of the interaction between traditional HPC and ML approaches and motivates the Learning Everywhere paradigm for HPC. We introduce the concept of effective performance that one can achieve by combining learning methodologies with simulation-based approaches, and distinguish between traditional performance as measured by benchmark scores. To support the promise of integrating HPC and learning methods, this paper examines specific examples and opportunities across a series of domains. It concludes with a series of open computer science and cyberinfrastructure questions and challenges that the Learning Everywhere paradigm presents.

Cell crawling is critical to biological development, homeostasis and disease. In many cases, cell trajectories are quasi-random-walk. In vitro assays on flat surfaces often described such quasi-random-walk cell trajectories as approximations to a solution of a Langevin process. However, experiments show quasi-diffusive behavior at small timescales, indicating that instantaneous velocity and velocity autocorrelations are not well-defined. We propose to characterize mean-squared cell displacement using a modified F\"urth equation with three temporal and spatial regimes: short- and long-time/range diffusion and intermediate time/range ballistic motion. This analysis collapses mean-squared displacements of previously published experimental data onto a single-parameter family of curves, allowing direct comparison between movement in different cell types, and between experiments and numerical simulations. Our method also show that robust cell-motility quantification requires an experiment with a maximum interval between images of a few percent of the cell-motion persistence time or less, and a duration of a few orders-of-magnitude longer than the cell-motion persistence time or more., Comment: 5 pages, plus Supplemental material

Biological cells are the prototypical example of active matter. Cells sense and respond to mechanical, chemical and electrical environmental stimuli with a range of behaviors, including dynamic changes in morphology and mechanical properties, chemical uptake and secretion, cell differentiation, proliferation, death, or migration. Modeling and simulation of such dynamic phenomena poses a number of computational challenges. A modeling language to describe cellular dynamics must be able to naturally represent complex intra and extra-cellular spatial structures, and coupled mechanical, chemical and electrical processes. In order to be useful to domain experts, a modeling language should be based on concepts, terms and principles native to the problem domain. A compiler must then be able to generate an executable model from this physically motivated description. Finally, an executable model must efficiently calculate the time evolution of such dynamic and inhomogeneous phenomena. We present a spatial hybrid systems modeling language, compiler and mesh-free Lagrangian based simulation engine which will enable domain experts to define models using natural, biologically motivated constructs and to simulate time evolution of coupled cellular, mechanical and chemical processes acting on a time varying number of cells and their environment., Comment: Symp. Theory of Modeling & Simulation (TMS/DEVS) Spring Sim 2017, Virginia Beach VA, USA, 2017

This paper presents libRoadRunner, an extensible, high-performance, cross-platform, open-source software library for the simulation and analysis of models \ expressed using Systems Biology Markup Language (SBML). SBML is the most widely used standard for representing dynamic networks, especially biochemical networks. libRoadRunner supports solution of both large models and multiple replicas of a single model on desktop, mobile and cluster computers. libRoadRunner is a self-contained library, able to run both as a component inside other tools via its C++ and C bindings andnteractively through its Python interface. The Python Application Programming Interface (API) is similar to the APIs of Matlab and SciPy, making it fast and easy to learn, even for new users. libRoadRunner uses a custom Just-In-Time (JIT) compiler built on the widely-used LLVM JIT compiler framework to compile SBML-specified models directly into very fast native machine code for a variety of processors, making it appropriate for solving very large models or multiple replicas of smaller models. libRoadRunner is flexible, supporting the bulk of the SBML specification (except for delay and nonlinear algebraic equations) and several of its extensions. It offers multiple deterministic and stochastic integrators, as well as tools for steady-state, stability analyses and flux balance analysis. We regularly update libRoadRunner binary distributions for Mac OS X, Linux and Windows and license them under Apache License Version 2.0. http://www.libroadrunner.org provides online documentation, full build instructions, binaries and a git source repository.

In wet liquid foams, slow diffusion of gas through bubble walls changes bubble pressure, volume and wall curvature. Large bubbles grow at the expenses of smaller ones. The smaller the bubble, the faster it shrinks. As the number of bubbles in a given volume decreases in time, the average bubble size increases: i.e. the foam coarsens. During coarsening, bubbles also move relative to each other, changing bubble topology and shape, while liquid moves within the regions separating the bubbles. Analyzing the combined effects of these mechanisms requires examining a volume with enough bubbles to provide appropriate statistics throughout coarsening. Using a Cellular Potts model, we simulate these mechanisms during the evolution of three-dimensional foams with wetnesses of $\phi=0.00$, $0.05$ and $ 0.20$. We represent the liquid phase as an ensemble of many small fluid particles, which allows us to monitor liquid flow in the region between bubbles. The simulations begin with $2 \times 10^5$ bubbles for $\phi = 0.00$ and $1.25 \times 10^5$ bubbles for $\phi = 0.05$ and $0.20$, allowing us to track the distribution functions for bubble size, topology and growth rate over two and a half decades of volume change. All simulations eventually reach a self-similar growth regime, with the distribution functions time independent and the number of bubbles decreasing with time as a power law whose exponent depends on the wetness.

We used X-ray tomography to characterize the geometry of all bubbles in a liquid foam of average liquid fraction $\phi_l\approx 17 %$ and to follow their evolution, measuring the normalized growth rate $\mathcal{G}=V^{-{1/3}}\frac{dV} {dt}$ for 7000 bubbles. While $\mathcal{G}$ does not depend only on the number of faces of a bubble, its average over $f-$faced bubbles scales as $G_f\sim f-f_0$ for large $f$s at all times. We discuss the dispersion of $\mathcal{G}$ and the influence of $V$ on $\mathcal{G}$., Comment: 10 pages, submitted to PRL

We have generated a novel form of shear banding in a 2D foam and measured the relative magnitude of drag forces on soap films at different lubrication layers. We injected air part way along a flowing bubble field in a narrow Hele-Shaw cell. The injected air inflates bubbles as they flow by, and these bubbles form a shear band down the middle of the Hele-Shaw cell. This channel appears to select a height that minimizes the total dissipation. Fitting data to a simple theoretical model, we show that the drag force on a soap film in contact with the wetting layer on a plate of glass is two orders of magnitude larger than the drag on a soap film in contact with another free soap film., Comment: Revised for submission, 05/15/07

This paper has been withdrawn, Comment: This paper has been withdrawn

The Cellular Potts Model (CPM) succesfully simulates drainage and shear in foams. Here we use the CPM to investigate instabilities due to the flow of a single large bubble in a dry, monodisperse two-dimensional flowing foam. As in experiments in a Hele-Shaw cell, above a threshold velocity the large bubble moves faster than the mean flow. Our simulations reproduce analytical and experimental predictions for the velocity threshold and the relative velocity of the large bubble, demonstrating the utility of the CPM in foam rheology studies., Comment: 10 pages, 3 figures. Replaced with revised version accepted for publication in JSTAT

Blood vessels form either when dispersed endothelial cells (the cells lining the inner walls of fully-formed blood vessels) organize into a vessel network (vasculogenesis), or by sprouting or splitting of existing blood vessels (angiogenesis). Although they are closely related biologically, no current model explains both phenomena with a single biophysical mechanism. Most computational models describe sprouting at the level of the blood vessel, ignoring how cell behavior drives branch splitting during sprouting. We present a cell-based, Glazier-Graner-Hogeweg-model simulation of the initial patterning before the vascular cords form lumens, based on plausible behaviors of endothelial cells. The endothelial cells secrete a chemoattractant, which attracts other endothelial cells. As in the classic Keller-Segel model, chemotaxis by itself causes cells to aggregate into isolated clusters. However, including experimentally-observed adhesion-driven contact inhibition of chemotaxis in the simulation causes randomly-distributed cells to organize into networks and cell aggregates to sprout, reproducing aspects of both de novo and sprouting blood-vessel growth. We discuss two branching instabilities responsible for our results. Cells at the surfaces of cell clusters attempting to migrate to the centers of the clusters produce a buckling instability. In a model variant that eliminates the surface-normal force, a dissipative mechanism drives sprouting, with the secreted chemical acting both as a chemoattractant and as an inhibitor of pseudopod extension. The branching instabilities responsible for our results, which result from contact inhibition of chemotaxis, are both generic developmental mechanisms and interesting examples of unusual patterning instabilities., Comment: Thoroughly revised version, now in press in PLoS Computational Biology. 53 pages, 13 figures, 2 supporting figures, 56 supporting movies, source code and parameters files for computer simulations provided. Supporting information: http://www.psb.ugent.be/~romer/ploscompbiol/ Source code: http://sourceforge.net/projects/tst/

The Cellular Potts Model (CPM) is a robust, cell-level methodology for simulation of biological tissues and morphogenesis. Both tissue physiology and morphogenesis depend on diffusion of chemical morphogens in the extra-cellular fluid or matrix (ECM). Standard diffusion solvers applied to the cellular potts model use finite difference methods on the underlying CPM lattice. However, these methods produce a diffusing field tied to the underlying lattice, which is inaccurate in many biological situations in which cell or ECM movement causes advection rapid compared to diffusion. Finite difference schemes suffer numerical instabilities solving the resulting advection-diffusion equations. To circumvent these problems we simulate advection-diffusion within the framework of the CPM using off-lattice finite-difference methods. We define a set of generalized fluid particles which detach advection and diffusion from the lattice. Diffusion occurs between neighboring fluid particles by local averaging rules which approximate the Laplacian. Directed spin flips in the CPM handle the advective movement of the fluid particles. A constraint on relative velocities in the fluid explicitly accounts for fluid viscosity. We use the CPM to solve various diffusion examples including multiple instantaneous sources, continuous sources, moving sources and different boundary geometries and conditions to validate our approximation against analytical and established numerical solutions. We also verify the CPM results for Poiseuille flow and Taylor-Aris dispersion.

We discuss a method to reconstruct an approximate two-dimensional foam structure from an incomplete image using the extended Potts mode with a pinned lattice we introduced in a previous paper. The initial information consists of the positions of the vertices only. We locate the centers of the bubbles using the Euclidean distance-map construction and assign at each vertex position a continuous pinning field with a potential falling off as $1/r$. We nucleate a bubble at each center using the extended Potts model and let the structure evolve under the constraint of scaled target areas until the bubbles contact each other. The target area constraint and pinning centers prevent further coarsening. We then turn the area constraint off and let the edges relax to a minimum energy configuration. The result is a reconstructed structure very close to the simulation. We repeated this procedure for various stages of the coarsening of the same simulated foam and investigated the simulation and reconstruction dynamics, topology and area distribution, finding that they agree to good accuracy., Comment: 31 pages, 20 Postscript figures Accepted in the Journal of Computational Physics

In a continuum description of materials, the stress tensor field $\bar{% \bar{\sigma}}$ quantifies the internal forces the neighbouring regions exert on a region of the material. The classical theory of elastic solids assumes that $\bar{\bar{\sigma}}$ determines the strain, while hydrodynamics assumes that $\bar{\bar{\sigma}}$ determines the strain rate. To extend both successful theories to more general materials, which display both elastic and fluid properties, we recently introduced a descriptor generalizing the classical strain to include plastic deformations: the ``statistical strain'', based on averages on microscopic details (``A texture tensor to quantify deformations'' M.Au., Y.J., J.A.G., F.G, companion paper, {\em Granular Matter}, same issue). Here, we apply such a statistical analysis to a two-dimensional foam steadily flowing through a constriction, a problem beyond reach of both theories, and prove that the foam has the elastic properties of a (linear and isotropic) continuous medium., Comment: 4 pages, 5 figures. To appear in Granular Matter

Under mechanical deformation, most materials exhibit both elastic and fluid (or plastic) responses. No existing formalism derived from microscopic principles encompasses both their fluid-like and solid-like aspects. We define the {\it statistical texture tensor} to quantify the intuitive notion of stored deformation. This tensor links microscopic and macroscopic descriptions of the material, and extends the definition of elastic strain., Comment: 4 pages, 2 figures. to appear in Granular Matter

In this paper we argue that the pattern of cell movements in the morphogenetic process known as convergent extension can be understood as a energy minimization process, provided the cell-cell adhesive energy has a certain kind of anisotropy. This single simple property is sufficient cause for the type of cell elongation, alignment, and intercalation of a cellular array that is the characteristic of convergent extension. We describe the type of anisotropy required. We show that the final aspect ratio of the array of cells is independent of the initial configuration and of the degree of cell elongation. We find how it depends on the anisotropy., Comment: Eight pages, four EPS figures, LaTeX: artilce class, graphics package, Addition: correction of three small typos, NEW: additional author; fewer, smaller figures; minor typos