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.

Introduction

With increased experimental information on atomic or near-atomic structure of biomolecules and intracellular components, there has been a growing need to incorporate such microscopic data (coming from X-ray crystallography, NMR spectroscopy or cryo-electron microscopy) into dynamical models of intracellular processes. A common approach is to use molecular dynamics (MD) simulations based on classical molecular mechanics. Such MD models are written as relatively large systems of ordinary or stochastic differential equations for the positions and velocities of individual atoms, which can also be subject to algebraic constraints (Leimkuhler and Matthews 2015; Lewars 2016). Although all-atom MD simulations of systems consisting of a million of atoms have been reported in the literature (Tarasova et al. 2017; Farafonov and Nerukh 2019), such simulations are restricted to relatively small computational domains, which are up to tens of nanometres long. It is beyond the reach of state-of-the-art computers to simulate intracellular processes which include transport of molecules over micrometers, because this would require simulations of trillions of atoms (Erban and Chapman 2019).

An example is modelling of calcium () dynamics. On one hand, at the macroscopic level, waves can propagate between cells over hundreds of micrometres and Kang and Othmer (2009) developed a model of waves in a network of astrocytes. It builds on previous modelling work by Kang and Othmer (2007) describing intracellular dynamics as a system of differential equations for concentrations of chemical species involved, including inositol 1,4,5-trisphosphate (), a chemical signal that binds to the receptor to release ions from the endoplasmic reticulum. On the other hand, at the atomic level, Hamada et al. (2017) recently solved -bound and unbound structures of large cytosolic domains of the receptor by X-ray crystallography and clarified the -dependent gating mechanism through a unique leaflet structure.

Although it is not possible to incorporate such a detailed information into modelling by using all-atom MD in the entire intracellular space, there is still potential to design multiscale (multi-resolution) models which compute dynamics with the resolution of individual ions. Dobramysl et al. (2016) implement such a methodology at the Brownian dynamics (BD) level to study puff statistics stemming from receptor channels. Denoting the position of an individual ion by , its diffusive BD trajectory is given bywhere D is the diffusion constant and are three independent Wiener processes. Since individual positions of ions are only needed in the vicinity of channel sites, Dobramysl et al. (2016) model diffusion of ions far away of the channel by a coarser model, utilizing the two-regime method developed by Flegg et al. (2012). This method enables efficient simulations with the BD level of resolution by coarse-graining the BD model in those parts of the simulation domain, where the coarse-grained model can be safely used without introducing significant numerical errors (Flegg et al. 2014, 2015; Robinson et al. 2015).

Although BD models or their multi-resolution extensions simulate individual molecules of chemical species involved, the binding of ions to channel sites or other interactions between molecules are only described using relatively coarse probabilistic approaches. For example, the BD model of Dobramysl et al. (2016) describes interactions in terms of reaction radii and binding probabilities as implemented by Erban and Chapman (2009) and Lipková et al. (2011). Atomic-level information is not included in BD models. In order to use this information, multi-resolution methodologies have to consider MD simulations in parts of the simulation domain. In the case of ions, such a multi-resolution scheme has been developed by Erban (2016), where an all-atom MD model of ions in water is coupled with a stochastic coarse-grained (SCG) description of ions in the rest of the computational domain.

The accuracy and efficiency of such multi-resolution methodologies depend on the quality of the SCG description of the underlying MD model. In this paper, we present and analyze a class of SCG models which can be used to fit non-Gaussian distributions estimated from all-atom MD simulations. While the velocity distribution of the coarse-grained particle can be well approximated by a Gaussian (normal) distribution in our MD simulations, this is not the case of the force distribution. Non-Gaussian force distributions have also been reported by Shin et al. (2010) and Carof et al. (2014) in their MD simulations of particles in Lennard-Jones fluids. Thus our SCG model is formulated in a way which incorporates a Gaussian distribution for the velocity and a non-Gaussian distribution for the force (acceleration).

Given an integer , a coarse-grained particle (for example, an ion) will be described by three-dimensional variables: its position , velocity and 2N auxiliary variables and , where . Denoting , , and , the time evolution of the SCG model is given bywhere is an increasing differentiable function, is its derivative, is its inverse, is a continuous function and are positive constants for and , 2, 3, 4. We note that some of our assumptions on can be relaxed as long as appearing in Eq. (4) can be suitably defined.

The SCG description (2)–(5) includes 2N functions and and 4N additional parameters , which can be all adjusted to fit properties of the detailed all-atom MD model. In particular the SCG model (2)–(5) can better match the MD trajectories of ions than the BD description given by Eq. (1), which only has one parameter, diffusion constant D, to fit to the MD results.

One of the shortcomings of Eq. (1) is that its derivation from the underlying MD model requires us to consider the limit of sufficiently large times. In particular, we need to discretize Eq. (1) with a relatively large time step, say a nanosecond, to use it as a description of the trajectory of an ion. Since the typical time step of an all-atom MD model is a femtosecond, it is difficult to design a multi-resolution scheme which would replace all-atom MD simulations by Eq. (1) in parts of the computational domain. The SCG model (2)–(5) can be used to fit not only the diffusion constant D but other important properties of all-atom MD models, which improves the accuracy of the SCG model at time steps comparable with the MD timestep.

SCG models can be constructed using a relatively automated procedure by postulating that an ion interacts with additional ‘fictitious particles’. Such a methodology has been applied to coarse-grained modelling of biomolecules by Davtyan et al. (2015, 2016) to improve the fit between an MD model and the dynamics on a coarse-grained potential surface. They use fictitious particles with harmonic interactions with coarse-grained degrees of freedom (i.e. they add quadratic terms to the potential function of the system and linear terms to equations of motions) and each fictitious particle is also subject to a friction force and noise. An application of such an approach to ions leads to systems of linear stochastic differential equations (SDEs) and can be used, after some transformation, to obtain a simplified version of the SCG model (2)–(5), where functions and are given as identities, i.e. for and . Using this simplifying assumption in the SCG model (2)–(5), we obtainThis is a linear system of SDEs with 4N parameters. It has been shown by Erban (2016) that such models can fit an increasing number of properties of all-atom MD simulations as we increase N. For example, the linear SCG model (6)–(9) can be used to fit the diffusion constant D and second moments of the velocity and the force for , while the velocity autocorrelation function can better be fitted for larger values of N, e.g. for . However, there are other properties of MD simulations which cannot be captured by linear models even if consider arbitrarily large N. They include, for example, all distributions which are not Gaussian. This motivates the introduction of general functions and in the SCG model (2)–(5).

Considering the SCG model (2)–(5) in its full generality, it can capture more interesting dynamics. However, coarse-grained models can only be useful if they can be easily parametrized. Thus in our analysis, we focus on choices of functions and which both improve the properties of the SCG description and do not complicate its analysis and parametrization. The rest of the paper is organized as follows. In Sect. 2, we consider the linear SCG model (6)–(9) for , which is followed in Sect. 3 with the analysis of the linear model for general values of N. To get some further insights into the properties of this model, we study its connections with the corresponding generalized Langevin equation. In Sect. 4, we consider the nonlinear SCG model (2)–(5) for . We consider specific choices of nonlinearity , for which the model can be solved in terms of incomplete gamma functions. This helps us to design three approaches to parametrize the nonlinear SCG model, which are applied to data obtained from MD simulations. We conclude with the analysis of the nonlinear SCG model (2)–(5) for general values of N in Sect. 5.

Linear model for and the generalized Langevin equation

We begin by considering the linear SCG model (6)–(9) for . To simplify our notation in this section, we will drop some subscripts and denote , , , and for , 2, 3, 4. Then Eqs. (6)–(9) read as followswhere X is (one coordinate of) the position of the coarse-grained particle V is its velocity, U is its acceleration, Z is an auxiliary variable, is white noise and , , 2, 3, 4, are positive parameters. In order to find the values of four parameters suitable for modelling ions, Erban (2016) estimates the diffusion constants D and three second moments , and from all-atom MD simulations of ions ( and ) in aqueous solutions. The four parameters of the SCG model (10)–(13) can then be chosen asThen the SCG model (10)–(13) gives the same values of D, , and as obtained in all-atom MD simulations.

Since the model (10)–(13) only has four parameters, we can only hope to get the exact match of four quantities estimated from all-atom MD. To get some insights into what we are missing, we will derive the corresponding generalized Langevin equation and study its consequences. The generalized Langevin equation can be written in the formwhere is a memory kernel and random term R(t) satisfies the generalized fluctuation-dissipation theorem, given below in Eq. (21). To derive the generalized Langevin equation (15), consider the two-variable subsystem (12)–(13) of the SCG model. Denoting where stands for transpose, Eqs. (12)–(13) can be written in vector notation as followswhere matrix and vectors , are given asLet us denote the eigenvalues and eigenvectors of B as and respectively. The eigenvalues of B are the solutions of the characteristic polynomial They are given bySince and are positive parameters,  we conclude that real parts of both eigenvalues are negative. In what follows, we will assume . Then we have two distinct eigenvalues and the general solution of the SDE system (16) can be written as followswhere is a constant vector determined by initial conditions and matrix is given asi.e. each column is a solution of the ODE system . Calculating the inverse of and considering long-time behaviour, Eq. (18) simplifies towhere memory kernel is given byand noise term R(t) is Gaussian with zero mean and the equilibrium correlation function satisfying the generalized fluctuation-dissipation theorem in the formUsing (17), memory kernel (20) can be rewritten aswhere . We note that the auxiliary coefficient is a square root of a real negative number for . However, formula (22) is still valid in this case: for it can be rewritten in terms of sine and cosine functions, taking into account that is pure imaginary, and

The memory kernel , given by Eq. (22), is plotted in Fig. 1a for different values of parameter . For simplicity, we use non-dimensionalized versions of our equations with dimensionless parameters and . We choose three different values of so that the values of are 1, and . In Fig. 1b, we plot the equilibrium velocity autocorrelation function which is defined asfor . More precisely, we plot which is normalized so that its value at is equal to 1. It is related to the memory kernel bywhere is the Laplace transform of the memory kernel and denotes Laplace inversion. Following Erban and Chapman (2019), we evaluate the right hand side of Eq. (23) as follows. Substituting Eq. (22) into (23), we obtainThe polynomial in the denominator, has positive coefficients. Since , it has one negative real root in interval , which we denote by . The other two roots ( and say) may be real or complex, but if they are complex they will be complex conjugates since p(s) has real coefficients. Assuming that the real part of each root is negative, we first find the partial fraction decomposition of the rational function in (24) aswhere are constants (which depend on , and ). Then we can rewrite Eq. (23) asThe results computed by (25) are shown in Fig. 1b. We note that although Eq. (25) may include complex exponentials, the resulting is always real. Since the diffusion constant, D, and the second moment of the equilibrium velocity distribution, , are related to bythe parametrization (14) guarantees that both the value of and the integral of are captured accurately. However, the simplified SCG description (10)–(13) is not suitable to perfectly fit the velocity autocorrelation function or the memory kernel for all values of . In order to do this, we have to consider the SCG model (6)–(9) for larger values of N as it is done in the following section.

Fig. 1

a Memory kernel given by Eq. (22) for , and three different values of , namely (solid line, ), (dashed line, ) and (dot-dashed line, ). b Normalized velocity autocorrelation function computed by using Eq. (25) for the same parameter values as in panel (a)

General linear SCG model and autocorrelation functions

Considering the linear SCG model (6)–(9) for general values of N, we can solve Eqs. (8)–(9) for each value of , 2, , N to generalize our previous result (19) aswhere kernel is given by [compare with (22)]withand noise term is Gaussian with zero mean and the equilibrium correlation function satisfyingSubstituting (26) to (7), we obtain the generalized Langevin equationwhereIn particular, we have 3N parameters to fit memory kernel , which can be estimated from all-atom MD simulations. There have been a number of approaches developed in the literature to estimate the memory kernel from MD simulations. Shin et al. (2010) use an integral equation with relates memory kernel with the autocorrelation function for the force and the correlation function between the force and the velocity. Estimating these correlation functions from long time MD simulations and solving the integral equation, they obtain memory kernel . Other methods to estimate the memory kernel, , of the corresponding generalized Langevin equation (29) have been presented by Gottwald et al. (2015) and Jung et al. (2017).

An alternative approach to parametrize the linear SCG model (6)–(9) is to estimate the velocity autocorrelation function, , from all-atom MD simulations. This can be done by computing how correlated is the current velocity (at time t) with velocity at previous times. Since Eqs. (10)–(13) are linear SDEs, we can follow Mao (2007) to solve them analytically, using eigenvalues and eigenvectors of matrices appearing in their corresponding matrix formulation. Using this analytic solution, Erban (2016) use an acceptance-rejection algorithm to fit the parameters of linear SCG model (6)–(9) for to match the velocity autocorrelation functions of ions estimated from all-atom MD simulations of and in the SPC/E water.

Since the parameter given by (28) is a square root of a real number, it can be both positive or purely imaginary. In particular, kernels given by Eq. (27) can include both exponential, sine and cosine functions as illustrated in Fig. 1a. Since memory kernel is given as the sum of in Eq. (30), typical memory kernels and correlation functions estimated from all-atom MD simulations can be successfully matched by linear SCG models for relatively small values of N. However, as shown by Mao (2007), analytic solutions of linear SDEs also imply that the process is Gaussian at any time , provided that we start with deterministic initial conditions. Thus the linear SCG model (6)–(9) for abtitrary values of N can only fit distributions which are Gaussian. This motivates our investigation of the nonlinear SCG model in the next two sections.

Nonlinear SCG model for

We begin by considering the nonlinear SCG model (2)–(5) for . As in Sect. 2, we simplify our notation by dropping some subscripts and denoting , , , , , and for , 2, 3, 4. Then Eqs. (2)–(5) read as followswhere X denotes (one coordinate of) the position of the coarse-grained particle, V is its velocity, U is its acceleration, Z is an auxiliary variable, is white noise, , for , 2, 3, 4, are positive parameters and functions and are yet to be specified.

Equation (31) describes the time evolution of the position, while Eqs. (32)–(34) admit a stationary distribution. We denote it by p(v, u, z). Then gives the probability that , and at equilibrium. The stationary distribution, p(v, u, z), of SDEs (32)–(34) can be obtained by solving the corresponding stationary Fokker-Planck equationwhich giveswhere C is the normalization constant, and functions G and H are integrals of functions g and h, respectively, which are given byWe note that for the special case where g and h are given as identities, i.e. for , the nonlinear SCG model (31)–(34) is equal to the linear SCG model (10)–(13) and functions G and H are . Then the stationary distribution (35) is product of Gaussian distributions in v, u and z variables. In particular, we can easily calculate the second moments of these distributions in terms of parameters . Estimating these moments from all-atom MD simulations, we can parametrize the resulting linear SCG model (10)–(13) as shown in Eq. (14). However, if we want to match a non-Gaussian force distribution, we have to consider nonlinear models. A simple one-parameter example is studied in the next section.

One-parameter nonlinear function

Consider that g is a function depending on one additional positive parameter as followswhere we use to denote the sign (signum) functionThe function defined by (37) only satisfies our assumptions on g for as it is not differentiable at for , but we will proceed with our analysis for any positive . Consider that function h is an identity, i.e. for , then Eqs. (31)–(34) reduce towhere we would have to be careful, if we used this model to numerically simulate trajectories for , because of possible division by zero for in Eq. (41). If , then we do not have such technical issues. Using Eq. (35), the stationary distribution is equal towhere the normalization constant C is given byIntegrating (43), we getwhere is the gamma function defined asLet . Integrating (43), we getUsing (45) for and , we obtain the following expression for kurtosisIn particular, the kurtosis is only a function of one parameter, . It is plotted in Fig. 2a as the blue solid line, together with the kurtosis obtained for a more general two-parameter SCG model studied in Sect. 4.2. We observe that the distribution of U is leptokurtic for and platykurtic for . If is equal to 1, then our SCG model given by Eqs. (31)–(34) reduces to the linear SCG model given by Eqs. (10)–(13), i.e. the stationary distribution is Gaussian and its kurtosis is 3. This is shown by the dotted line in Fig. 2a.

Fig. 2

a Kurtosis given by Eq. (59) as a function of parameter for three different values of parameter . The result for (blue solid line) corresponds to the case of one-parameter function g, defined by (37), where the kurtosis is given by (46). b Distribution of U estimated from a long-time MD simulation (blue circles) compared with the results obtained by the linear SCG model (10)–(13) (black dotted line), nonlinear SCG models (31)–(34) with one-parameter function g, defined by (37), fitting and (red dot-dashed line) and and (green dashed line), and the nonlinear SCG model (31)–(34) with two-parameter function g defined by (52), matching all three moments , and (cyan solid line) (color figure online)

Since Eq. (46) only depends on parameter , we can use the kurtosis of the acceleration distribution (which is equal to the kurtosis of the force distribution) esimated from MD simulations to find the value of parameter To calculate the kurtosis, we estimate the fourth moment in addition to the second moment, , used before in our estimating proceduce (14) for the linear model. In particular, we not only get Eq. (46) for calculating the value of parameter , but also a restriction on other parameters , and . Using (45) for , it can be stated as followswhere we have used properties of the gamma function, including and Euler’s reflection formula, , to simplify the right hand side. We note that in the Gaussian case, , the right hand side of Eq. (47) further simplifies towhich is indeed the formula for the second moment of U given by the linear SCG model (10)–(13). Equation (47) provides one restriction on four remaining parameters, , and , which need to be specified. This can be done by estimating three additional statistics from MD simulations, as in the case of the linear SCG model (10)–(13) in Eq. (14). Indeed, the stationary distributions of V and Z are Gaussian with mean zero. Their second moments and the diffusion constant, D, for the nonlinear SCG model (31)–(34) can be calculuted asTherefore, assuming that D, , are obtained from MD simulations and is given by (47), we can calculate parameters byWe note that in the Gaussian case, , we can substitute Eq. (48) for and the parametrization approach (50)–(51) simplifies to Eq. (14) used in the case of the linear SCG model (10)–(13). In the next subsection, we generalize formula (37) to a two-parameter function and show that the parametrization approach (50)–(51) is still applicable to the case of more general SCG models.

Two-parameter nonlinear function

Consider that g is a function depending on two positive parameters and as followswhere function is defined by (38). In particular, our expression for function g is equal to the formula (37) for sufficiently large values of |y|. As discussed in the previous section, if we used formula (37), there would be some issues for y close to zero [for example, the division by zero for and in Eq. (41)], so our generalized formula (52) replaces (37) with a linear function for smaller values of |y|. On the face of it, it looks that there could also be some issues with the generalized formula (52), because it is not strictly increasing for . However, function (52) is increasing and invertible away of this region with its inverse given byMoreover, what we really need in Eqs. (31)–(34) is which can be defined as the following continuous functionwhere the removable discontinuity at has disappeared because we have defined . Integrating (52) and substituting (53), we getwhere G is the integral of function g defined by (36). Consider again that h is an identity, i.e. for . Then the stationary distribution (35) is again Gaussian in V and Z variables with their second moments given by Eq. (49). Let us denote the marginal stationary distribution of U byUsing (35) and (54), we havewhere is the normalization constant given byLet us defineIntegrating (55), we get, for any ,where function is defined byand (resp. ) is the upper (resp. lower) incomplete gamma function defined bySubstituting and in Eq. (57), we getThis formula for the kurtosis is visualized in Fig. 2a as a function of parameter for three different values of parameter . We note that the case corresponds to the case studied in Sect. 4.1. If , then Eq. (56) implies . Since and , where is the standard gamma function given by (44), we can confirm that Eq. (59) converges to our previous result (46) as

Substituting into (58), we obtain Consequently, using in Eq. (57), we obtainUsing in Eq. (57), we getConsequently, if we use MD simulations to estimate not only the second and fourth moments, and , but also the first absolute moment , we can substitute the estimated MD values into Eqs. (59) and (61) to obtain two equations for two unknowns and . Solving these two equations numerically, we can get and . Then we can use (56) and (60) to get the original parameters and byMoreover, Eq. (56) also implies the following restriction on other parameters , and This restriction is equivalent to restriction (47). Therefore, assuming again that D, , are obtained from MD simulations and is given by (63), we can calculate parameters , , and by Eqs. (50)–(51).

We note that the two additional parameters and can be used to satisfy both Eqs. (59) and (61), while in Sect. 4.1 we could only use one equation (Eq. (46) for kurtosis) to fit one parameter . However, in the case of one-parameter function (37), we could (instead of fitting the kurtosis) match the quantity with MD simulations, i.e. we could replace Eq. (46) by Eq. (61) simplified to the one-parameter case corresponding to function (37). Passing to the limit in Eq. (61) and using Euler’s reflection formula, , we obtain that the one-parameter nonlinearity (37) implies the following formulaThus, in Sect. 4.1, we could use and estimated from long-time MD simulations to calculate the left hand side of Eq. (64), which could then be used to select parameter . Other parameters could again be chosen by Eqs. (50)–(51).

Application to MD simulations

In Sects. 4.1 and 4.2, we have presented three approaches to fit nonlinear SCG models which have non-Gaussian force distributions to data obtained from MD simulations. In this section, we apply them to the results obtained by an illustrative MD simulation of a Lennard-Jones fluid, where we consider a box of 512 atoms which interact with each other through the Lennard-Jones force terms for parameters given for liquid argon (Rahman 1964), i.e. particles interact in pairs according to the Lennard-Jones potential , where , and r being the distance between particles. We use standard NVT simulations where the temperature () is controlled using the thermostat of Nosé (1984) and Hoover (1985) and the number of particles ( in a cubic box of side ) is kept constant by implementing periodic boundary conditions.

Using a long time MD simulation (time series of lentgh ), we estimate three moments , and as averages over all three coordinates, i.e.where is the acceleration of one specific atom (tagged particle) to which our SCG model is applied. Rounding all computational results to three significant figures, we obtain , and .

In Fig. 2b, we plot the equilibrium MD distribution of the acceleration (average over all three coordinates) using blue circles. The resulting distribution is leptokurtic (with positive excess kurtosis). Its kurtosis has been estimated as . The numerical values on the u-axis in Fig. 2b are expressed in [nm ps]. Since the acceleration, U, is proportional to the force exerted on the tagged particle (with the scaling factor equal to the atomic mass of argon), the plot of the acceleration distribution in Fig. 2b can also be interpreted as the plot of the force distribution, which has the same kurtosis, provided that we suitably rescale the units on the u-axis.

If we only attempt to fit the value of , we could parametrize the linear SCG model (10)–(13), which leads to the Gaussian acceleration distribution (plotted as the black dotted line in Fig. 2b). Using the one-parameter nonlinear function (37) from Sect. 4.1, we can use Eq. (46) to find parameter so that the nonlinear SCG model gives the same kurtosis as observed in MD simulations (). The resulting distribution is given as the red dot-dashed line in Fig. 2b. It matches both second and fourth moments, and .

Using all-atom MD simulations, we can not only estimate the kurtosis, but other dimensionless ratios of moments of U. For example, we obtain . This estimate can be substituted in Eq. (64), which provides an alternative approach to obtain the value of parameter of the one-parameter nonlinear function (37). Using and solving Eq. (64) numerically, we obtain The resulting distribution, which matches and , is plotted as the green dashed line in Fig. 2b. We note that the parameter is dimensionless, because both Eqs. (46) and (64) only depend on dimensionless quantities estimated from MD simulations. Since both distributions (for and ) are given by (43), they are unbounded for u close to zero. This motivates the choice of our two-parameter nonlinear function g used in Sect. 4.2.

Substituting and in Eqs. (59) and (61) and solving them numerically, we obtain and . Substituting into (62), we get the two parameters of model from Sect. 4.2 as and . The resulting distribution, given by Eq. (55), is plotted in Fig. 2b as the cyan solid line. We observe that the distribution is now bounded. It is a piecewise defined function which is Gaussian for the values of u satisfying , which removes the singularity at . At the same time, the distribution given by Eq. (55) matches all three moments estimated from MD simulations, , and . As we can see in Fig. 2b, this distribution does not perfectly fit the acceleration distribution estimated from MD simulations. If our aim is to obtain a SCG model which better fits the whole distribution, we can use SCG models for larger values of N as we will discuss in the next section.

Nolinear SCG model for general values of N

We have already observed in Sects. 2 and 3 that the linear SCG model (6)–(9) can match the MD values of a few moments for , while we need to consider larger values of N to match the entire velocity autocorrelation function. Considering the nonlinear SCG model (2)–(5), we have two options to capture more details of the non-Gaussian force distribution observed in MD simulations. We could either keep , as in Sect. 4, and introduce additional parameters into nonlinearity , or we could consider larger values of N. In Sect. 4, we have shown that by going from one-parameter to two-parameter function g, we improve the match with MD results. In this section, we will discuss the second option: we will use larger values of N.

Consider equations corresponding to the i-coordinate, , 2, 3, of the nonlinear SCG model (2)–(5). Let us denote the stationary distribution of Eqs. (3)–(5) byThen gives the probability that , and , for , 2, , N, at equilibrium. The stationary distribution can be obtained by solving the corresponding stationary Fokker-Planck equationOur analysis in Sect. 4.1 shows that parameters , and appear on the left hand side of Eq. (47) as a suitable fraction, which in the Gaussian case corresponds to the second moment of the acceleration (see Eq. (48)). Considering general N, we define this fraction as new parametersand we again assume that the second moment of the velocity distribution, , can be estimated from long-time MD simulations. In order to find the stationary distribution, we will require that parameters , , and satisfy (compare with Eq. (49) for )Then the stationary distribution, obtained by solving (65), is given bywhere C is the normalization constant and functions and are integrals of functions and , respectively, which are given byFollowing (37), we assume that and each is a function of one additional positive parameter , 2, , N, given asThen we haveThen the stationary distribution (66) is Gaussian in and variables and we can integrate (66) to calculate the marginal distribution of byConsequently,where the normalization constant is given byIntegrating (68), we can calculateasThe acceleration of the coarse-grained particle is given byUsing the symmetry of (68), odd moments of are equal to zero. In particular, and for , 2, , N. Consequently,which givesSubstituting Eq. (69) for moments on the right hand side of Eq. (72), we can express the kurtosis of in terms of 2N parameters and , where , 2, , N. For example, if we choose the values of dimensionless parameters equal to given numbers and define new parametersthen Eq. (69) implies that is a linear function of and is a quadratic function of . Equations (70) and (71) can then be rewritten as the following system of two equations for where and are known constants, which will depend on our initial choice of values of . Thus, using , we still have an opportunity to not only fit the second and fourth moments of the force distribution, but other moments as well. For example, the 6-th moment, , would include the linear combination of the third powers of . We could also fit other properties of the force distribution estimated from MD simulations. For example, we could generalize one-parameter nonlinearities (67) to two-parameter nonlinear functions, as we did in Eq. (52). Then we could match the value of the distribution at , if our aim was to get a better fit of the MD acceleration distribution obtained in the illustrative example in Fig. 2b. Another possible generalization is to consider nonlinear functions , provided that we estimate more statistics on the auxiliary variable Z from MD simulations.

Discussion and conclusions

We have presented and analyzed a family of SCG models given by Eqs. (2)–(5), which can be parametrized to fit properties of detailed all-atom MD models. A special choice of functions and in Eqs. (2)–(5) leads to the linear SCG model (6)–(9) which is used in a multiscale (multi-resolution) method developed by Erban (2016) as an intermediate description between all-atom MD simulations and BD models. The linear SCG model is studied in more detail in Sects. 2 and 3, where we highlight that 4N parameters of this model can match some statistics estimated from all-atom MD simulations with increased accuracy as we increase N, but there are also statistics which cannot be matched for any value of N. They include non-Gaussian force distributions.

In Sects. 2 and 3, we show that the linear SCG model (6)–(9) corresponds to the generalized Langevin equation with the stochastic driving force being Gaussian. Such systems have been analysed since the work of Kubo (1966). One approach to match non-Gaussian MD force distributions could be to use the non-Gaussian generalized Langevin equation which was analyzed by Fox (1977) using methods of multiplicative stochastic processes. However, if we want to generalize the linear SCG model (6)–(9) while keeping its structure as a relatively low-dimensional system of SDEs, then it can be done by introducing nonlinear functions and as shown in Eqs. (2)–(5). The advantage of the presented approach is that we can directly replace the linear model by Eqs. (2)–(5) in multiscale methods which use all-atom MD simulations in parts of the computational domain and (less detailed) BD simulations in the remainder of the domain. Coupling MD and BD models is a possible approach to incorporate atomic-level information into models of intracellular processses which include transport of molecules between different parts of the cell (Erban 2014, 2016; Gunaratne et al. 2019).

The nonlinear SCG model (2)–(5) is studied in Sect. 4 for . Describing the nonlinearity as the one-parameter function given by (37), we can use its dimensionless parameter to match the kurtosis of the force distribution estimated from all-atom MD simulations. Although the one-parameter case is easy to analyze in terms of the gamma function, it has some undesirable properties for small forces. If , we can obtain large terms in the dynamical equation (41) for small values of U; this corresponds to the zero value of stationary probability distribution (43) for . If , we have small terms in the dynamical equation (41), but the stationary probability distribution (43) is unbounded for . In Sect. 4.2, we have shown that these issues can be avoided if the two-parameter nonlinear function (52) is used instead of the one-parameter function (37). The resulting equations are solved in terms of incomplete gamma functions. In Sect. 5, we study the nonlinear model for general values of N where each is a one-parameter nonlinearity given by Eq. (67). However, we could also consider two-parameter functions , like we did in Eq. (52) for , to improve the properties of the SCG model for general values of N.