Coarse-grained Continuum Model Research Articles

Reliable and efficient parameter estimation using approximate continuum limit descriptions of stochastic models

Stochastic individual-based mathematical models are attractive for modelling biological phenomena because they naturally capture the stochasticity and variability that is often evident in biological data. Such models also allow us to track the motion of individuals within the population of interest. Unfortunately, capturing this microscopic detail means that simulation and parameter inference can become computationally expensive. One approach for overcoming this computational limitation is to coarse-grain the stochastic model to provide an approximate continuum model that can be solved using far less computational effort. However, coarse-grained continuum models can be biased or inaccurate, particularly for certain parameter regimes. In this work, we combine stochastic and continuum mathematical models in the context of lattice-based models of two-dimensional cell biology experiments by demonstrating how to simulate two commonly used experiments: cell proliferation assays and barrier assays. Our approach involves building a simple statistical model of the discrepancy between the expensive stochastic model and the associated computationally inexpensive coarse-grained continuum model. We form this statistical model based on a limited number of expensive stochastic model evaluations at design points sampled from a user-chosen distribution, corresponding to a computer experiment design problem. With straightforward design point selection schemes, we show that using the statistical model of the discrepancy in tandem with the computationally inexpensive continuum model allows us to carry out prediction and inference while correcting for biases and inaccuracies due to the continuum approximation. We demonstrate this approach by simulating a proliferation assay, where the continuum limit model is the well-known logistic ordinary differential equation, as well as a barrier assay where the continuum limit model is closely related to the well-known Fisher-KPP partial differential equation. We construct an approximate likelihood function for parameter inference, both with and without discrepancy correction terms. Using maximum likelihood estimation, we provide point estimates of the unknown parameters, and use the profile likelihood to characterise the uncertainty in these estimates and form approximate confidence intervals. For the range of inference problems considered, working with the continuum limit model alone leads to biased parameter estimation and confidence intervals with poor coverage. In contrast, incorporating correction terms arising from the statistical model of the model discrepancy allows us to recover the parameters accurately with minimal computational overhead. The main tradeoff is that the associated confidence intervals are typically broader, reflecting the additional uncertainty introduced by the approximation process. All algorithms required to replicate the results in this work are written in the open source Julia language and are available at GitHub.

A fast Chebyshev method for the Bingham closure with application to active nematic suspensions

Continuum kinetic theories provide an important tool for the analysis and simulation of particle suspensions. When those particles are anisotropic, the addition of a particle orientation vector to the kinetic description yields a 2d−1 dimensional theory which becomes intractable to simulate, especially in three dimensions or near states where the particles are highly aligned. Coarse-grained theories that track only moments of the particle distribution functions provide a more efficient simulation framework, but require closure assumptions. For the particular case where the particles are apolar, the Bingham closure has been found to agree well with the underlying kinetic theory; yet the closure is non-trivial to compute, requiring the solution of an often nearly-singular nonlinear equation at every spatial discretization point at every timestep. In this paper, we present a robust, accurate, and efficient numerical scheme for evaluating the Bingham closure, with a controllable error/efficiency tradeoff. To demonstrate the utility of the method, we carry out high-resolution simulations of a coarse-grained continuum model for a suspension of active particles in parameter regimes inaccessible to kinetic theories. Analysis of these simulations reveals that inaccurately computing the closure can act to effectively limit spatial resolution in the coarse-grained fields. Pushing these simulations to the high spatial resolutions enabled by our method reveals a coupling between vorticity and topological defects in the suspension director field, as well as signatures of energy transfer between scales in this active fluid model.

Combined Force-Torque Spectroscopy of Proteins by Means of Multiscale Molecular Simulation

Assessing the structural properties of large proteins is important to gain an understanding of their function in, e.g., biological systems or biomedical applications. We propose a method to examine the mechanical properties of proteins subject to applied forces by means of multiscale simulation. Both stretching and torsional forces are considered, and these may be applied independently of each other. As a proof of principle, we apply torsional forces to a coarse-grained continuum model of the antibody protein immunoglobulin G using fluctuating finite element analysis and use it to identify the area of strongest deformation. This region is essential to the torsional properties of the molecule as a whole because it represents the softest, most deformable domain. Zooming in, this part of the molecule is subjected to torques and stretching forces using molecular dynamics simulations on an atomistically resolved level to investigate its torsional properties. We calculate the torsional resistance as a function of the rotation of the domain while subjecting it to various stretching forces. From this, we assess how the measured twist-torque profiles develop with increasing stretching force and show that they exhibit torsion stiffening, in qualitative agreement with experimental findings. We argue that combining the twist-torque profiles for various stretching forces effectively results in a combined force-torque spectroscopy analysis, which may serve as a mechanical signature for a biological macromolecule.

Multi-scale coarse-grained continuum models for bifurcation and transient analysis of coupled homogeneous-catalytic reactions in monoliths

We use the Lyapunov-Schmidt (L-S) reduction technique of bifurcation theory to develop reduced order models for monolith reactors with catalytic reactions in the washcoat layer and homogeneous reactions in the flow channel, and with thermal effects. The L-S technique eliminates the two transverse coordinates and converts the detailed models expressed in terms of three spatial coordinates and time to multi-mode models in axial coordinate and time. The reduced order models are expressed in terms of phase averaged multiple concentration (temperature) modes and interfacial fluxes which are related through various transfer coefficients that are local scale, flow and property dependent. Accurate expressions are provided for estimating the local transfer coefficient matrices for multi-component systems. Various limiting cases of the multi-mode reduced order models are discussed and compared to other models in the literature. The usefulness of the reduced order models is illustrated by two examples: the first is the oxidative coupling of methane in a laboratory scale monolith reactor, the second is the determination of the ignition-extinction behavior of catalytically assisted combustion of propane. Other possible applications and extensions are also discussed.

Mean field approach of dynamical pattern formation in underdamped active matter with short-ranged alignment and distant anti-alignment interactions

Many active matter systems, especially on the microscopic scale, are well approximated as overdamped, meaning that any inertial momentum is immediately dissipated by the environment. On the other hand, especially for macroscopic active systems but also for many mesoscopic systems the time scale of translational inertial motion can become large enough to be relevant for the dynamics. This raises the question how collective dynamics and the resulting states in active matter are influenced by inertia. Therefore, we propose a coarse-grained continuum model for underdamped active matter based on a mean field description for passive systems. Furthermore, as an example, we apply the model to a system with interactions that support an alignment on short distances and an anti-alignment on longer length scales as known in the context of pattern formation due to orientational interactions. Our numerical calculations of the under- and overdamped dynamics both predict a structured laning state. However, activity induced convective flows that are only present in the underdamped model destabilize this state when the anti-alignment is weakened, leading to a collective motion state which does not occur in the overdamped limit. A turbulent transition regime between the two states can be characterized by strong density fluctuations and the absence of global ordering.

Sliding Filament and Fixed Filament Mechanisms Contribute to Tension of the Fission Yeast Cytokinetic Ring

A fundamental challenge in cell biology is understanding how cells generate contractile force from actomyosin machineries. The contractile ring is a vital actomyosin machine that generates tension and constricts and divides the cell during cytokinesis. Following discovery of the contractile ring decades ago, a sliding filament mechanism for ring tension was proposed similar to that in muscle (Schroeder, 1972). However, contractile rings lack the highly ordered periodic sarcomeric architecture of muscle, and it remains unknown how tension emerges from the disordered ring organization. Here, we developed a highly coarse-grained continuum model of the fission yeast cytokinetic ring, severely constrained by experimentally measured component amounts. The model incorporates the two myosin-II isoforms, Myo2 in membrane-anchored protein complexes called nodes, and Myp2 clusters assumed unanchored. We found that the model does indeed predict a sliding filament mechanism, but somewhat different to that envisaged by the pioneering studies of Schroeder and others. The sliding filament mechanism is a spatially and temporally homogeneous version of that in muscle: pairs of oppositely oriented actin filaments are made tense as they are pulled toward one another by myosin-II clusters through which they slide. The mechanism relies on anchoring of actin filament barbed ends to the plasma membrane, which enables filaments to become tense when pulled. A second fixed filament component is generated by chains of like-oriented actin filaments encircling the ring, independent of lateral anchoring. Cellular actomyosin contractile machines are intrinsically vulnerable to contractile instabilities. We also studied how such instabilities are controlled. We find that tension loss due to runaway catastrophic aggregation of the contractile components Myo2 and Myp2 is prevented by fast turnover and sufficiently low myosin cluster mobilities, enforced by lateral anchoring and crowding.

Coarse graining the dynamics of immersed and driven fiber assemblies

An important class of fluid-structure problems involve the dynamics of ordered arrays of immersed, flexible fibers. While specialized numerical methods have been developed to study fluid-fiber systems, they become infeasible when there are many, rather than a few, fibers present, nor do these methods lend themselves to analytical calculation. Here, we introduce a coarse-grained continuum model, based on local-slender body theory, for elastic fibers immersed in a viscous Newtonian fluid. It takes the form of an anisotropic Brinkman equation whose skeletal drag is coupled to elastic forces. This model has two significant benefits: (1) the density effects of the fibers in a suspension become analytically manifest, and (2) it allows for the rapid simulation of dense suspensions of fibers in regimes inaccessible to standard methods. As a first validation, without fitting parameters, we achieve very reasonable agreement with 3D Immersed Boundary simulations of a bed of anchored fibers bent by a shear flow. Secondly, we characterize the effect of density on the relaxation time of fiber beds under oscillatory shear, and find close agreement to results from full numerical simulations. We then study buckling instabilities in beds of fibers, using our model both numerically and analytically to understand the role of fiber density and the structure of buckling transitions. We next apply our model to study the flow-induced bending of inclined fibers in a channel, as has been recently studied as a flow rectifier, examining the nature of the internal flows within the bed, and the emergence of inhomogeneous permeability. Finally, we extend the method to study a simple model of metachronal waves on beds of actuated fibers, as a model for ciliary beds. Our simulations reproduce qualitatively the pumping action of coordinated waves of compression through the bed.

Multiscale modelling of charged impurities in two-dimensional materials

Charged impurities influence functional properties of two-dimensional materials and a detailed theoretical understanding of charged defects is required to enable a rational design of defect-engineered nanomaterials for applications in ultrathin devices. To achieve this goal, we have developed multiscale approaches that combine atomistic first-principles theories, such as density-functional theory, with coarse-grained continuum models, such as effective mass models. This allows us to model large supercells which are required to accurately describe the slow decay of the screened defect potential and the defect-induced changes in the electronic properties of the two-dimensional host material. I will describe the results of our multiscale calculations for charged defects in doped graphene and in transition-metal dichalcogenide monolayers which have revealed novel mechanisms for controlling and tuning the electronic structure of two-dimensional materials.

Computer Simulation of Protein Materials at Multiple Length Scales: From Single Proteins to Protein Assemblies

Computer simulation of protein materials and their dynamic or mechanical behavior is of high significance, as proteins perform their functions through their structural changes in response to a force (or stimulus). The computer simulation enables the detailed insight into the structure and behavior of proteins at atomistic resolution, which is inaccessible with experimental toolkits such as single-molecule experiments. With the advancement of computing resources, the computer simulation has recently played a vital role as a virtual microscopy in understanding and characterizing the structure and behaviors of protein materials at atomistic resolution. For examples, computer simulations allow for gaining insight into how some protein domains can exhibit the remarkable mechanical properties and functions. In this article, we would like to address the current state-of-arts of computer simulations that have been employed for studying the structure and properties (or behaviors) of proteins at multiple length scales ranging from single proteins to protein assemblies. Specifically, we summarize various computational modeling techniques, ranging from atomistic models to coarse-grained models and continuum models, applicable for modeling protein structures at multiple length scales from single proteins to protein assemblies. This paper discusses how such various computational modeling/simulation techniques can be employed for studying the structure and properties of protein materials at multiple length scales. This paper sheds light on computer simulations, which are able to unveil the hidden, complex mechanisms related to the structure and properties of protein materials at multiple length scales.

Computer Simulation Studies of Polymer Adsorption and Aggregation – From Flexible to Stiff

This mini-review gives a brief overview of our recent Monte Carlo computer simulations of polymer adsorption and aggregation based on coarse-grained continuum models of theta polymers focusing on excluded volume, short-range attraction and bending stiffness. The simulations are performed in generalized ensembles (multicanonical and parallel tempering) and analyzed with canonical and microcanonical methods. For flexible polymers, besides plane adsorbing surfaces also the curved inner wall of a spherical cage has been considered. For semiflexible polymers, the dependence on bending stiffness is quite intricate and was recently found to lead to phases with stable knots of unique characteristics. For the special case of P3HT chains interacting with a gold surface, we have employed a chemically realistic coarse-grained model and compared the simulation results with scanning tunneling microscopy (STM) experiments under ultra-high vacuum conditions. The paper concludes with recent results on the aggregation process of semiflexible polymers. Our data shows that the bending stiffness plays a key role as to whether the polymer system forms amorphous aggregates or twisted bundle-like motifs.

On the acceleration of spatially distributed agent-based computations: A patch dynamics scheme

In recent years, individual-based/agent-based modeling has been applied to study a wide range of applications, ranging from engineering problems to phenomena in sociology, economics and biology. Simulating such agent-based models over extended spatiotemporal domains can be prohibitively expensive due to stochasticity and the presence of multiple scales. Nevertheless, many agent-based problems exhibit smooth behavior in space and time on a macroscopic scale, suggesting that a useful coarse-grained continuum model could be obtained. For such problems, the equation-free framework [16–18] can significantly reduce the computational cost. Patch dynamics is an essential component of this framework. This scheme is designed to perform numerical simulations of an unavailable macroscopic equation on macroscopic time and length scales; it uses appropriately initialized simulations of the fine-scale agent-based model in a number of small “patches”, which cover only a fraction of the spatiotemporal domain. In this work, we construct a finite-volume-inspired conservative patch dynamics scheme and apply it to a financial market agent-based model based on the work of Omurtag and Sirovich [22]. We first apply our patch dynamics scheme to a continuum approximation of the agent-based model, to study its performance and analyze its accuracy. We then apply the scheme to the agent-based model itself. Our computational experiments indicate that here, typically, the patch dynamics-based simulation needs to be performed in only 20% of the full agent simulation space, and in only 10% of the temporal domain.

Adsorption of self-avoiding tethered membranes: A Monte Carlo simulation study

Adsorption of self-avoiding tethered membranes of hexagonal orientation on a solid impenetrable plane is studied by means of Monte Carlo computer simulations of a coarse-grained continuum model, varying the membrane linear size L and the strength epsilon of the short-range attractive wall potential. A second-order adsorption transition is found to take place at a critical strength of the adsorption potential epsilon(c), as predicted earlier for binding manifolds in the so-called strong fluctuating regime. By means of finite-size scaling analysis for membranes of size 5<or=L<or=30, containing 61<or=N<or=2611 monomers, we find that the critical crossover exponent for adsorption phi approximately 0.60+/-0.01. Thus the fraction m of adsorbed segments at epsilon(c) is found to scale as m proportional to N(phi-1)=N(-0.4). The membrane thickness lambda(min) decreases with growing strength epsilon of the adhesive potential as lambda(min) proportional to mid R:epsilon/epsilon(c)-1mid R:(-psi), where psi approximately 0.58+/-0.02. The monomer density profiles of adsorbed membranes decay exponentially with the distance z from the substrate rho(z) proportional to exp[-(z/xi)], where the correlation length xi proportional to mid R:epsilon/epsilon(c)-1mid R:(-psi) with psi approximately 0.70+/-0.01. We also investigate the kinetics of adsorption of a polymerized membrane in the regime of strong adsorption and find that the order parameter variation with elapsed time during the adsorption process is given by a power law m(t) proportional to t(omega), where omega approximately 1.0, regardless of the strength of the adsorbing potential epsilon>>epsilon(c). The characteristic time for complete adsorption in this regime scales as tau(ads) proportional to L(2). Regarding the strength epsilon of the substrate potential, tau(ads) is found to diminish linearly as the respective equilibrium value of the order parameter increases. A simple analytic model yields also tau(ads) proportional to L(2) while suggesting that the process of adsorption cannot be identified by simple "unrolling."

Dislocation density-based constitutive model for the mechanical behaviour of irradiated Cu

Performance degradation of structural steels in nuclear environments results from the formation of a high number density of nanometre-scale defects. The defects observed in copper-based alloys are composed of vacancy clusters in the form of stacking fault tetrahedra and/or prismatic dislocation loops that impede the motion of dislocations. The mechanical behaviour of irradiated copper alloys exhibits increased yield strength, decreased total strain to failure and decreased work hardening as compared to their unirradiated behaviour. Above certain critical defect concentrations (neutron doses), the mechanical behaviour exhibits distinct upper yield points. In this paper, we describe the formulation of an internal state variable model for the mechanical behaviour of such materials subject to these (irradiation) environments. This model has been developed within a multiscale materials-modelling framework, in which molecular dynamics simulations of dislocation–radiation defect interactions inform the final coarse-grained continuum model. The plasticity model includes mechanisms for dislocation density growth and multiplication and for irradiation defect density evolution with dislocation interaction. The general behaviour of the constitutive (homogeneous material point) model shows that as the defect density increases, the initial yield point increases and the initial strain hardening decreases. The final coarse-grained model is implemented into a finite element framework and used to simulate the behaviour of tensile specimens with varying levels of irradiation-induced material damage. The simulation results compare favourably with the experimentally observed mechanical behaviour of irradiated materials.

Coarse-grained continuum model for molecular diffusion in a lipid monolayer

Abstract Recent experimental studies of lateral diffusion of lipid molecules in bilayers using single-molecule detection techniques (Schutz et al., 1997. Biophysical J. 73:1073–1080), gave rise to speculation about confinement effects due to anomalous diffusion. Here we report a computer simulation study of the lateral diffusion of conformationally disordered lipid monolecules in a monolayer structure. No evidence for anomalous diffusion effects are found in the homogeneous monolayer studied. The classical Cohen-Turnbull theory is found to provide a good description of the simulated lateral diffusion coefficients at high densities, above a certain threshold c ∗ . At densities c ∗ , where the diffusion of the flexible lipids is governed by their “soft core” repulsion, we found an entropically activated diffusion according to D ∼ exp (− γ √ c ) over a large concentration regime 0.03 c

Theoretical modeling of Langmuir monolayers

We study coarse-grained continuum models for Langmuir monolayers by self-consistent field theory and by Monte Carlo simulations. Amphiphilic molecules are represented by stiff chains of monomers with one end grafted to a planar surface. In particular, we discuss the origin of successive fluid–fluid transitions, the possible origin of tilt order and the factors which determine the direction of tilt.