Asymmetric Periodic Boundary Conditions for All-Atom Molecular Dynamics and Coarse-Grained Simulations of Nucleic Acids.

Periodic boundary conditions are commonly applied in molecular dynamics simulations in the microcanonical (NVE), canonical (NVT), and isothermal-isobaric (NpT) ensembles. In their simplest application, a biological system of interest is placed in the middle of a solvation box, which is chosen 'sufficiently large' to minimize any numerical artifacts associated with the periodic boundary conditions. This practical approach brings limitations to the size of biological systems that can be simulated. Here, we study simulations of effectively infinitely long nucleic acids, which are solvated in the directions perpendicular to the polymer chain, while periodic boundary conditions are also applied along the polymer chain. We study the effects of these asymmetric periodic boundary conditions (APBC) on the simulated results, including the mechanical properties of biopolymers and the properties of the surrounding solvent. To get some further insights into the advantages of using the APBC, a coarse-grained worm-like chain model is first studied, illustrating how the persistence length can be extracted from the local properties of the polymer chain, which are less affected by the APBC than some global averages. This is followed by all-atom molecular dynamics simulations of DNA in ionic solutions, where we use the APBC to investigate sequence-dependent properties of DNA molecules and properties of the surrounding solvent.

Symmetries of many-body systems imply distance-dependent potentials.

Considering an interatomic potential U(q), where q=[q_{1},q_{2},⋯,q_{N}]∈R^{3N} is a vector describing positions q_{i}∈R^{3}, it is shown that U can be defined as a function of the interatomic distance variables r_{ij}=|q_{i}-q_{j}| provided the potential U satisfies some symmetry assumptions. Moreover, the potential U can be defined as a function of a proper subset of the distance variables r_{ij}, provided N>5, with the number of distance variables used scaling linearly with the number of atoms N.

Chemical Systems with Limit Cycles.

The dynamics of a chemical reaction network (CRN) is often modeled under the assumption of mass action kinetics by a system of ordinary differential equations (ODEs) with polynomial right-hand sides that describe the time evolution of concentrations of chemical species involved. Given an arbitrarily large integer [Formula: see text], we show that there exists a CRN such that its ODE model has at least K stable limit cycles. Such a CRN can be constructed with reactions of at most second-order provided that the number of chemical species grows linearly with K. Bounds on the minimal number of chemical species and the minimal number of chemical reactions are presented for CRNs with K stable limit cycles and at most second order or seventh-order kinetics. We also show that CRNs with only two chemical species can have K stable limit cycles, when the order of chemical reactions grows linearly with K.

On Stretching, Bending, Shearing, and Twisting of Actin Filaments I: Variational Models.

Mechanochemical simulations of actomyosin networks are traditionally based on one-dimensional models of actin filaments having zero width. Here, and in the follow up paper (arXiv, DOI 10.48550/arXiv.2203.01284), approaches are presented for more efficient modeling that incorporates stretching, shearing, and twisting of actin filaments. Our modeling of a semiflexible filament with a small but finite width is based on the Cosserat theory of elastic rods, which allows for six degrees of freedom at every point on the filament's backbone. In the variational models presented in this paper, a small and discrete set of parameters is used to describe a smooth filament shape having all degrees of freedom allowed in the Cosserat theory. Two main approaches are introduced: one where polynomial spline functions describe the filament's configuration, and one in which geodesic curves in the space of the configurational degrees of freedom are used. We find that in the latter representation the strain energy function can be calculated without resorting to a small-angle expansion, so it can describe arbitrarily large filament deformations without systematic error. These approaches are validated by a dynamical model of a Cosserat filament, which can be further extended by using multiresolution methods to allow more detailed monomer-based resolution in certain parts of the actin filament, as introduced in the follow up paper. The presented framework is illustrated by showing how torsional compliance in a finite-width filament can induce broken chiral symmetry in the structure of a cross-linked bundle.

On standardised moments of force distribution in simple liquids.

The force distribution of a tagged atom in a Lennard-Jones fluid in the canonical ensemble is studied with a focus on its dependence on inherent physical parameters: number density (n) and temperature (T). Utilising structural information from molecular dynamics simulations of the Lennard-Jones fluid, explicit analytical expressions for the dependence of standardised force moments on n and T are derived. Leading order behaviour of standardised moments of the force distribution are obtained in the limiting cases of small density (n → 0) and low temperature (T → 0), while the variations in the standardised moments are probed for general n and T using molecular dynamics simulations. Clustering effects are seen in molecular dynamics simulations and their effect on these standardised moments is discussed.

Coarse-graining molecular dynamics: stochastic models with non-Gaussian force distributions.

Incorporating atomistic and molecular information into models of cellular behaviour is challenging because of a vast separation of spatial and temporal scales between processes happening at the atomic and cellular levels. Multiscale or multi-resolution methodologies address this difficulty by using molecular dynamics (MD) and coarse-grained models in different parts of the cell. Their applicability depends on the accuracy and properties of the coarse-grained model which approximates the detailed MD description. A family of stochastic coarse-grained (SCG) models, written as relatively low-dimensional systems of nonlinear stochastic differential equations, is presented. The nonlinear SCG model incorporates the non-Gaussian force distribution which is observed in MD simulations and which cannot be described by linear models. It is shown that the nonlinearities can be chosen in such a way that they do not complicate parametrization of the SCG description by detailed MD simulations. The solution of the SCG model is found in terms of gamma functions.

Limiting stochastic processes of shift-periodic dynamical systems.

A shift-periodic map is a one-dimensional map from the real line to itself which is periodic up to a linear translation and allowed to have singularities. It is shown that iterative sequences x n+1 = F(x n ) generated by such maps display rich dynamical behaviour. The integer parts ⌊ x n ⌋ give a discrete-time random walk for a suitable initial distribution of x 0 and converge in certain limits to Brownian motion or more general Lévy processes. Furthermore, for certain shift-periodic maps with small holes on [0,1], convergence of trajectories to a continuous-time random walk is shown in a limit.

Multiscale Stochastic Reaction-Diffusion Algorithms Combining Markov Chain Models with Stochastic Partial Differential Equations.

Two multiscale algorithms for stochastic simulations of reaction-diffusion processes are analysed. They are applicable to systems which include regions with significantly different concentrations of molecules. In both methods, a domain of interest is divided into two subsets where continuous-time Markov chain models and stochastic partial differential equations (SPDEs) are used, respectively. In the first algorithm, Markov chain (compartment-based) models are coupled with reaction-diffusion SPDEs by considering a pseudo-compartment (also called an overlap or handshaking region) in the SPDE part of the computational domain right next to the interface. In the second algorithm, no overlap region is used. Further extensions of both schemes are presented, including the case of an adaptively chosen boundary between different modelling approaches.

Multi-resolution dimer models in heat baths with short-range and long-range interactions.

This work investigates multi-resolution methodologies for simulating dimer models. The solvent particles which make up the heat bath interact with the monomers of the dimer either through direct collisions (short-range) or through harmonic springs (long-range). Two types of multi-resolution methodologies are considered in detail: (a) describing parts of the solvent far away from the dimer by a coarser approach; (b) describing each monomer of the dimer by using a model with different level of resolution. These methodologies are then used to investigate the effect of a shared heat bath versus two uncoupled heat baths, one for each monomer. Furthermore, the validity of the multi-resolution methods is discussed by comparison to dynamics of macroscopic Langevin equations.

Multi-resolution polymer Brownian dynamics with hydrodynamic interactions.

A polymer model given in terms of beads, interacting through Hookean springs and hydrodynamic forces, is studied. A Brownian dynamics description of this bead-spring polymer model is extended to multiple resolutions. Using this multiscale approach, a modeller can efficiently look at different regions of the polymer in different spatial and temporal resolutions with scalings given for the number of beads, statistical segment length, and bead radius in order to maintain macro-scale properties of the polymer filament. The Boltzmann distribution of a Gaussian chain for differing statistical segment lengths gives a diffusive displacement equation for the multi-resolution model with a mobility tensor for different bead sizes. Using the pre-averaging approximation, the translational diffusion coefficient is obtained as a function of the inverse of a matrix and then in closed form in the long-chain limit. This is then confirmed with numerical experiments.

Noise control for molecular computing.

Synthetic biology is a growing interdisciplinary field, with far-reaching applications, which aims to design biochemical systems that behave in a desired manner. With the advancement in nucleic-acid-based technology in general, and strand-displacement DNA computing in particular, a large class of abstract biochemical networks may be physically realized using nucleic acids. Methods for systematic design of the abstract systems with prescribed behaviours have been predominantly developed at the (less-detailed) deterministic level. However, stochastic effects, neglected at the deterministic level, are increasingly found to play an important role in biochemistry. In such circumstances, methods for controlling the intrinsic noise in the system are necessary for a successful network design at the (more-detailed) stochastic level. To bridge the gap, the noise-control algorithm for designing biochemical networks is developed in this paper. The algorithm structurally modifies any given reaction network under mass-action kinetics, in such a way that (i) controllable state-dependent noise is introduced into the stochastic dynamics, while (ii) the deterministic dynamics are preserved. The capabilities of the algorithm are demonstrated on a production-decay reaction system, and on an exotic system displaying bistability. For the production-decay system, it is shown that the algorithm may be used to redesign the network to achieve noise-induced multistability. For the exotic system, the algorithm is used to redesign the network to control the stochastic switching, and achieve noise-induced oscillations.

Reaction time for trimolecular reactions in compartment-based reaction-diffusion models.

Trimolecular reaction models are investigated in the compartment-based (lattice-based) framework for stochastic reaction-diffusion modeling. The formulae for the first collision time and the mean reaction time are derived for the case where three molecules are present in the solution under periodic boundary conditions. For the case of reflecting boundary conditions, similar formulae are obtained using a computer-assisted approach. The accuracy of these formulae is further verified through comparison with numerical results. The presented derivation is based on the first passage time analysis of Montroll [J. Math. Phys. 10, 753 (1969)]. Montroll's results for two-dimensional lattice-based random walks are adapted and applied to compartment-based models of trimolecular reactions, which are studied in one-dimensional or pseudo one-dimensional domains.

Hard-sphere interactions in velocity-jump models.

Group-level behavior of particles undergoing a velocity-jump process with hard-sphere interactions is investigated. We derive N-particle transport equations that include the possibility of collisions between particles and apply different approximation techniques to get expressions for the dependence of the collective diffusion coefficient on the number of particles and their diameter. The derived approximations are compared with numerical results obtained from individual-based simulations. The theoretical results compare well with Monte Carlo simulations providing the excluded-volume fraction is small.

Steric Effects Induce Geometric Remodeling of Actin Bundles in Filopodia.

Filopodia are ubiquitous fingerlike protrusions, spawned by many eukaryotic cells, to probe and interact with their environments. Polymerization dynamics of actin filaments, comprising the structural core of filopodia, largely determine their instantaneous lengths and overall lifetimes. The polymerization reactions at the filopodial tip require transport of G-actin, which enter the filopodial tube from the filopodial base and diffuse toward the filament barbed ends near the tip. Actin filaments are mechanically coupled into a tight bundle by cross-linker proteins. Interestingly, many of these proteins are relatively short, restricting the free diffusion of cytosolic G-actin throughout the bundle and, in particular, its penetration into the bundle core. To investigate the effect of steric restrictions on G-actin diffusion by the porous structure of filopodial actin filament bundle, we used a particle-based stochastic simulation approach. We discovered that excluded volume interactions result in partial and then full collapse of central filaments in the bundle, leading to a hollowed-out structure. The latter may further collapse radially due to the activity of cross-linking proteins, hence producing conical-shaped filament bundles. Interestingly, electron microscopy experiments on mature filopodia indeed frequently reveal actin bundles that are narrow at the tip and wider at the base. Overall, our work demonstrates that excluded volume effects in the context of reaction-diffusion processes in porous networks may lead to unexpected geometric growth patterns and complicated, history-dependent dynamics of intermediate metastable configurations.

Coupling all-atom molecular dynamics simulations of ions in water with Brownian dynamics.

Molecular dynamics (MD) simulations of ions (K+, Na+, Ca2+ and Cl-) in aqueous solutions are investigated. Water is described using the SPC/E model. A stochastic coarse-grained description for ion behaviour is presented and parametrized using MD simulations. It is given as a system of coupled stochastic and ordinary differential equations, describing the ion position, velocity and acceleration. The stochastic coarse-grained model provides an intermediate description between all-atom MD simulations and Brownian dynamics (BD) models. It is used to develop a multiscale method which uses all-atom MD simulations in parts of the computational domain and (less detailed) BD simulations in the remainder of the domain.

Convergence of methods for coupling of microscopic and mesoscopic reaction-diffusion simulations.

In this paper, three multiscale methods for coupling of mesoscopic (compartment-based) and microscopic (molecular-based) stochastic reaction-diffusion simulations are investigated. Two of the three methods that will be discussed in detail have been previously reported in the literature; the two-regime method (TRM) and the compartment-placement method (CPM). The third method that is introduced and analysed in this paper is called the ghost cell method (GCM), since it works by constructing a "ghost cell" in which molecules can disappear and jump into the compartment-based simulation. Presented is a comparison of sources of error. The convergent properties of this error are studied as the time step Δt (for updating the molecular-based part of the model) approaches zero. It is found that the error behaviour depends on another fundamental computational parameter h, the compartment size in the mesoscopic part of the model. Two important limiting cases, which appear in applications, are considered: (i) Δt → 0 and h is fixed; (ii) Δt → 0 and h → 0 such that √Δt/h is fixed. The error for previously developed approaches (the TRM and CPM) converges to zero only in the limiting case (ii), but not in case (i). It is shown that the error of the GCM converges in the limiting case (i). Thus the GCM is superior to previous coupling techniques if the mesoscopic description is much coarser than the microscopic part of the model.

Tensor methods for parameter estimation and bifurcation analysis of stochastic reaction networks.

Stochastic modelling of gene regulatory networks provides an indispensable tool for understanding how random events at the molecular level influence cellular functions. A common challenge of stochastic models is to calibrate a large number of model parameters against the experimental data. Another difficulty is to study how the behaviour of a stochastic model depends on its parameters, i.e. whether a change in model parameters can lead to a significant qualitative change in model behaviour (bifurcation). In this paper, tensor-structured parametric analysis (TPA) is developed to address these computational challenges. It is based on recently proposed low-parametric tensor-structured representations of classical matrices and vectors. This approach enables simultaneous computation of the model properties for all parameter values within a parameter space. The TPA is illustrated by studying the parameter estimation, robustness, sensitivity and bifurcation structure in stochastic models of biochemical networks. A Matlab implementation of the TPA is available at http://www.stobifan.org.

Noise-induced multistability in chemical systems: Discrete versus continuum modeling.

The noisy dynamics of chemical systems is commonly studied using either the chemical master equation (CME) or the chemical Fokker-Planck equation (CFPE). The latter is a continuum approximation of the discrete CME approach. It has recently been shown that for a particular system, the CFPE captures noise-induced multistability predicted by the CME. This phenomenon involves the CME's marginal probability distribution changing from unimodal to multimodal as the system size decreases below a critical value. We here show that the CFPE does not always capture noise-induced multistability. In particular we find simple chemical systems for which the CME predicts noise-induced multistability, whereas the CFPE predicts monostability for all system sizes.

Multiscale reaction-diffusion simulations with Smoldyn.

UNLABELLED: Smoldyn is a software package for stochastic modelling of spatial biochemical networks and intracellular systems. It was originally developed with an accurate off-lattice particle-based model at its core. This has recently been enhanced with the addition of a computationally efficient on-lattice model, which can be run stand-alone or coupled together for multiscale simulations using both models in regions where they are most required, increasing the applicability of Smoldyn to larger molecule numbers and spatial domains. Simulations can switch between models with only small additions to their configuration file, enabling users with existing Smoldyn configuration files to run the new on-lattice model with any reaction, species or surface descriptions they might already have. AVAILABILITY AND IMPLEMENTATION: Source code and binaries freely available for download at www.smoldyn.org, implemented in C/C++ and supported on Linux, Mac OSX and MS Windows.

Stochastic Turing patterns: analysis of compartment-based approaches.

Turing patterns can be observed in reaction-diffusion systems where chemical species have different diffusion constants. In recent years, several studies investigated the effects of noise on Turing patterns and showed that the parameter regimes, for which stochastic Turing patterns are observed, can be larger than the parameter regimes predicted by deterministic models, which are written in terms of partial differential equations (PDEs) for species concentrations. A common stochastic reaction-diffusion approach is written in terms of compartment-based (lattice-based) models, where the domain of interest is divided into artificial compartments and the number of molecules in each compartment is simulated. In this paper, the dependence of stochastic Turing patterns on the compartment size is investigated. It has previously been shown (for relatively simpler systems) that a modeler should not choose compartment sizes which are too small or too large, and that the optimal compartment size depends on the diffusion constant. Taking these results into account, we propose and study a compartment-based model of Turing patterns where each chemical species is described using a different set of compartments. It is shown that the parameter regions where spatial patterns form are different from the regions obtained by classical deterministic PDE-based models, but they are also different from the results obtained for the stochastic reaction-diffusion models which use a single set of compartments for all chemical species. In particular, it is argued that some previously reported results on the effect of noise on Turing patterns in biological systems need to be reinterpreted.