A Bayes-inspired theory for optimally building an efficient coarse-grained folding force field

2.1. Setting up the system

In the IsRNA model, we use two CG beads (P and S) to represent the phosphate and sugar moiety on the backbone, respectively. For the base, purines (pyrimidines) have three (two) CG beads, and the model uses a total of ten beads. Non-local, base pairing interactions are realized with two pairwise interactions, so we can model hydrogen bonding configurations and energies for various base pairing configurations. See Fig. 1 for details of the mapping from all-atom structure to CG representation.

We begin our approach with the most basic assumptions about RNA polymer interactions. Namely, we assume a non-interacting reference state and extract statistical potentials of bond stretching connectivity and Lennard-Jones (LJ) volume exclusion terms so that the backbone energy is

and we use those as the starting CG force field, generating the ref₀ conformational ensemble. As is shown in Fig 2A, this information is sufficient to reproduce the θ₂ bond angle distribution through simulation, which is already determined by the chain connectivity and volume exclusion information.

Our goal is to accurately estimate the folding energy function for RNA interactions, which are captured by bond distance b, bond angle θ, torsion ϕ, and distance-based r non-local base interaction energies. In summary, our total energy function can be written as

Although this is written as a sum of energy functions, the correlations are implicitly included through simulations, which allows us to ascertain what information about a new folding coordinate is already included in the estimated CG force field.

2.2. Supporting the IsRNA theory with statistical mechanics

Here, we provide an example to illustrate how we can use statistical mechanics to build the total force field by consecutively updating our prior energy function to enforce compatibility with an observed folding coordinate distribution. Using $E_{\text{ref }}^{(i)}$ as the reference state force field from inclusion of the ith folding coordinate, we can run a CGMD simulation and extract the simulated distribution p_ref(θ_j) = p_sim(θ_j)|_ref of the jth bond angle θ_j. From the PDB database, we also know the observed distribution p_obs(θ_j), and we can calculate the difference between the observed and reference energy functions using

where we use the dimensionless k_BT = 1 in this section for simplicity. Accounting for the difference between the observed distribution and the CGMD-generated reference distribution, we define the energy function for the next simulation as

Physically, this ensures that our energy function is only adding the missing information about θ_j into the force field. Explicitly, we can show that this energy function will reproduce the observed distribution. The marginal reference distribution of θ_j generated via simulation can be written in terms of the reference joint energy function with previously added terms integrated out:

Similarly, the simulated distribution of θ after including information about the observed distribution can be written

which shows that simulations using $E_{\text{ref}}^{(i+1)}(b,r,\theta ,\varphi )$ as the CGMD force field can reproduce the marginal probability for some bond angle θ_j after information about that coordinate is included in the force field. Indeed, once we include a coordinate, we find that the simulated distribution does match the observed distribution for that folding coordinate (see Fig. 2A–D).

To emphasize our previous point, although the force field is an additive sum of the different terms, the local variables (b, θ, ϕ) are still correlated through the non-local distance (r), which is a function of (b, θ, ϕ) in a highly coupled form. Similarly, the E(θ) and E(ϕ) energy terms are also coupled through non-local interactions. Without inclusion of the non-local terms, some of the local terms may be physically uncorrelated.

2.3. Generalizing and numerically stabilizing the model

Up until here, we have simplified the theoretical arguments in the previous work [18]. Now, we will pivot to expand the IsRNA approach in greater theoretical and practical detail to provide a general method for parameterizing a CG polymer force field. In general, the joint probability distribution cannot be uniquely determined from the set of marginal probability distributions that compose it, which is why we cannot simply multiply together all of the observed marginal distributions we find in the structure database to calculate the joint distribution and the total force field. However, by simulating the conditional probability distributions, we may estimate the joint probability distribution. In theory, we can broaden our method to use any collective variables of the system, and we are not limited to simple folding coordinates, such as distances and angles. However, the more time-consuming a collective variable is to calculate, the less efficient simulations using that variable will be, so simple distances and angles are appealing building blocks for a CG force field. We take the general perspective that the marginal distribution of a collective variable extracted from a reference simulation as described in the previous section can be viewed as the marginal distribution, conditional on the previously included distributions of collective variables. Then, the estimated energy function after each consecutive step is extracted from the estimated joint distribution of all the considered collective variables, conditional on the already included distributions of collective variables.

Our energy function is considered to converge when the estimated joint distribution after a step is sufficient to reproduce via simulation all of the considered observed marginal distributions, including those marginal distributions that have not been explicitly included in the joint distribution. For instance, assume we have a set of m calculated collective variables (x₀, … , x_m). Assume that so far in our procedure, we have explicitly included j (< m) collective variables to produce an energy function E(x₀, … , x_j). If the marginal distributions observed from the structural database agree with the marginal distributions extracted from the simulation results, the energy function converged, and our force field is adequately parameterized, having captured all of the available information in the considered distributions. While here we use marginal probability distributions, the explicit inclusion of some correlation effects by consecutively including some joint distributions of pairs of collective variables may be beneficial to explore, if enough data is available.

Again using statistical mechanics, we show here that the model provides a sound method for parameterizing a predictive force field. After running the simulation with reference energy parameters for backbone bond length and volume exclusion to provide a starting reference state with E = 0, we obtain the simulated distribution for the first collective variable to be added, p_sim(x₀)|₀. We can write the energy function function from the observed structures as

where p(x₀) is observed from solved structures. Again, we use k_BT = 1 for convenience. This energy function provides constraints for the next simulation step and is analytically sufficient to reproduce the p(x₀) distribution via simulation.

which is the result we obtained in Eq. 6. Here, we used the identity g(x₀) = Ω · p_sim(x₀)|₀, where Ω is the total number of conformations generated by the simulation, and g(x₀) is the degenerate distribution function obtained from simulation.

Generalizing the estimated joint distribution after consecutive inclusion of j marginal distributions, we find

where j = n + 1, P (x_j) = p_est(x₀, … , x_j) are the joint distributions estimated from inclusion of j = n + 1 marginals, and E_n = E_est(x₀, … , x_n) is the estimated total energy function that provides constraints in the simulation step. The question now is whether we can analytically show that this estimation of the joint probability will provide appropriate energy constraints for simulated reproduction of the marginal probabilities observed in the database. To evaluate this, we consider the case where j = 1 (i.e. the inclusion of the second observed distribution into the energy function). Generating a simulated distribution of x₁ on the basis of constraints provided by Eq. 7, ${p_{\text{sim}}\left(x_1\right)|}_{E_0}$, we include the observed marginal distribution p(x₁) using

where ideally

everywhere, according to Eqs. 6 and 8. However, retaining this in our expression allows us to correct any small numerical errors that would grow with propagation as we build our estimated joint distribution, maintaining compatibility with the previously incorporated marginal distributions at each step. This can also be viewed as including effects from cross conditional probabilities into the estimated joint distribution. At this stage, our joint energy function has the form

where as we alluded to earlier, the cross-conditional correction to the energy ${E\left(x_0\right)|}_{x_0,x_1}\approx 0$ everywhere, which implies that in the next step ${p_{\text{sim}}\left(x_0\right)|}_{E_1}=p\left(x_0\right)$. Since we have explicitly shown our decomposition of the joint into quasi-marginals that account for cross conditionals, the argument we made in Eqs. 6 and 8 holds, and it can be shown that the simulated distribution of x₀ remains true to the observed marginal distribution from the structure database(see Figs. 3 and 4 for examples).

2.4. Maximizing entropy to optimize model construction

Although the consecutive procedure detailed above will lead to an estimated joint distribution of our collective variables that is consistent with all of the extracted marginal distributions—regardless of the order that the collective variables are included—the convergence and dimensionality of the energy function can be optimized on the set of all considered collective variables using the principle of maximum entropy, capturing all of the information with the fewest possible variables. In order to speed up convergence of the energy function, the IsRNA model incorporates the collective variable that requires the most change to the reference energy function in order to reproduce the observed distribution of the variable via simulation. This well-defined approach leads to an intuitively reasonable order of collective variable inclusion. First, we extract bond length and volume exclusion terms. After that, we include increasingly non-local terms, starting with bond angles, moving to torsion angles, and ending with pairwise distances. Viewing our problem from a Bayesian perspective, our objective is to optimally update from a prior distribution $Q={q_{\text{est}}\left(x_j\right)|}_{E_n}$ to a posterior distribution $P={P_{\text{est }}\left(x_j\right)|}_{E_j}$, given that the posterior lies within a certain family of possible distributions, $p={p_{\text{est }}\left(x_0,\dots ,x_j\right)|}_{E_j}$. In our case, the posterior depends on which collective variable we choose to include next (x_j), and we want to include the collective variable that maximizes our information gain at each step. The selected posterior P is that which maximizes the entropy S[p, q]. Since prior information is valuable, the functional S[p, q] is chosen so that beliefs are minimally updated to the extent required by the new information.

The relative entropy between a prior distribution q(x) and a posterior distribution p(x) is given by

This is a measure of surprise that is also called the Kullback-Leibler divergence or the information gain. For our purposes, the term “information gain” helps us intuitively understand how we can use this concept to parameterize a CG force field. Including this concept in our iterative, consecutive procedure will allow us to reach a convergent energy function with the fewest marginal distributions included in our calculation, which has three main benefits. First, using fewer collective variables to parameterize the force field will accelerate simulations. Second, parameterizing the force field in this fashion requires less computational resources. Third, the concept of maximum entropy is well established, and the direct connection between IsRNA and the theory provides support for the validity of IsRNA. In essence, we check after each simulation to see which candidate observed marginal distribution disagrees the most with the simulated marginal distribution that was constrained by the energy function comprised of all of the previously incorporated marginal and cross conditional distributions. The one that disagrees the most gets incorporated into the joint distribution, providing the minimum information that the model must include to gain compatibility with all of the considered marginal distributions.

As can be seen in Eq. 12, our energy constraints supplied by new information at step j can be written as

where j = n + 1. Identifying the candidate p(x_j) as p and the previously known information about the distribution of x_j, ${p_{\text{sim}}\left(x_j\right)|}_{E_n}$, as q, we see that the principle of maximum entropy agrees with the incorporation of x_j into our energy function as long as no other x_j can be found that increases the amount of information we can gain in this step.

More directly, we can take the ratio of the estimated joint probability before and after incorporating the new information to see how this is true. Using Eq. 9

the joint probability from the previous step cancels in the ratio and leaves us with our proposed information gain times the product of cross conditional probability corrections, which we can take for the moment to be approximately unity (i. e. not contributing much to the information gain, according to Eqs. 6 and 8). We therefore find

which is the main ratio contributing to the maximum entropy (information gain) formula.

Although we may choose any order of collective variables to include in our energy function (see Figs. 2C–D), this method provides us with a method to minimize the number of variables required to accurately simulate structures in agreement with the solved structure ensemble. The principle of maximum entropy method also optimizes the parameterization process—removing the need to guess and check the parameterization order—and supports our choice of marginal energy function. As we can see above, the principle of maximum entropy is consistent with our statistical mechanics approach to update the estimated joint probability distribution and corresponding energy function.

2.5. Applying the model to fold RNA

Finally, the derived CG force field was applied to study RNA folding behaviors and 3D structure predictions. From coiled states, several RNA molecules were de novo folded by applying the derived IsRNA force field through simulated annealing Molecular Dynamics (MD) in the modified LAMMPS [31] package. In detail, annealing MD simulations were carried out on the structures, scaling the temperature from 350K to 200K over 500 ns of simulation time with an integration time step of Δt=1fs. Then, an additional 500 ns of simulation (1μs simulation time in total for each RNA molecule) was performed at constant temperature (200K) to sufficiently sample the conformational space and snapshots were recorded every 100 ps from the trajectories. To obtain the predicted 3D structure, structures in the 10% lowest potential energy within the collected snapshots were clustered and the centroid structure of the largest cluster was chosen. As shown in Fig. 5, for the 19-nucleotide duplex with a cytosine bulge [32] (PDB ID: 1DQF) and the 12-nucleotide CUUG hairpin loop [33] (PDB ID: 1RNG), the 3D structures predicted by the IsRNA model have heavy-atom root-mean-square deviations (RMSDs) of 2.2 Å and 2.5 Å, respectively. Additionally, the folding behaviors of these RNA molecules are characterized by funnel-like plots of the RMSDs with potential energies. Furthermore, simulation in the proposed IsRNA model successfully elucidated the folding pathway for an RNA pseudoknot by predicting a series of kinetically important intermediates and incorporating information from a nanopore experiment and master equation analysis [34]. Taken together, those results demonstrate that the derived CG force field in the IsRNA model could serve as a potent tool to study RNA folding behaviors and to predict RNA 3D structures.