RNA regulation mechanisms study using physics-informed machine learning

Overview of physics-informed ML paradigm

The physics-informed ML paradigm (see Fig. 1) mainly comprises the following steps. First, this approach employs an neural relational inference (NRI)-informed ML model to examine the potential interactions of structural domains during the dynamic processes of RNA complexes (see Fig. 1a, b) and identify key functional regions combined with RNet²⁴ (see Fig. 1c). The NRI model learns the network dynamics by minimizing the error between the reconstructed and simulated trajectories, then infers edges between nucleotides and residues as latent variables. These learned embeddings capture the key roles of essential elements in conformational transitions, providing insights into RNA regulation mechanisms. Second, the approach uses a regulatory strategy that evaluates the plausibility of potential small-molecular inhibitors, generating mechanistic insights and guiding rational therapeutics design. Specifically, we use ZHMol-RLinter²⁵ for searching and identification of possible candidate small-molecular inhibitors (see Fig. 1d). Together, this integration bridges fundamental understanding of RNA regulation with practical identification of candidate small-molecule inhibitors. To validate the effectiveness of this paradigm, we applied it to two representative RNA complex systems: (i) the P-TEFb/Tat/TAR system, which plays a critical role in the transcriptional activation of the HIV-1 proviral genome; (ii) the aminoacyl-tRNA synthetase (aaRS)/tRNA system, which is essential for genetic translation. Furthermore, we conducted a search to identify potential inhibitors, aiming for precise functional modulation of RNA target activity.

Dynamical regulatory identification in the P-TEFb/Tat/TAR complex

The P-TEFb/Tat/TAR complex is widely regarded as a promising therapeutic target with significant potential for the development of anti-HIV strategies²⁶. We conducted five independent MD simulation trajectories, a total of \(1{{{\rm{\mu }}}}{{{\rm{s}}}}\) simulations, to systematically reveal the dynamic behavior of the P-TEFb/Tat/TAR system. We performed a backbone RMSD analysis to assess the structural stability of the P-TEFb/Tat/TAR complex. The RMSD values were averaged across the five independent trajectories. As shown in Fig. 2b, the system reached equilibrium after approximately 50 ns, as indicated by the stabilization of RMSD values. The secondary structure plot (see Supplementary Fig. 1a) shows that the overall fold is well-preserved over the 200 ns trajectory, indicating stable secondary structure elements throughout the simulation. The Rg-RMSD scatter plot (see Supplementary Fig. 1b) further supports that the system explores a compact and restricted conformational space after equilibration. Together, these analyses demonstrate that the simulations have achieved a stable convergence.

As demonstrated in our previous study²⁷, the key regulatory region of the system is localized at the interaction interface between TAR and the P-TEFb/Tat complex. Specifically, TAR hijacks Tat’s tail to overcome transcriptional pausing (see Fig. 2a). To investigate the regulatory mechanism further, we applied the NRI model using MD trajectories of the Tat/TAR complex. Across all Tat/TAR ensembles sampled of 50 steps, the model accurately reconstructs the trajectories with a mean squared error (MSE) of 0.006 between the truth and reconstruction RMSF (see Supplementary Fig. 2). We also derived the distribution of learned edges between residues (see Fig. 2d) and then constructed a domain interaction map (see Fig. 2e) by grouping adjacent residues/nucleotides into blocks according to their secondary structure (see Fig. 2c). The learned edges frequently occur between the Tat protein and TAR RNA domains, indicating that the Tat/TAR interactions play a crucial role in functional regulation, demonstrating a high connection between TAR RNA and the Tat protein (see Fig. 2f, g). Tat protein exhibits a strong directional preference towards the base-paired regions (R1), bulge loop (R2), and hairpin loop (R4) of TAR RNA, while among the three regions initially considered, only R1 and R2 exhibit this directional preference. This suggests that R1 and R2 act as critical nodes in the Tat/TAR interaction and may serve as potential regulatory targeting sites.

We further employ the network-guided ML method, RNet, to identify the potential functional binding sites of RNA molecules. We computed the binding probabilities of functional sites for TAR RNA to consider RNA’s flexibility based on 5000 frames derived from the P-TEFb/Tat/TAR MD simulation trajectories. As shown in Fig. 3, the curve illustrates the average predicted binding probabilities throughout the entire MD trajectory frames. The TAR RNA sequence displays two peaks (A7-G10 and C23-C25) that correspond to regions with high binding potential (shown in the blue shaded region of Fig. 3a). These regions precisely align with the bulge loop (R2), indicating that the bulge loop domain (nucleotides highlighted in maroon in Fig. 3b) shows the highest functional binding site probability, emphasizing its potential as a primary target for small-molecule inhibitor regulation. Additionally, we ranked the functional binding site probabilities for each nucleotide, with the top five being A7, G10, U24, U8, and A11. Among these, A7, G10, U24, and U8 are found within the loop region, while A11 is located near the loop region. These nucleotides are likely to play a crucial role in the structural dynamics and functional regulation of TAR RNA. Targeting these nucleotides could interfere with the interactions between Tat and TAR, preventing HIV-1 infection from transcriptional elongation. In fact, the multiple sequence alignment of representative HIV-1 subtypes highlights strong conservation in the U8 of the bulge loop (see Supplementary Fig. 3a). Structurally, a comparison between HIV-1 (PDB code: 6MCE) and HIV-2 (PDB code: 1AKX) TAR reveals bulge loop regions in both the secondary (see Supplementary Fig. 3b, c) and tertiary structures (see Supplementary Fig. 3d, e). Functionally, previous research indicates that nucleotides in the bulge are conserved and essential for Tat interaction. Specifically, U8 is fully conserved across all HIV isolates and is the only base that cannot be replaced by any of the other three natural bases (A, C, and G). Specifically, U8 forms hydrogen bonds with residues in TAR, creating a critical tertiary structure for Tat binding²⁸. Importantly, studies have shown that methylation at U8’s N3 position disrupts high-affinity TAR binding²⁹.

Computational validation through removal of the bulge loop

Given the potential of the bulge loop region, we conducted a controlled experiment by removing this bulge loop structural element from the TAR RNA and performing MD simulations (see the “Methods” section for additional details). As illustrated in Fig. 4a, we present the P-TEFb/Tat/TAR system with the bulge loop region removed, referred to as P-TEFb/Tat/TAR-Delta. To quantitatively assess the structural stability of the P-TEFb/Tat/TAR-Delta complex, we also carried out an RMSD analysis. The RMSD values were calculated and then averaged across five independent trajectories to ensure statistical reliability. As shown in Fig. 4c, the system achieved convergence at approximately 100 ns. However, the RMSD values exhibited significantly greater fluctuations compared to those of the P-TEFb/Tat/TAR system with the bulge loop region present. To further investigate the molecular basis of these fluctuations, we conducted RMSF analysis on the backbone atoms of the P-TEFb/Tat/TAR complex. Removing the bulge loop in TAR RNA led to significant fluctuations in the nucleotides, particularly near the bulge loop region, while other areas remained relatively stable (see Supplementary Fig. 4a, b). These findings suggest that the absence of the TAR bulge loop increases structural flexibility near the Tat binding sites, potentially disrupting the Tat/TAR interface, which is essential for the complex’s functional regulation.

To further investigate the mechanism after removing the bulge region, we applied the NRI model using MD trajectories of the Tat/TAR-Delta complex within the P-TEFb/Tat/TAR-Delta system. From all Tat/TAR ensembles sampled over 50 steps, the model achieves highly accurate trajectory reconstruction with a mean squared error (MSE) of 0.004 between the truth and reconstruction RMSF (see Supplementary Fig. 5). We derived the distribution of learned edges between residues and nucleotides (see Fig. 4d) and constructed a domain interaction map (see Fig. 4e) by grouping adjacent residues into blocks based on the secondary structure (see Fig. 4b). Compared to the P-TEFb/Tat/TAR system, the influence of learned edges between Tat and TAR becomes weaker (see Fig. 4g, h). To quantitatively characterize these changes, we measured the differences in interaction strengths between P-TEFb/Tat/TAR-Delta and P-TEFb/Tat/TAR (see Fig. 4f). The heatmap shows a reduction in nearly all interaction strengths. This finding further supports the crucial role of the bulge region as a mediator of the Tat/TAR interface in mechanistic regulation.

Regulatory application through ML-based inhibitor identification

Based on our findings, we can evaluate potential small-molecule inhibitors that target the RNA to disrupt the interaction between Tat and TAR (see Fig. 5a). We applied our ZHMol-RLinter methods, which can identify the binding probabilities between loop motifs and inhibitors. This approach demonstrated that loop motifs are highly likely to interact with small-molecule inhibitors, making them primary targets for inhibitor binding analysis. Additionally, the secondary structure of TAR RNA has two characteristic loop regions: the bulge loop and the hairpin loop (see Fig. 5b).

Previously, our research showed five potential inhibitors (110FA, 115FA, F07#13, AM6538, and DB00594) that may target P-TEFb/Tat/TAR³⁰. We applied ZHMol-RLinter to predict the binding preferences between the loop motifs and inhibitors. The analysis with ZHMol-RLinter revealed distinct binding preferences among the selected inhibitors. All five inhibitors demonstrated higher binding probabilities for the bulge loop region compared to the hairpin loop region (see Fig. 5c, d). 110FA exhibited the strongest binding probability toward the bulge loop (see Fig. 5d), suggesting it may be the most promising candidate for targeted TAR RNA inhibition.

We performed molecular docking between 110FA and TAR RNA (labeled as P-TEFb/Tat/TAR-L1), followed by MD simulations (see “Methods” for details). Due to the competitive binding of the small-molecule inhibitor at the RNA binding site, the ARM region of Tat remains free (see Fig. 5a). To evaluate the stability of the simulations upon inhibitor binding, we calculated the backbone RMSD across five 200 ns independent trajectories, amounting to a cumulative simulation time of \(1{{{\rm{\mu }}}}{{{\rm{s}}}}\). As shown in Fig. 5e, both systems stabilized after 50 ns of simulation, indicating the convergence of the trajectories. While the RMSD values of the P-TEFb/Tat/TAR-L1 complex were consistently much higher than those of the P-TEFb/Tat/TAR complex during the simulations. This increased RMSD suggests that the binding of inhibitor 110FA enhances the structural fluctuations within the P-TEFb/Tat/TAR complex. To further investigate structural fluctuations, we conducted the backbone RMSF analysis on the P-TEFb/Tat/TAR-L1 complex. As shown in Supplementary Fig. 6a, the system displayed overall structural stability, except for TAR RNA. Significant fluctuations were observed in TAR RNA, especially near the bulge loop region (Supplementary Fig. 6b). These findings further confirm that binding of the inhibitor 110FA induces considerable structural destabilization between TAR RNA and the Tat protein, potentially disrupting the Tat/TAR interface.

Dynamical regulatory identification in the aaRS/tRNA complex

We further extend our analysis to the second system, namely the aaRS/tRNA complex, which is essential for the aminoacylation of tRNA, a critical step in protein synthesis³¹. In this process, aaRS charge tRNAs with their cognate amino acids, ensuring accurate translation of the genetic code³². Specifically, the amino acid is transferred to the 3’ end of the cognate tRNA, and the synthetase distinguishes among a large pool of cellular tRNAs by recognizing particular nucleotides called identity elements. This precise recognition ensures each tRNA is charged with the correct amino acid, preserving the accuracy and fidelity of protein synthesis. As a result, inhibiting tRNA aminoacylation has been confirmed as an effective antimicrobial strategy³³.

We obtained the aaRS/tRNA system (see Supplementary Fig. 7a, right) trajectory from previous research³⁴. To explore the underlying regulatory mechanisms, we applied the NRI model to the MD trajectory of tRNA. Across all sampled tRNA ensembles of 50 steps, the model accurately reconstructed the trajectories, achieving an MSE of 0.0004 between the ground truth and reconstructed RMSF (see Supplementary Fig. 8). We also derived the distribution of learned edges between nucleotides (see Supplementary Fig. 7b) and built a domain-level interaction map (see Supplementary Fig. 7c) by grouping adjacent nucleotides into blocks based on the secondary structure (see Supplementary Fig. 7a, left). The learned edges among the D-loop, variable loop, and anticodon arm suggest that interactions between these regions are essential in the aminoacylation process of tRNA (see Supplementary Fig. 7d, e). Additionally, we calculated the shortest pathways from nucleotides in the anticodon loop to those in the acceptor arm using the learned edges, representing key allosteric communication routes within the tRNA molecule (see Supplementary Fig. 7f). The pathway mainly passes through the variable loop, indicating their relative importance in improving global connectivity and facilitating interactions that enhance allosteric signaling.

We further use the network-guided ML method, RNet, to predict potential functional binding sites for tRNA molecules. We calculated the binding probabilities of functional sites for tRNA, accounting for RNA’s flexibility based on 3200 frames from the aaRS/tRNA MD simulation trajectories. As shown in Supplementary Fig. 9a, the curve represents the average predicted binding probabilities across the entire MD trajectory frames. The tRNA sequence shows four distinct peaks (A7-G10, U19-G23, A46-G48, and G56-U59), which correspond to regions with high binding potential (indicated by the blue shaded areas in Supplementary Fig. 9a). These regions align precisely with the D-loop, T-loop, and variable loop domains (nucleotides highlighted in maroon in Supplementary Fig. 9b), emphasizing their potential as primary targets for small-molecule inhibitors. Indeed, the critical binding regions identified by RNet match the key interaction sites that sustain the L-shaped stability of tRNA, a structural requirement for its functional activity³⁵.

Experimental evidence further supports the findings that small-molecule inhibitors specifically bind to the variable loop and D-loop regions of tRNA, disrupting the structural integrity needed for efficient aminoacylation. And the aminoglycoside antibiotic Neomycin B has been reported to inhibit the in vitro aminoacylation of E. coli tRNA^Phe 36.

Performance evaluation of the physics-informed ML models

To further validate the physics-informed ML methods, we conducted a comparative analysis of these tools against traditional approaches. One key benefit of using the NRI ML method in our framework is that it can learn hidden interaction edges directly from the MD simulation data, rather than relying on past correlation-based network analysis measures. First, the NRI model doesn’t just construct static interaction networks. Instead, it learns how dynamic interaction edges work in a way that lets it generate accurate MD trajectories. This type of validation demonstrates that the learned edges contain sufficient dynamical information, a benefit that conventional network-based approaches can’t match. The Pearson correlation coefficient between the learned node weights and residue-level RMSF values is 0.59 (see Supplementary Fig. 10a), indicating a strong connection between the learned interaction strengths and changes in structural stability. On the other hand, traditional network-based methods, such as correlation or contact frequency analysis, including the dynamical cross-correlation matrix (DCCM, the details of DCCM can be found in Supplementary Note 1), can’t directly construct trajectories. The DCCM only has a Pearson correlation of −0.34 (see Supplementary Fig. 10b), which isn’t markedly compared to NRI. Further, we compared the NRI model with DCCM using network shortest-path analysis. Specifically, we constructed networks where nodes represent nucleotides or residues, and edges are defined by spatial distances less than 20 Å, with edge weights assigned from either NRI or DCCM. Then, we applied shortest-path analysis, a standard method in network studies, to calculate how often paths go through RNA nucleotides. The comparison showed that NRI (see Supplementary Fig. 10c) has sharper and more distinct peak patterns than DCCM (see Supplementary Fig. 10d), with most peaks located near regulatory bulge regions. This highlights the advantage of NRI. Overall, the comparison results show that (i) NRI gathers more detailed dynamic information and (ii) it provides more accurate and notable identification of key nucleotides.

We also compared the performance of the binding site prediction method RNet with other methods, including RNAsite, RBind, and Rsite^37,38,39. RNAsite is also a ML-based method. RBind is a physics-based method grounded in complex network theory that analyzes binding sites by leveraging network properties like degree and closeness centrality. In contrast, Rsite is an approach based on Euclidean distance that identifies binding sites by calculating the spatial distances between nucleotides and other molecular components. As shown in Supplementary Fig. 11b, Rsite failed to accurately predict the bulge loop region, with its predicted binding sites scattered across various positions on the RNA. The precision of its predictions for the bulge loop region is 42.9%. This indicates that relying solely on simple Euclidean distance metrics is insufficient for accurately identifying specific functional areas. On the other hand, due to RBind’s strict cutoff criteria, it was unable to detect the binding sites effectively with 0% precision for the bulge loop region (see Supplementary Fig. 11a). The precision of RNAsite in predicting binding sites on the bulge loop is only 0.333 (see Supplementary Fig. 11c). In contrast, RNet, which utilizes ML to capture complex network features, successfully identified high-probability binding sites near the bulge loop region, achieving a precision of 80% in the Top 5 predictions, as shown in Supplementary Fig. 11d, representing a 37.1% improvement compared to Rsite and 46.7% to RNAsite.

Additionally, we evaluated the ability of ZHMol-RLinter in the pipeline to identify RNA-small molecule binding preferences using a set of experimentally determined RNA-small molecule PDB structures (see Supplementary Table 1). We started with a non-redundant set of 31 structures from the published benchmark RL98^25,40. Subsequently, we added five recently released RNA-small molecule complexes from the PDB (post after April 2024). These were combined into a curated test set named RSM36, which was used to assess the pipeline’s performance in identifying corresponding small-molecule binding sites. On RSM36, the pipeline achieved an accuracy of 0.63, with a precision of 0.62 and a recall of 0.62 in predicting small molecules that bind to correct RNA motifs. When using a more relaxed criterion for determining whether the small molecule binds to the correct RNA chain, the accuracy increased to 0.74, the precision to 0.77, and the recall rose dramatically to 0.94 (see Supplementary Fig. 12). These results show that the pipeline can reliably identify small molecules binding at relevant sites, providing valuable insights for discovering RNA-targeting inhibitors.