1. Introduction
A solvent environment is an essential factor one needs to consider when performing molecular dynamics (MD) and structure modeling. With an explicit solvent, molecular dynamics aims to model the movements of individual water molecules predicted by calculations. Although this methodology sounds quite reasonable for modern theoretical approaches, it suffers from considerable and somewhat prohibitive computational cost during implementation. Such obstacles impose limits on the studies of biomolecules, since they are often susceptible to sharp conformational transitions akin to the protein folding process. Therefore, this approach demonstrates rather low performance when calculating the free energies of multiple states of the same multi-atomic systems. Solving this issue requires accounting for the high-volume degrees of freedom demonstrated by the molecules in the solution.
A model system that implemented the solvent implicitly can fulfill the abovementioned task. Briefly, such an approach replaces discrete water molecules with an ideal space continuum, emulated by some combination of dielectric and “hydrophobic” properties, assumed to resemble those of physical water. Not only did this allow us to approximate the general solvent effect by the method of molecular dynamics, but also allows us to assess the pool of solvent degrees of freedom implicitly. The latter is important as soon as this approach allows interactions between the molecules of solute and the solvent to be described with the coordinates of the solute molecules. These interactions are accounted for by computable functions to calculate the solvation free energy that changes on the path of the molecule being transferred from a vacuum with zero dielectric constant (e = 0) into a solvent with moderately high dielectrics (e ≫ 6), including water (e = 80). This solvent model has practical advantages over water represented explicitly and incorporated by the MD simulation. An implicit solvent model demands a much lower computing cost compared to large molecular systems, i.e., when the solution is taken explicitly. Notably, the enhanced sampling of the conformational space is possible because the viscosity of the solvent is omitted. Such an approach can promote the evaluation of solute conformations because solution molecule rearrangements are neglected due to the implicit solvent. Regarding the method of molecular dynamics, the stage of equilibrating the solvent molecules is removed, which is required for MD simulation in explicit water. Finally, the implicit solution model provides a rather instantaneous dielectric response from the solvent, which could be important to monitor flexible structures with labile conformations. Additionally, pattern-averaging the free energy values removes an excessive number of local energy minima that arise from minor shifts in the solvent’s conformations.
Bringing these advantages in effect depends on the kind of software provided by the implicit solvent framework. Generally, the implementation of computationally consumable approximations aims to reveal the electrostatic part of the solvation free energy, which requires further trade-offs between the accuracy and the speed of the modeled dynamics since the stage of molecular docking spliced with MD highly demands algorithmic simplicity. An example of such an approach is the generalized Born (GB) mode [
1]. Compared to the Poisson–Boltzmann (PB) model, the GB method is analytical and is able to approximate the calculated electrostatic part of the solvation free energy. The GB approximation model is relatively simple yet capable of providing an effective way to compute the long-range electrostatic interactions over the molecular structure.
The integral equation method is known in the theory of fluids as the reference interaction site model (RISM) approximation, an approach to studying molecular solvation [
2]. At infinite dilution approximation, the solution density can be neglected and the RISM method enables one to find either solution–solvent or atom-to-atom correlation functions from which the solvation thermodynamics, in particular, the Gibbs energy value can be calculated [
3]. Apart from the atom-to-atom interacting potentials, the “prime load” of the RISM theory is the solvent molecules geometry, which is depicted with a matrix of interatomic distances. The correlation functions and solvation thermodynamics, derived as the RISM equation solutions, correspond to the unique state, i.e., the “rigid” form of the dissolved molecule, including one of the “frozen” states of the MD trajectory that are calculated during the molecular dynamics run. However, the molecules at room temperature often demonstrate a high level of versatility, including during the semi-equilibrium phase of the MD trajectory, where the molecules are still capable of intensively changing their conformations. Fluctuations in molecular geometry influence the solvation thermodynamics and can be assessed by accounting for the atom’s mutual influences, which are defined for each state of geometry, followed by averaging the values estimated for the Gibbs energy. The MD trajectory facilitates a set of the geometry configurations, which are chained by the “frozen” states of the molecular structure that are generated in course of the MD simulation [
4].
In RISM theory, there is an alternative way of considering the conformational versatility of a solute, which is the focus of this study. For the case of infinite dilution, we develop a modified system of RISM equations, which is averaged over the entire MD trajectory only once, RISM with average matrix (RISM-AM). The solutions of this system resemble the solutions of standard correlation functions averaged over the trajectory. Indeed, the proposed method for estimating the mean value of the Gibbs energy from the data on the semi-equilibrium scope of the MD trajectory does not require solving RISM equations for each frame (time point) of the dynamics. However, it is not quite clear a priori whether and how the thermodynamic quantities estimated using the averaged RISM equation correspond to the actual means of thermodynamics.
The Gibbs energy and excess chemical potential are fundamentally important thermodynamic quantities. According to the Ornstein–Zernike theory, the method of integral equations applied to simple liquids allows one to introduce a term for the excess chemical potential, which can be calculated explicitly [
4,
5]. This formula was obtained in general form analytically by performing thermodynamics integration for the energy. Although it has been generalized quite well for the evaluation of polyatomic molecules using RISM [
6], nevertheless, the simpler term for excess chemical potential can be discarded. For this purpose, the nonlinear part of the RISM equation is represented by quite a simple analytical term of the so-called “closure equation” [
4]; otherwise, thermodynamic integration requires numerical calculations.
Unlike in standard RISM equations, the mean matrix allows one to numerically integrate the mean of the entire MD trajectory only once. This approach is substantially simpler than solving the conventional RISM equations for each current time period in an MD experiment, so the free energy integration step is repeated at each time point along the MD trajectory in order to eventually obtain the mean value. Therefore, the new methodology can provide mean value equations with different functions to estimate the Gibbs energy. To not forget is that the closure term takes into account the entropic contribution to the Gibbs energy, which manifests itself in a particular function for correcting the repulsion of the bridge functional added to the equation [
4,
7].
The formula for the excess chemical potential uses a correction term derived from the “first principle” thermodynamic perturbation theory (TPT). Instead of time-consuming numerical integration that requires huge computing, the Gibbs energy can be estimated by evaluating the approximated value according to the TPT.
This functionality of RISM-AM is an advantage of this method, which, in addition, requires further development and full implementation of this approach.
To evaluate the applicability of RISM-AM for estimating the Gibbs energy, the corresponding values are calculated using the perturbation theory formulas by numerical integration and analyzed in comparison with the performance of other methods.
2. Methodology
2.1. RISM Equations for an Infinitely Diluted Solution
The equation for an infinitely diluted solution [
3], when the partial density of the dissolved substance is neglected, relates the direct,
ciα(
r), and full,
hiα(
r), correlation functions, where
r denoted the coordinate, by assigning to the atoms of the solvent molecule (α = 1,…,
M) distributed around the atom of the solute
i (i = 1,…,
N):
and
where
k is the wave vector. In Equation (1), the expression in the parentheses reflects the static susceptibility of the solvent. It involves a matrix
of intramolecular correlation functions describing the rigid bonds between the atoms of the solvent molecule, the matrix
of complete correlation functions of the pure solvent with the density ρ. The upper indexes
v and
u stand for solvent and for solute, respectively, in Equation (1).
The elements of the matrix of intramolecular correlation functions (
Wu) stand for the rigid bonds between the atoms in dissolved molecules [
2]:
where
lij is the distance between the atoms
i and
j, and
δij is the Kronecker delta.
The matrix of intramolecular correlation functions also contains diagonal static susceptibility elements related to the static structure factor. This means that the expression defines a constant and can be either calculated or taken from the experimental data. The dimension of the matrix Wu equals to the number of atoms in the solution molecule, N, and the dimension of the matrix Wv is equal to the number of atoms in the solvent molecule, M. Therefore, the correlation function matrices are the N × M matrices.
In addition to Equation (1), related to the correlation functions in coordinate space to provide the solvability of the system of correlation functions, there is the so-called closure equation, which can be either in the form of a hypernetwork closure (HNC) [
8]:
or as its partially linearized form (PLHNC) [
9]:
While γ
iα(
r) =
hiα(
r) −
ciα(
r) is a nondirect correlation function,
uiα(
r) are pair-wise interaction potentials between the atom
i of the solution molecule and the atom α of the solvent molecule, where
β = 1/
kT with
k being the Boltzmann constant and
T the temperature, is the reciprocal temperature, let us assume:
2.2. Molecular Dynamics and Integral Equations
The RISM equations allow us to obtain correlation functions for a specific geometry of a dissolved molecule, characterized by a given matrix of interatomic distances. Since a molecule in a solution undergoes significant conformational transitions, it would be incorrect to limit the matrix to a single set of interatomic distances. One way to take into account the versatility of a molecule on an equilibrium MD trajectory is to average the values of correlation functions for a set of multiple geometry configurations.
Assume the MD trajectory is a set of frames, each defined by the instant interatomic distances in the dissolved molecule {
lij(
t)}, with
t = 1,…,
NT, where
t stands for the value, which depends on the molecule geometry, and
NT denotes the number of solute configurations on molecular trajectory. Having denoted the static susceptibility of the solvent as V (see Equation (13) in Ref. [
10] for water solvent) and
h = vec(H
uv), and
c = vec(C
uv), the operation
vec(
x) expresses the matrix transformation as a column-vector, where the columns in matrix
x are consequently written one under another; hence, Equation (1) is applicable to every
t configuration in general, such as
where ⊗ is the direct matrix product (Kronecker product) of matrices V and
combines the elements of the matrix of intramolecular correlation functions (see Equation (2)). Hence, this notation of the matrix A
t corresponds to the set of distances {
lij (
t)}. Solutions of Equation (5), assessing the set of
t geometries, allow for the mean of the full correlation functions vector to be determined:
There is a shortcoming in Equation (6): in order to average solutions to the equations, it is necessary to solve the RISM equations for numerous snap-shot configurations of the molecule according to the information provided by the MD trajectory data. This obstacle can be avoided by estimating the mean vector of full correlation functions:
At this point, let us introduce the mean vector of the direct correlation functions formally:
Since the only solution for the mean value of interest is
an average-matrix of the RISM equations could be written as:
We call the solutions to Equation (7) pseudo-average correlation functions. The convenience of Equation (7) is quite straightforward: since a single RISM equation must be solved, the Ā matrix (8) already contains averaged intramolecular correlations for the entire MD trajectory. However, the question remains: what are the differences between the pseudo-average functions and the actual average correlation functions obtained using standard RISM? This difference can be estimated as the difference between the results obtained using Equations (6) and (7). A direct transformation of the difference expression,
yields the expression
assuming the difference between the average values arisen from Equation (7) as the pseudo-mean correlation functions, and the full correlation functions are expressed with the values, averaged on MD trajectory. At the same time, the differences, originated from the use of the linear Equation (7), can be reduced by taking the direct full correlation functions as averaged probe functions. It is worth mentioning that the term for the difference (9) equals to zero if all the matrices A
t are equal to the average matrixes been equal each other. From Equation (9) it follows that the difference Δ is as small as the difference between the intramolecular correlation matrices calculated for the two neighboring points on the MD trajectory. This statement is true for the semi-equilibrium portion of the MD trajectory. Indeed, let us divide the differences arisen in result of averaging intramolecular matrices for the MD trajectory with Equation (9) employed for the particular momentum states on the trajectory:
.
When molecular system subjects to some thermal fluctuations near the equilibrium state, there is no reason to neglect the points, at which the values of ΔĀ
s are small. Therefore, the value of ΔĀ
s can be close to zero, as determined by the direct correlation function:
The linear part of the RISM equation assumes the error is estimated precisely by the expression (7) for the average-matrix. This is a time-consuming task, since estimating the contribution of the nonlinear part of the RISM equation seems somewhat problematic. Therefore, the closure Equations (3) and (4) can be related to the error problem. Anyway, it appears that preliminary calculations are to be made using MD data and to compare the value of Δ with the dispersion of the average correlation function. When the error (Δ) is smaller than the dispersion, the MD data can be used to calculate the mean value by solving the RISM-AM. In such cases, the proposed approach is an applicable tool since the static susceptibility of the solvent matrix V remains constant, while the solute intramolecular correlation matrix Wu can vary. If the functions subjected to significant fluctuations, then the matrix A to be averaged, but the solvent related matrix V does not change.
Finally, the average-matrix of the RISM equations can be defined via the following expression:
2.3. Intramolecular Correlation Matrix
Let us dwell upon the calculation of the mean matrix of intramolecular correlations. Let us consider a simple case of a matrix
defined by interatomic distances in a solute molecule averaged over the trajectory 〈
lij〉 according to the formulation (2) for rigid bonds:
where
Definition (10) has an advantage since the parameter dependence of functions on the scalar value of average interatomic distance. This not only allows the data set to be compressed, but also significantly reduces the storage space for information associated with the data in the w matrix. The matrix memory saving mode solves the main practical problem, since this allows for the numerical implementation of the method. This is important for analysis of large molecules, which typically requires almost prohibitively high computational power. An analysis of the dependence of the parameters (10) shows that the memory saving mode can be achieved using two quite simple techniques. First, inside the molecular globule, there are many approximately equal interatomic distances, since the lengths of covalent bonds are largely limited by the nature of the chemical elements. Once these distances are rounded off to predefined limits, the whole set of interatomic distances can be also reduced to a compact set of the nonequivalent values [
4]. (By rounding these distances to given limits, the entire set of interatomic distances can also be reduced to a compact set of nonequivalent values [
4].) In the next step, intramolecular correlations of distant atoms can be relaxed since the atoms influence each other as a low-force event. As mentioned earlier, correlation relaxation allows for the use of a short form of intramolecular correlation functions, suitable for in silico studies of large molecular systems, including DNA (deoxyribonucleic acid) [
10]. Therefore, calculating the average molecular geometry becomes a task similar to operating an object with rigid geometry. Both cases can be straightforwardly implemented in a software package for use with input data represented by any interatomic distance matrix. Although the expression (10) looks simple, it poorly performs in experiments. Often, the correlation functions did not coincide well enough with the average values, contaminated by some artifacts like the parasitical peaks spaced apart distantly. Thermodynamics deals mostly with the mean values for the cases of rigid small molecules. With all the average values mentioned above, this approach is considered not to be of high use for the analysis of versatile molecular structures.
A different approach can be taken to average intramolecular correlation functions, using histograms of interatomic distance distribution,
, obtained from the corresponding MD trajectory:
Such a way of averaging the values looks more practical for the system of a multistate molecule. This statement is illustrated on
Figure 1, where the upper graph displays the distance distribution between two groups of atoms. The histogram shows that most distances indicate that the atoms are tightly grouped within a 2 Å area. However, in the scope of conformational versatility, there are distances as big as 7 Å, indicated by the corresponding peaks in the histogram. It is understandable that all these wave vector functions are significantly different, with the exception of the numerical volumes around the zero value. Consequently, solutions to equations containing matrices of this type can differ significantly.
Of high practical importance is the issue of using functions averaged with the use of histograms. The averages allow for efficient calculations because the short-range functions in the intramolecular function matrix converge better by solving the RISM equations iteratively. Furthermore, the time required to compute the solutions of the RISM-AM equations for an MD trajectory is comparable to the cost of solving the RISM equations for all the individual (snap-shot) configurations of a molecule. The method (11) does not allow us to determine an equivalent functions to implement full data compression. However, the nature of short-range average intramolecular correlations allows the functions to be cut off at the maximum of the wave vector
k*ij:
while the accuracy power, ξ, of this operation can be set arbitrarily and thereby predetermined. Factor ξ can influence both matrix content and the solution of the equations, since the approximated matrix of ω
ij(
k) elements filled with another approximate values of ω
ij(
k)
H(
k −
k*ij), where
H(
x) states for the Heaviside step function. The choice of the ξ-value magnitude implies that the matrix should be rather a definite one, so the errors of the cutting-off procedure do not exceed the accuracy of the RISM equation solutions, calculated numerically with steady iterations. Additionally, the effect of cutting-off transformation to be reduced even more by setting the boundary for the wave vectors
k*ij of the ω
ij (
k*ij) added correctly to the term (12).
2.4. Thermodynamics
To compare the free energies of the dissolved molecules, one should find the Gibbs free energy, Δ
G, as a sum of the excess chemical potential Δ
μ, and molecular mechanical energy,
EMM:
In case the energy was calculated for a panel of snap-shot configurations, the array of data on configurations has to be averaged:
The excess chemical potential Δ
μ value equals the work performed on the particle being introduced into the medium. On the path, the particle penetrated the medium can be emulated by gradually weakening the interaction potentials with the scaling parameter
λ:
where
g(
r,λ) =
hiα(
r,λ) + 1 determines the atom-to-atom radial functions of distribution of the solvent atom α around the atom
i of the solvent, while the RISM equations solutions are weakened with the
λ parameter considering the inter-atomic interaction potentials,
u(
r,λ). In this vein,
Section 3 deals with the correlation functions, assuming the thermodynamic potentials calculated either with mean, pseudo-mean or ordinary functions.
In some cases, the integration over the
λ-factor in Equation (14) can be performed analytically. Thus, in the theory of Gaussian density fluctuations (GF), the expression for the excessive potential exists in the form [
11]
Another kind of the potential was defined by the Singler–Chandler (SC) expression for the excess chemical potential [
5,
11,
12,
13], and has a different form for the HNC:
As shown [
6], the closure,
treated according to the theory of liquids is related to the term,
where
Bαγ(
r) is a bridge functional, presented by the infinite integrant. It is noteworthy that in RISM the functionality of the bridge is unknown; therefore, the terms for a closure are based either on fluid theory or other assumptions. From above it follows that if the integration of bridge functional γ is possible, the result contains simple functions. For instance, if
Bαγ (γ
′) = 0, then the second part of the term vanishes, so the entire term becomes just the SC expression [
13]. By analogy, the term for the PLHNC kind of closure is defined similarly [
9]:
Additionally, the partial-wave method (PW) can be used to form another expression for the excess chemical potential:
where the elements of matrix Z are:
More elaborate comments on these expressions can be found in the original papers [
14,
15]. The authors introduced the repulsion correction to the closure equation in order to balance the entropy component that could influence the energy [
6]. Therefore, the modified hyper net chain equation corrected with the repulsive bridge correction (RBC) as HNC/RBC takes the form
while the partially linearized hyper net chain equation (PLHNC/RBC) is written as
In addition, the correction index
biα can be defined as
where * denotes convolution in coordinate space, while the repulsive potential
uRBCiμ can be represented by either the repulsive
r−12-term of the Lennard–Jones potential (R12):
or the corresponding part of the Weeks–Chandler–Andersen (WCA) potential:
For the RISM equations with a repulsion correction, the excess chemical potential cannot be revealed in explicit form analytically, though it can be integrated numerically by solving the RISM system for values of λ on the interval from 0 to 1. To date, despite all the progress made in making complex, highly specialized graphics processing unit GPU-empowered computing platforms, integration still remains a laborious, time-consuming, and a costly stage of experiment that tends to consume quite a large amount of computing power, depending progressively on the size of the molecular system under investigation. The advantage is that the term (14) can be modified and then used to average the solutions of Equation (7) along the MD trajectory.
2.5. Calculations
To compare the data assessed with pseudo-mean correlation functions (7) with real mean values (6), we used the following MD trajectories calculated for the peptide hormone oxytocin: peptide in vacuum (A), peptide in explicit water (B), and peptide in implicit solvent according to the GB method (C). The correlation functions were constructed for 1250 peptide configurations selected randomly from each trajectory. These correlation functions were used to assess the mean values of the Gibbs energy calculated using several different terms to find the excess chemical potentials (15)–(18) and marked GF, SC, and PW, respectively. An approximate expression (23), taking into account various expressions for the repulsion corrections (21) and (22), allowed us to find the Gibbs energy, designated as RBC1 and RBC2, respectively. If the data set is representative, it is possible to compare the real mean values with the solutions to the mean equation. Parameter ξ emulated the precision level, which was preset during the cutting-off procedure performed on the intramolecular correlation functions (see Equation (12)). We arbitrarily set the ξ value to 0.01. Histograms of interatomic distance distributions were calculated for 12,500 states taken from the MD trajectory. The mean RISM equation was always solved for the same grid of parameters and the potentials. It is noteworthy that the numerical integration over the parameter λ in Equation (14) with respect to the repulsion corrections (21) and (22) was carried out with an adaptive step. The integration accuracy was constantly controlled by repeating the calculation of the selected steps. The values obtained by the integration are designated as INT1 and INT2, respectively.
4. Discussion
There is small difference between the two approaches of the RISM, either the averaged over the MD trajectory and the RISM-AM, featured with the repulsive correction by the Lennard–Jones potential TPT(R12). It was previously noticed that the total changes in the Gibbs energy values estimated with RISM-AM did not exceed the limit seen due to temperature fluctuations. The ionic peptide complexes provided with an opportunity to gather an information about the objects of different size contained a spare charges. At the same time, the relative differences in the Gibbs free energy assessed by the RISM-AM exceeded the differences observed for the oxytocin. Nevertheless, the main trends (patterns?) of changes in the free energy during the formation of peptide complexes retained similar, no matter what the functionals were used by two versions of the RISM method.
These results support the authors intention to apply the average matrix RISM method to study the thermodynamics and the process of profilament complex formation by the ionic peptides. This method is prospective since allows for thermodynamic integration and analytical consideration of repulsive corrections. The RISM-AM used to estimate the Gibbs free energies of the RADA16, KADA16, and RAEA16 peptide complexes, allowed to assess the free energy values by taking into account the repulsive corrections using numerical thermodynamic integration. Further on, the corresponding data are notified as INT(R12) and INT(WCA).
In addition, the results of continual approaches used to account for the solution were compared with the performance of the RISM. The generalized Born and the solution of the Poison–Boltzmann equation were used to calculate the Gibbs free energy over the same 10 ns trajectories with an integration interval of 2 ps [
20]. The free energy averaged values were obtained by different methods presented in
Table 6. It can be seen the trend, the larger peptide complex, the lower free energy corresponded to the peptide monomer involved in the assembly. Truly remarkable is the difference found between the free energy values determined by GF, PW, and GB or PB approaches. Strangely, the terms of SC, TPT(R12), and TPT(WCA) led to very confusing results contradicted to the internal logic of the process of peptide self-assembly (protofilament formation). At the same time, the results demonstrated by the SC term previously revealed the free energy values were similar to the experimental data obtained for the interaction of DNA with a small dye molecule [
11].
What seems to be important is the method of thermodynamic integration empowered with the repulsion correction, allows to assess the free energy values comparable to the results obtained by the methods based on continual models of solution [
10].
Also, the free energies of interacting parts in the octamer peptide complex were estimated by different methods (
Table 7). The results, supplied with uniform logic of the process of peptides self-assembling, were gathered with GF and PW terms or thermodynamic integration supplied with the repulsion correction by the Lennard–Jones potential—INT(R12). Curiously enough, the usage of the SC, TPT(R12), and TPT(WCA) terms yielded data seeming to display inverted energy trends on the path of growing protofilament formation. The only explanation we have on that fact is that the abovementioned terms for the repulsion corrections may poorly account for the non-polar interactions and thusly undermine the total effects of the interacting charges.
More sensible results obtained with the thermodynamic integration supplied with the Lennard–Jones repulsion correction term (21), are notified as INT(R12) in
Table 7. Importantly, the thermodynamic integration supplied with the WCA repulsion correction part (22) and marked as INT(WCA) in
Table 7, fails to provide the growing protofilament structure via self-assembling process. Since the integration with the repulsion correction by Lennard–Johns term (21) yields better results, the latter points at yet-unidentified features of the WCA term applied to flexible molecular structures with spare charges (22).
5. Conclusions
The results of comparing the RISM-AM and the standard RISM approaches suggest that the mean (7) and averaging method (11) are both applicable to study rigid molecular structures, and ionic peptide complexes with multiple states of freedom as well. Pseudo-mean correlation functions were found coincided well with the actual means at the level of accuracy, typical to the statistical-error-accompanied stage of averaging in the standard RISM.
Consequently, the corresponding thermodynamics obtained with pseudo-average functions differed from the real mean values by no more than the averaging error. At the same time, the costs of computing the pseudo-mean values were found of being approximately less by two orders of magnitude, than the actual means averaged through the intermediate values along the MD trajectory.
Even when the molecular trajectory is non-static, i.e., the molecule undergoes considerable conformational transitions, the RISM-AM method is still applicable, but only after the correct quasi-stationary region located and selected on the calculated MD trajectory. Such a move eliminates the precondition of working with pseudo-mean correlation functions. Generally, the RISM-AM approach appears to be much faster than conventional RISM equations implemented on the same computing platform of hardware.
Since the RISM equation with an average matrix is fast, it can resolve the problem of long lasted session of calculations, when the standard RISM used on moderate hardware. For instance, when using the pseudo-mean correlation functions, the integration of thermodynamics for given solvating temperature dependences was carried out using Big Data arrays. However, it should be kept in mind that the difference between the actual mean and pseudo-mean RISM values for the particular system may be significantly bigger than the dispersion interval.
The question still remains about how consistent these relationships are for both pseudo-means and actual means. This issue is likely to be independently assessed and considered on a case-by-case basis.