Unraveling the Mechanism of Drug Binding to SARS-CoV-2 RNA
arXiv:2604.14906v1 [physics.bio-ph] 16 Apr 2026
Pseudoknot with Thermodynamics-Driven Machine Learning Mariia Ivonina∗,† and Jakub Rydzewski∗,‡ †RIKEN Center for Interdisciplinary Theoretical and Mathematical Sciences (iTHEMS), RIKEN, Wako 351-0198, Japan ‡Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziądzka 5, 87-100 Toruń, Poland E-mail: [email protected]; [email protected]
1
Abstract The SARS-CoV-2 RNA pseudoknot is a promising target for antiviral intervention, as it regulates the efficiency of −1 programmed ribosomal frameshifting (−1 PRF), a mechanism that is essential for viral protein synthesis. The pseudoknot represents a viral RNA sequence composed of helical stems that adopts two long-lived topologies, threaded and unthreaded. Ligand-induced distortion of this fold is thought to underlie the susceptibility of −1 PRF to small-molecule inhibitors. Resolving these distortions from unbiased molecular dynamics (MD) requires collective variables (CVs) that isolate the slowest dynamic modes of the RNA–ligand system from the high-frequency fluctuations. Here, we use spectral map (SM), a thermodynamics-driven machinelearning method, to learn such CVs directly from all-atom MD trajectories of the SARS-CoV-2 RNA pseudoknot in complex with the −1 PRF inhibitor merafloxacin and two related analogs. We examine both threaded and unthreaded pseudoknot topologies and consider the neutral and ionized ligand forms relevant at physiological pH. Our learned free-energy landscapes show that ligand-induced destabilization is topology-selective: merafloxacin and its analogs destabilize the S2 stem in the threaded pseudoknot, whereas in the unthreaded pseudoknot, destabilization shifts to the S1 and S3 stems. Spectral analysis of the learned CVs further clarifies ligand effects: for example, the zwitterionic form of merafloxacin uniquely imposes slow dynamics on the otherwise spectrally featureless unthreaded pseudoknot. Additionally, the degree to which each ligand reshapes the free-energy landscape corresponds with its experimentally reported antiviral potency. This modulation is mainly governed by kinetics rather than thermodynamics, and the neutral and zwitterionic forms of merafloxacin differ qualitatively in their mechanisms within the same RNA topology. Overall, these results clarify how pseudoknot topology, ligand type, and protonation state shape the slow conformational dynamics of viral RNA and establish physiological protonation as an essential factor for modeling RNA-targeted drug action.
2
1
Introduction
Programmed −1 ribosomal frameshifting (−1 PRF) is an essential regulatory mechanism used by many RNA viruses, including coronaviruses and HIV, to control the synthesis of proteins required for viral replication. 1 In this process, the translating ribosome slips backward by one nucleotide on the viral RNA and enters an alternative reading frame, altering the encoded viral protein sequence. An efficient −1 PRF is required to maintain the correct balance of viral proteins produced, making the RNA elements responsible for frameshifting particularly attractive antiviral targets. The frameshifting event is typically stimulated by two cis-acting RNA elements: a heptanucleotide slippery sequence of the form X XXY YYZ and a downstream RNA pseudoknot. The pseudoknot acts as a mechanical roadblock that pauses ribosomal elongation and increases the probability of ribosomal shift into the −1 frame at the slippery sequence. Because the efficiency of this process depends on specific RNA conformations, the pseudoknot emerges as an appealing target for antiviral drug development. Unlike protein-targeting strategies that intervene after viral gene expression has begun, targeting viral RNA can suppress propagation at an earlier stage, before viral proteins are fully produced and assembled. 2–5 Following the outbreak of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in late 2019, substantial effort has been devoted to developing vaccines and antiviral therapies to reduce viral transmission and mortality. In SARS-CoV-2, −1 PRF is estimated to occur in more than 50% of all translation events, 6 making the frameshift-stimulatory pseudoknot a particularly promising therapeutic target and highlighting the need to understand the structural basis of its function. Experimental techniques such as X-ray crystallography, 7,8 cryo-EM 9–11 and NMR, 12–14 as well as molecular dynamics (MD) simulations, 15–18 captured the conformational heterogeneity of the pseudoknot in solution, showing that the SARSCoV-2 frameshift signal adopts a three-stem H-type pseudoknot that can populate multiple conformers, including 5′ -end threaded and unthreaded topologies. Overall, this indicates that, while the pseudoknot retains relatively stable secondary structure, it can undergo 3
substantial tertiary rearrangements. Several inhibitors that affect −1 PRF in coronaviruses via interaction with the pseudoknot have been identified, demonstrating the feasibility of RNA-targeted antiviral modulation. 19 MTDB, a small-molecule inhibitor, retains activity against virus replication due to conservation of the pseudoknot structure across SARS-CoV-1 and SARS-CoV-2. 20–22 Merafloxacin, a fluoroquinolone antibacterial, remains an effective −1 PRF inhibitor against naturally occurring SARS-CoV-2 mutants. 23 Nafamostat, a serine protease inhibitor, demonstrates broad-spectrum activity against multiple coronavirus-derived frameshift signals, with effects comparable to those of merafloxacin. 22,24 Geneticin, an aminoglycoside antibiotic, exhibits antiviral activity against multiple SARS-CoV-2 variants. 25 Aminoquinazoline derivatives also act as inhibitors of −1 PRF through interaction with the SARS-CoV-2 pseudoknot structure. Additional compounds with similar pseudoknot-targeting mechanisms have also been identified. 26–28 Despite this progress, the dynamical basis of ligand recognition by the SARS-CoV-2 pseudoknot remains incompletely understood. Prior computational results have shown that ligand binding can strongly depend on RNA topology. For example, 5′ -end threading can create clefts that are absent in simpler folds, 20 and some ligands remain bound to the pseudoknot on the nanosecond timescale. 28 Nevertheless, the mechanism remains unresolved: it is still unclear whether a ligand stabilizes a unique binding mode, migrates among multiple transient interaction sites, or actively reshapes the RNA conformational ensemble. Given that RNA binding pockets are often shallow, dynamic, and topology-dependent, a detailed description of ligand recognition therefore requires analysis of how the RNA free-energy landscape (FEL) changes upon ligand association. 29–31 A further consideration specific to RNA-targeted ligands is the charge state of the bound compound. Because the RNA backbone is polyanionic, electrostatic contacts with positively charged or zwitterionic moieties often dominate recognition and can alter affinity, selectivity, and association kinetics by orders of magnitude. 5 For example, this is evident in merafloxacin:
4
like other fluoroquinolones, it carries both an acidic carboxylate and a basic amine, and is predominantly zwitterionic at physiological pH. 32 Thus, its inhibiting functions might be associated with strong electrostatic effects. Conformational changes relevant to RNA function span a broad range of timescales, from picoseconds for puckering and bond rotations, to nanoseconds for base-pair opening, and to much longer timescales for large-scale conformational switching. 33–36 Experimental methods provide only limited access to this microscopic dynamics, yet such information is essential for a mechanistic understanding of RNA-mediated processes. MD simulations offer an atomistic view of RNA motions. However, quantitatively characterizing the thermodynamics and kinetics of RNA conformational transitions requires collective variables (CVs) that can identify and disentangle slow dynamic modes. 37–39 For RNA folding, CVs must distinguish between long-lived stem-connected and stem-disrupted states, locate transition regions, and capture motions associated with barriers of several kB T . To address this challenge, we employ a thermodynamics-driven machine-learning (ML) framework 40,41 to construct CVs for SARS-CoV-2 pseudoknot folding. Specifically, we use spectral map (SM), 42–44 a recently developed method for performing dimensionality reduction using neural networks (NNs) and identifying slow dynamical modes without any prior knowledge of the underlying reaction coordinate. SM constructs CVs by maximizing timescale separation, thus directly targeting slow dynamics. This characteristic distinguishes SM from conventional descriptors and makes it particularly effective for characterizing ligand-induced conformational transitions in RNA pseudoknots, where the coupling between ligand association and large-scale RNA rearrangement is difficult to quantify. In this work, we investigate how the SARS-CoV-2 frameshift-stimulatory pseudoknot interacts with the −1 PRF inhibitor merafloxacin and two related analogs proposed by Sun et al. 23 We analyze how ligand type, protonation state, and pseudoknot topology jointly determine RNA structural stability and drive conformational transitions. To achieve this, we perform extensive all-atom MD simulations of both threaded and unthreaded pseudoknot
5
models in ligand-free and ligand-bound states, and analyze the resulting trajectories with SM to identify the slowest modes and trace the minimum free-energy pathways in the FELs of each RNA–ligand system. We find that ligand-induced destabilization is selective in the region of the pseudoknot it targets: the threaded topology predominantly loses native contacts at the S2 stem, whereas the unthreaded topology disrupts native contacts in the S1 and S3 stems. The NN-learned slow modes reveal a qualitative dynamical distinction between the two topologies, and binding of the merafloxacin zwitterion imposes a clear structure on the otherwise featureless unthreaded pseudoknot, indicating that the inhibitor actively slows the dynamics rather than stabilizing a preexisting RNA conformer. The degree to which each ligand reshapes the learned FEL follows the experimentally reported antiviral ordering, 23 and this modulation is governed primarily by kinetic rather than thermodynamic control: free-energy differences between metastable states are small, while transition barriers span more than an order of magnitude across the systems studied. Finally, the neutral and zwitterionic forms of merafloxacin yield qualitatively different FELs on the same RNA topology, confirming that physiological protonation is a load-bearing variable in the mechanistic picture of RNA–ligand recognition.
2
Thermodynamics-Driven Machine Learning
SM is a thermodynamics-driven ML algorithm, grounded in a specific definition of what makes a CV slow. In unbiased MD, trajectories sample the equilibrium distribution p(x) ∝ exp(−βU (x)) of the underlying potential U (x). As a result, pairwise proximity between configurations encodes thermodynamic similarity. Based on this connection, the densitynormalized Gaussian kernel used below converges, for a large sample and small bandwidth, to the generator of overdamped Langevin dynamics on U . Equivalently, it converges to the Fokker–Planck operator associated with the equilibrium ensemble. 45,46 Its eigenvalues λk
6
correspond directly to physical relaxation timescales, τk = −∆t/ ln λk . Thus, the dominant eigenmodes approximate the system’s slowest relaxation processes. These describe transitions that cross the highest free-energy barriers, defining the separation of metastable basins in the FEL.
2.1
Dimensionality Reduction with an Encoder
SM employs an NN encoder fθ (x) that learns slow CVs directly from MD data. 42–44 SM can capture nonlinear relationships in the data while optimizing for dynamical slowness through a spectral loss function. The encoder fθ : Rn → Rd (see Fig. 1) maps high-dimensional descriptors x = (x1 , . . . , xn ) to a reduced space of d CVs, z = fθ (x) = (z1 , . . . , zd ), where d ≪ n and θ denotes the NN weights. Spectral Map
FEL F (z)
x1 w11
Descriptors x
x2
CVs z h1 z1
wm1
x3
h2 w1n
x4
.. .
···
.. .
.. . zd
hm wmn
Lλ + L D + L L
xn Encoder z = fw (x)
Figure 1: Architecture of spectral map (SM). The encoder fθ maps high-dimensional descriptors x to a low-dimensional latent space of CVs z. The total loss combines three terms: Lλ (spectral loss to capture slow dynamics), LD (decorrelation loss to ensure CVs represent distinct modes), and LL (graph Laplacian regularizer for smoothness). The learned CVs z are used to estimate the corresponding free-energy landscape (FEL) F (z). Here, the CVs can be used to describe transitions between metastable states (A, B, and C): z1 between A/C and B, and z2 between A/B and C. These transitions are tied to the slowest dynamic modes in the system.
7
2.2
Markov Transition Matrix
To learn slow CVs, we train the encoder by maximizing the separation between the slowest (dominant in magnitude) and fastest eigenvalues of a Markov transition matrix M (zk , zl ) constructed in z space. The algorithm to construct M proceeds as follows: 42–44
1. We compute a Gaussian kernel G(zk , zl ) = exp − 1ε ∥zk − zl ∥2 , where ∥zk −zl ∥2 denotes the squared Euclidean distance between CV samples and ε is the bandwidth parameter. 2. To correct for non-uniform sampling density, we employ an anisotropic diffusion kernel: 46 G(zk , zl ) K(zk , zl ) = q , ϱ(zk )ϱ(zl ) where ϱ(zk ) =
(1)
l G(zk , zl ) is a kernel density estimate.
P
3. We compute the Markov transition matrix M by row-normalizing K:
M (zk , zl ) = D(zk )−1 K(zk , zl ),
where D(zk ) =
(2)
n K(zk , zn ). The matrix M models a discrete Markov chain in the CV
P
space, which approximates the Fokker–Planck operator of the equilibrium ensemble. 4. We compute the symmetric conjugate Ms = D1/2 M D−1/2 . This similarity transformation preserves the eigenvalues of M while yielding a symmetric matrix amenable to symmetric eigendecomposition. 5. We perform an eigendecomposition Ms Φ = λΦ, obtaining eigenvalues λ0 = 1 > λ1 ≥ · · · ≥ λN −1 , which are used to define the spectral loss.
2.3
Spectral Loss
The transition matrix is then used to define a loss function that SM optimizes during training to obtain slow CVs. In this work, we use the negative of cumulative spectral weight of the 8
k leading eigenvalues λ0 , . . . , λk−1 as the core loss: Pk−1
λi , i=0 λi
Lλ = − PNi=0 −1
(3)
where N is the batch size. As the dominant eigenvalues of M correspond to the longest timescales in the system dynamics, this loss term ensures that the slowest degrees of freedom are captured by the encoder. Unlike the spectral gap ∆λ = λk−1 − λk used in our previous work, 42–44 the cumulative spectral weight is less sensitive to the choice of k. An incorrect guess for the number of metastable states can significantly affect CVs learned via the spectral gap, whereas the cumulative weight provides a smoother loss. In addition to Lλ , we employ two auxiliary losses. First, we use a decorrelation loss:
LD =
X 1 R2 , d(d − 1) k̸=l kl
(4)
where Rkl is the correlation coefficient between CVs zk and zl , and the factor d(d − 1) counts the off-diagonal pairs. The squared correlation ensures a positive definite penalty. This loss encourages decorrelation among the learned CVs, ensuring that each CV represents a distinct slow mode. Second, a graph Laplacian regularizer is added:
LL =
1X Ms (zk , zl )∥zk − zl ∥2 . 2 kl
(5)
This term penalizes large differences in CV values between states with high transition probability, enforcing smoothness of the learned representation with respect to the underlying dynamics. The total loss combines these three components: L = Lλ + LD + LL .
2.4
Orthogonality
To stabilize training and prevent scale collapse, we enforce orthogonality on the weight matrices in the hidden layers. This constraint prevents gradient explosion during backprop9
agation and ensures that the NN does not collapse the input space. With the orthogonality constraints, the encoder is forced to preserve the scale of its representation, which directly affects the spectral properties of the learned CVs. The training is focused on arranging CV space so that a few leading eigenvalues are genuinely large (reflecting slow dynamics), rather than trivially collapsing the scale to reach the minimum of the spectral loss. For a weight matrix W ∈ Rm×n connecting neighboring linear layers with n inputs and m outputs, we impose W ⊤ W = In if m ≥ n (tall matrix) and W W ⊤ = Im if m < n (wide matrix), where I denotes an identity matrix. The orthogonality is enforced using the Cayley parametrization via Householder reflectors. 47
3
Computational Details
3.1
RNA Models and Ligands
The initial structures of the SARS-CoV-2 pseudoknot used for MD simulations were retrieved from the Nucleic Acid Knowledgebase (NAKB) (https://nakb.org/). We used two cryo-EM pseudoknot structures, 6XRZ 9 and 8VCI. 11 These two models differ in their fold topology: in 6XRZ, the 5′ RNA strand is threaded through the junction between stem 1 (S1) and stem 3 (S3), whereas in 8VCI the 5′ end remains outside this junction (Fig. 2). Hereafter, we refer to these RNA conformers as the threaded (RT) and unthreaded (RU) forms, respectively. To make both conformers comparable in this study, we extracted 67 nucleotides from 6XRZ and 8VCI, corresponding to the main elements of the pseudoknot: stems 1, 2, and 3 (S1, S2, S3) and loops 1 and 2 (L1, L2). To preserve the native folding of the pseudoknot at its boundaries, we included seven nucleotides upstream at the 5′ end (GGGUUU) and four nucleotides downstream at the 3′ end (UUUG). In total, 78 nucleotides constitute the SARS-CoV-2 pseudoknot body in both threaded and unthreaded conformations simulated in this study, as shown in Fig. 2. 10
(a)
(b)
(c)
U U U G 3’ 20 A U 13546 G C C G C G S2 C G 70 C G G U U G C A U A A U C S1 U A 30 A L2 G C U U A A 10 G C G C C G G U 1 5’ G G G U U U G U C G 13469 G C 40 G U S3 C G 60 A U A C G A U G C G A C U L1 A G C U A
(d)
U U U G 3’ A U G C C G C G S2 C G 70 C G G U U G C A U A A U C L2 S1 U A 30 A G C U A 10 G C G C 1 5’ G G G U U U G C G U U A G U C G G C 40 G U S3 C G 60 A U A C G A U G C G A C U L1 A G C U A 20
50
50
Figure 2: Secondary structures of the SARS-CoV-2 −1 ribosomal frameshift-stimulatory pseudoknot derived from crystallographic data. (a) Threaded (RT) conformation, 9 in which the 5′ -end passes through the junction formed by stems 1 (S1), 2 (S2), and 3 (S3), thereby producing a topological tertiary fold. (b) 3D representation of the RT model. (c) Unthreaded (RU) conformation, 11 in which this topological fold is absent. (d) 3D representation of the RU model. Nucleotides are colored according to their secondary-structure elements (stems: S1 in green; S2 in orange; S3 in purple; and loops L1, L2 in dark gray). The 5′ and 3′ termini are depicted in light gray.
11
To investigate ligand recognition by the SARS-CoV-2 pseudoknot, we focused on antibiotics from the fluoroquinolone family. These compounds share a characteristic core structure comprising a bicyclic ring system bearing a carboxylic acid group at the 3-position and a carbonyl group at the 4-position. We selected merafloxacin (MERA) together with two structurally related analogs, denoted A1 and A2, which differ in their C7 substituents, following the nomenclature of the original experimental study by Sun et al. 23 Merafloxacin was previously identified as an inhibitor of −1 PRF in coronaviruses and was shown to restrict viral replication in cell-based assays, making it a relevant reference compound for RNA-targeted antiviral design. 22,23,48 Compounds A1 and A2 were included as closely related analogs with weak −1 PRF inhibition to compare how changes in ligand structure influence RNA binding and their effects on the RT and RU structures. To account for physiological conditions, we considered the relevant protonation states of each ligand at pH≈7.0. Based on the reported dissociation constants of fluoroquinolones, 32 MERA was modeled in both its neutral (MERA0 ) and zwitterionic (MERA± ) forms. In MERA± , the quinolone carboxyl group is deprotonated, while the side-chain amine is protonated. By contrast, A1 and A2 were modeled in their neutral (A10 , A20 ) and anionic (A1− , A2− ) forms, consistent with the presence of a quinolone carboxyl group and the absence of a side-chain amine group (Tab. 1). The three selected ligands in their neutral and ionized forms are shown in Fig. 3. The RNA–ligand complexes were prepared using docking with AutoDock Vina 49,50 and vinardo 51 scoring function. For each ligand, docking calculations were carried out against both the RT and RU models of the pseudoknot, and the three top-ranked poses were retained as starting structures for MD simulations (see Fig. S8).
3.2
Molecular Dynamics
We performed MD simulations for RT and RU without any ligands, as well as for each RNA– ligand complex with the ligands present in the neutral and ionized states (Tab. 1). The 12
HN HO
F N
N OH
F O
F
F N
H 2N
N
O
OH
F
O
O
MERA0
N
N OH
F O
O
A10
O
A20
H 2N HO
F N
O
F O
F
F
N
N
N O
F O
O
MERA±
H 2N
N O
N O
F
O
A1-
O
O
A2-
Figure 3: Ligands considered in this study: (left) Merafloxacin in the neutral (MERA0 ) and zwitterionic (MERA± ) forms; (center) Analog 1 molecule in its neutral (A10 ) and anionic (A1− ) forms; (right) Analog 2 molecule in its neutral (A20 ) and anionic (A2− ) forms.
Table 1: The complete list of the systems considered in this work. Abbreviation RT RT-MERA0 RT-MERA± RT-A10 RT-A1− RT-A20 RT-A2− RU RU-MERA0 RU-MERA± RU-A10 RU-A1− RU-A20 RT-A2−
RNA conformation Threaded Threaded Threaded Threaded Threaded Threaded Threaded Unthreaded Unthreaded Unthreaded Unthreaded Unthreaded Unthreaded Unthreaded
13
Ligand — Merafloxacin Merafloxacin (zwitterion) Analog 1 Analog 1 (anion) Analog 2 Analog 2 (anion) — Merafloxacin Merafloxacin (zwitterion) Analog 1 Analog 1 (anion) Analog 2 Analog 2 (anion)
simulations were carried out in GROMACS 2024.2 using the CHARMM36 force field 52 (July 2022 release) together with the CHARMM-modified TIP3P water model. Ligand topologies and parameters were assigned using the CHARMM General Force Field (CGenFF). 53 After system preparation, each system underwent 1 ns of NVT equilibration followed by 1 ns of NPT equilibration. Temperature was controlled at 310 K with the velocity-rescaling thermostat, 54 and pressure was maintained at 1 atm with stochastic cell rescaling. 55 Unrestrained production runs were then performed for 0.4–3 µs with a 2 fs time step. We used a dodecahedron box with a 15 Å margin. The resulting trajectories were analyzed using the MDTraj library. 56 To assess the structural stability and conformational dynamics of the pseudoknot in the MD simulations, we computed the root-mean-square deviation (RMSD) of the full structure, per-residue root-mean-square fluctuation (RMSF), and the fraction of native contacts, Q, for each of the three stems over time. These quantities were calculated from mass-centered trajectories after structural alignment. The fraction of native contacts, Q, was evaluated following the definition of Best, Hummer, and Eaton. 57 Q was calculated only for the stem regions of the pseudoknot in order to specifically monitor the preservation of the native helical core. For a given conformation x, Q(x) was computed as a smooth contact-based measure relative to the reference native structure, using all pairs of heavy atoms (k, l) belonging to stem residues separated by more than three positions in sequence and whose reference distance satisfied x0kl < 0.45 nm. The contact contribution of each pair was weighted according to:
Q(x) =
n −1 1X 0 1 + eα(xkl −γxkl ) , c k>l
(6)
where c is the number of native contacts, xkl is the heavy-atom pairwise distance between atoms k and l, x0kl is the distance in the reference structure, and the parameters were set to α = 50 nm−1 and γ = 1.5, as recommended for all-atom simulations.
14
3.3
Machine Learning
We define the descriptors by calculating pairwise distances between 1-position nitrogen atoms in nucleobases. For each system, we take 100 descriptors with highest variance. Each dataset consists of samples collected every 50 ps (every 5th frame from trajectories saved at 10 ps intervals). We employ NNs with two hidden layers of 75 nodes each and ELU activation functions. Each network contains approximately 13,000 trainable parameters. The output comprises two CVs. Training uses batches of 1000 samples, with 20% of the data reserved for validation. Training proceeds until early stopping is triggered when the validation loss reaches 10−4 precision with a patience of 20 epochs. We use the Adam optimizer with a learning rate of 10−3 and default parameters. In all cases, we enforce orthogonality on NN parameters. For the spectral loss (Eq. 3), we use k = 3 dominant eigenvalues calculated from the symmetric conjugate of the Markov transition matrix (Eq. 2). Following standard practice in diffusion map methods, we set the Gaussian kernel bandwidth ε to the median of pairwise Euclidean distances. To assess the quality of the training protocol, we analyzed the stability of the loss function across epochs. All systems converge within ∼100 epochs with agreement between the training and validation losses, ruling out overfitting (Fig. S1). The spectral loss Lλ (negative cumulative spectral weight) dominates the total loss, while the decorrelation (LD ) and Laplacian (LL ) losses remain near zero throughout, indicating the two learned CVs are easily kept decorrelated and smooth (Fig. S2 and Fig. S3). Additionally, the Gaussian kernel bandwidth ε, estimated based on the median of pairwise distances between samples in CV space (Eq. 1), converges to a narrow range (0.15–0.25) for all systems (Fig. S4), confirming that orthogonal weight constraints in the encoder prevent scale collapse, and thus degeneration of the spectral characteristics of the learned CVs. For ML, we use our development version of the spectre package, implemented using the PyTorch library 58 and the PyTorch Lightning framework (https://github.com/ 15
Lightning-AI/lightning/).
A development version of spectre is available on Code-
berg (https://codeberg.org/jkrd/spectre/) with a documentation at https://jkrd. codeberg.page/spectre/.
3.4
Free-Energy Landscapes
We compute free-energy landscapes as F (z) = −β −1 log P (z), where β = 1/kB T is the inverse temperature with the Boltzmann constant denoted as kB and temperature as T . We calculate the marginal probability density in CV space by kernel density estimation (KDE). In KDE, we use Gaussian kernels with bandwidth estimated from Scott’s rule. We estimate freeenergy differences as ∆FA→B = −β −1 log
R
−βF (z) / B dz e−βF (z) , where the integrals A dz e
R
go over states A and B. For analysis, we identify minimum free-energy pathways (MFEPs) between each minimum by performing path optimization on a free-energy landscape estimated via KDE. In contrast to standard string optimization algorithms, 59,60 rather than computing forces on-the-fly from MD, we keep the free-energy landscape fixed. The force on each interior path point com(t+1)
bines perpendicular gradient descent and spring forces zi
(t)
(t)
= zi + α −∇⊥ F (zi ) + κFs ,
(t)
where zi is the position of the i-th path point at iteration t. The force combines two terms: ∇⊥ F (zi ) is the perpendicular component of the free-energy gradient and the spring force (t)
(t)
(t)
Fs = zi+1 − 2zi + zi−1 with spring constant κ = 0.1 that maintains equal spacing between adjacent points. We use a learning rate α = 0.05. Around each path point, we select 10 closest samples to determine structural changes along each MFEP.
4
Results and Discussion
4.1
Dynamics of the Ligand-Free Pseudoknot Topologies
Before analyzing the ligand-containing trajectories, we first characterized the dynamics of ligand-free RT and RU pseudoknots. Both RT and RU remained globally folded throughout 16
the simulation (Fig. 4a). Despite the topological difference associated with 5′ -end threading, the two models exhibited a similar pattern of residue fluctuations. In particular, the largest fluctuations occur in the S3 (purple) and L2 (gray) regions, whereas the S1 (green) and S2 (orange) regions remain relatively stable (Fig. 4b–c), consistent with previous studies. 17 The high flexibility of S3 is further supported by the fraction of native contacts, where QS3 decreases to ∼ 0.2 in both RT and RU (Fig. 4f), while QS1 and QS2 remain closer to 1 (Fig. 4d–e). Below, we treat these ligand-free trajectories as the reference against which ligand-induced perturbations were evaluated.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 4: Comparison of structural stability and conformational dynamics of the ligandfree threaded (RT) and unthreaded (RU) SARS-CoV-2 pseudoknot models. (a) RMSD of the overall pseudoknot structure. (b,c) RMSF values for RT and RU, respectively. Colors denote the pseudoknot regions: 5′ - and 3′ -end tails (white), S1 (green), S2 (orange), S3 (purple), and L1–L2 (gray). (d–f) Fraction of native contacts, Q, calculated for S1, S2, and S3, respectively.
17
4.2
Ligand-Induced Destabilization Is Topology-Selective
To assess the influence of ligand binding on pseudoknot stability, we performed MD simulations of duration 0.4–1 µs for RT and RU forms while adding each ligand shown in Fig. 3. For each RNA–ligand pair, three separate MD simulations were initiated from the three top-ranked docking poses shown in Fig. S8 to assess the sensitivity of the observed RNA– ligand dynamics to the initial ligand placement and to test the reproducibility of ligandinduced structural responses. Across the full dataset of trajectories (Fig. S9–Fig. S13), most ligand-containing simulations remained structurally similar to the corresponding ligand-free trajectories, exhibiting only modest changes. However, a subset of systems demonstrated clear ligand-induced distortions of the pseudoknot conformation. The most pronounced effects were observed for RT-MERA0 (top-1 pose), RT-MERA± (top-1 pose), RT-A10 (top-1 pose), RU-MERA± (top-2 pose), and RU-A10 (top-1 pose); the pose number in parentheses identifies the corresponding docking-derived initial ligand position. To further characterize their long-timescale behavior, these simulations were extended to 3 µs for RT-MERA0 , RT-MERA± , and RT-A10 , and to 2 µs for RU-MERA± and RU-A10 . Among the RT–ligand systems, MERA0 produced the most pronounced structural perturbation (Fig. 5a). This effect is localized primarily in the S2 region, which becomes markedly more flexible upon ligand binding (Fig. 5b). MERA± and A10 also enhance fluctuations in S2, although the effect is more moderate. In agreement with these observations, the fraction of native contacts for S2, QS2 , drops significantly in all RT–ligand trajectories (Fig. 5d), whereas QS1 and QS3 (Fig. 5c,e) remain higher and do not fall to zero. Importantly, the S2 stem is partially restored after ∼1500 ns in the RT-MERA± simulation, while in RT-MERA0 and RT-A10 it drops down to zero quickly and remains disrupted for the remainder of the simulation time. These results suggest that ligand binding may selectively weaken a structurally vulnerable S2 region of the threaded pseudoknot, thereby causing breakage of the secondary structure. For the unthreaded RNA–ligand systems, both RU-MERA± and RU-A10 exhibit pro18
(a)
(c)
(b)
(d)
(e)
Figure 5: Structural stability and conformational dynamics of the threaded (RT) topology of the SARS-CoV-2 pseudoknot for the RNA–ligand systems showing the largest conformational distortions. (a) RMSD of the overall pseudoknot structure for RT–ligands. (b) RMSF values for the RT–ligand trajectories. Colors denote the pseudoknot regions: 5′ - and 3′ -end tails (white), S1 (green), S2 (orange), S3 (purple), and L1–L2 (gray). (c)–(e) Q, calculated for S1, S2, and S3, respectively. nounced deviations in RMSD over the course of the MD simulations (Fig. 6a). While the threaded ligand-containing systems show the dominant perturbation localized in the S2 region, the unthreaded complexes, however, display substantial distortions primarily in the S1 and S3 regions (Fig. 6b). This behavior matches what is seen in the Q values: in the RU-MERA± trajectory, QS1 eventually decreases to zero, and the native contacts of S1 do not recover. Similarly, in RU-A10 , QS1 also shows a substantial decrease (Fig. 6c). S2 remains comparatively stable in all RU-containing systems, with QS2 staying close to that of the ligand-free RU reference (Fig. 6d). As in the threaded RT topology, S3 remains the most flexible stem region. QS3 decreases markedly during the simulation, although it does not collapse completely to zero (Fig. 6e). In summary, these results show that ligand-induced destabilization of the unthreaded pseudoknot primarily affects the S1 and S3 regions. This
19
differs from the S2-centered response seen in the threaded topology.
(a)
(c)
(b)
(d)
(e)
(b)
Figure 6: Structural stability and conformational dynamics of the unthreaded (RU) topology of the SARS-CoV-2 pseudoknot for the RNA–ligand systems showing the largest conformational distortions. (a) RMSD of the overall pseudoknot structure for RU–ligands. (b) Perresidue RMSF profile for the RU–ligand trajectories. Colors denote the pseudoknot regions: 5′ - and 3′ -end tails (white), S1 (green), S2 (orange), S3 (purple), and L1–L2 (gray). (c)–(e) Q, calculated separately for S1, S2, and S3, respectively. The RMSD, RMSF, and Q analyses show that substantial pseudoknot reorganization occurs only in a minority of the ligand-containing MD trajectories. Thus, pseudoknot destabilization is not a generic consequence of ligand proximity, but instead depends on a specific combination of ligand type, protonation state, and pseudoknot topology. The many nonperturbing trajectories therefore serve as an important internal control, indicating that the pronounced conformational changes observed in the representative systems reflect selective RNA—ligand recognition events rather than systematic simulation drift. To further analyze pseudoknot conformational ensembles in these systems, we examined their trajectories using SM and compared them with those of the ligand-free RT and RU references. For contrast, we also included the RT-A20 system, which did not show significant 20
structural changes, in the ML algorithm. Eigenvalues
0
50 100 Epoch
0.0
1 0 −1
0.5
λk
0.5 0.0
0
50 100 Epoch
0.0
0
50 100 Epoch
0 1 2 3 4 5 6
k
0.5
0.5 0.0
1.0 0
50 100 Epoch
(f)
1 0 −1
1.0 0
(d)
50 100 Epoch
0 1 2 3 4 5 6
k
0.5 0.0
1 0 −1
0.5
0
50 100 Epoch
0.0
RU 1.0
1 0 −1
RU-A10
L
1 0 −1
L
(e)
(c)
0.0
RT-MERA± 1.0
1.0
RT-MERA0
(g)
1 0 −1
RU-MERA± 1.0
0
50 100 Epoch
0 1 2 3 4 5 6
k
0.5 0.0
(h)
1 0 −1
L
0.5
(b)
L
1.0
L
1 0 −1
L
(a)
L
λk
1.0
RT-A20
L
RT-A10
RT
0
50 100 Epoch
0 1 2 3 4 5 6
k
Figure 7: Eigenvalues λk corresponding to the converged CVs after SM learning. Dots and error bars show means and standard deviations, respectively, calculated over the last 50 epochs. Insets show the training loss L during 100 epochs.
4.3
Slow Modes Distinguish Threaded and Unthreaded Topologies
The eigenspectra of the converged CVs reveal a fundamental dynamical distinction between the threaded and unthreaded pseudoknot topologies (Fig. 7). All RT systems show clear spectral gaps (∆λ = λk − λk−1 ), having the first relevant eigenvalue λ1 ≈ 0.5 and the second λ2 ≈ 0.4. After that, the remaining eigenvalues drop rapidly to 0. This indicates three well-separated slow modes and, correspondingly, three dominant metastable states: S0 , S1 , and S2 . In contrast, ligand-free RU shows a gradual eigenvalue decay with no clear spectral gap (Fig. 7f), and the cumulative spectral weight saturates slowly (Fig. S5). This indicates that, for RU, there are no distinct slow modes; instead, all conformational transitions are relatively fast, as the system is essentially a broad superstate composed of smaller states separated by negligible free-energy barriers. Perhaps most strikingly, introducing the merafloxacin zwitterion to the unthreaded pseudoknot (RU-MERA± ) recovers a clear spectral gap resembling that of the RT systems (Fig. 7h). This demonstrates that MERA± imposes 21
slow dynamical behavior akin to that of the RT systems on the otherwise spectrally featureless unthreaded pseudoknot topology, suggesting that the inhibitor has a strong effect on the dynamics of the unthreaded system.
4.4
Ligands Reshape the Pseudoknot Free-Energy Landscape
Based on FELs constructed from the learned CVs, we find that the ligands can significantly reshape the conformational dynamics of the pseudoknot (Fig. 8). In ligand-free RT, we observe three metastable basins in a symmetric triangular arrangement with energy barriers between them of around ∼12–15 kJ/mol (Fig. 8a). Additionally, the FEL allows the systems to transition between all states, following a cyclic S0 → S1 → S2 → S0 pathway, with no discernible intermediate. In this case, stem S2 fluctuates the most. Adding MERA0 to the threaded form significantly reshapes the RT FEL (Fig. 8b). RT-MERA0 is characterized by an asymmetric landscape with one dramatically lowered barrier (FS†1 →S2 ≈ 3.9 kJ/mol, see Tab. 2). This facilitates rapid distortions in the stem where MERA0 currently binds. These occur without major conformational changes in RT. Both S1 and S2 , which correspond to this fast transition, are the most disrupted conformations relative to ligand-free RT. In both states, the ligand is bound to the S3 stem, and in S2 , the ligand is also bound to the S2 stem (Fig. 8b). This can be linked to the overall S2 stem disruption upon ligand binding. In relatively undistorted RT-A20 (Fig. 8h), we can see one very high barrier (FS†0 →S1 ≈ 29.3 kJ/mol, see Tab. 2), kinetically trapping the system. This observation is consistent with A2 being experimentally inactive. Most ligand conformations in RT-A20 diffuse around the pseudoknot rather than binding at a specific location. The RNA states here are qualitatively similar to those of ligand-free RT, which underscores that A20 has the weakest effect on RT. In ligand-free RU (Fig. 8e), the system is located in a flat, featureless landscape, having barriers F † < 4 kJ/mol (see Tab. 2). However, as explained above, in RU-MERA± , we can see that the ligand transforms the RU landscape into three well-separated basins with substantial barriers (13–18 kJ/mol, see Tab. 2), dramatically slowing conformational dynamics (Fig. 8f). 22
Figure 8: Free-energy landscapes F (z) reconstructed in CV space z = (z1 , z2 ) learned by SM for each system. Each state is labeled by S0 , S1 , and S2 . The rotation of each landscape is set so that the state labels are ordered clockwise. The ensembles corresponding to the metastable states depict 10 structures closest to the minima (white dots).
23
Moreover, in RU ligand destabilization shifts to the S1 and S3 stem regions. Additionally, in 6 out of 8 analyzed systems, there is no direct S0 → S2 minimum free-energy pathway, and thus transitions must pass through S1 . Only ligand-free RT and RT-A20 support a direct cyclic path (Fig. 8a,h). This hub-like role of S1 suggests it may represent a structurally intermediate conformation that gates access between the two endpoint states, further increasing the time required for the endpoint states to interchange. Table 2: Free-energy differences (∆F ) and free-energy barriers (F † ) between the identified states (S). Missing barrier values (marked by “—”) indicate missing direct path from S0 to S2 , and a transition goes through an intermediate S1 (see Fig. 8). See Fig. S7, for a figure showing free-energy barriers along MFEPs for each system. System RT
RT-MERA0
RT-MERA±
RT-A10
RU
RU-MERA±
RU-A10
RT-A20
4.5
Transition S0 → S 1 S1 → S2 S0 → S2 S0 → S 1 S1 → S 2 S0 → S 2 S0 → S 1 S1 → S 2 S0 → S 2 S0 → S 1 S1 → S 2 S0 → S2 S0 → S 1 S 1 → S2 S0 → S2 S0 → S 1 S1 → S 2 S0 → S 2 S0 → S 1 S1 → S 2 S0 → S2 S0 → S 1 S1 → S2 S0 → S2
∆F [kJ/mol] -4.67 0.12 -2.22 4.50 -0.85 1.06 -1.76 0.52 0.43 3.24 -2.69 0.8 0.05 -3.48 -4.53 0.27 1.99 2.84 2.02 -1.08 0.41 -0.84 -1.65 2.08
F † [kJ/mol] 12.43 14.91 13.99 18.63 3.92 — 12.50 15.58 — 11.93 9.63 — 0.51 3.69 — 17.74 13.12 — 11.21 0.93 — 29.33 15.68 13.50
Kinetic Control Reproduces Experimental −1 PRF Inhibitory Activity
Interestingly, the ligand perturbation effect on pseudoknot dynamics learned by SM agrees with the experimental antiviral activity. The degree of restructuring of the conformational landscape follows the order MERA > A1 > A2. This matches the degree of experimentally observed −1 PRF inhibition efficacy reported by Sun et al. 23 CVs learned by SM provide 24
a mechanistic basis for this ranking. Effective inhibitors are those that most strongly reshape the FEL and alter the slow-mode spectrum of the pseudoknot. Additionally, kinetics, not thermodynamics, governs the conformational preferences of the RNA-ligand systems (Tab. 2). Free-energy differences between metastable states are generally small (|∆F | < 5 kJ/mol in most cases). Thus, the states are roughly isoenergetic. However, barriers vary over an order of magnitude, from F † ≈ 0.5 kJ/mol in ligand-free RU to F † ≈ 29 kJ/mol in RT-A20 . The drug mechanism is therefore primarily kinetic, as the ligands modulate transition rates rather than equilibrium populations.
4.6
Merafloxacin Protonation Dictates the Binding Mechanism
An interesting result is that significant differences exist between FELs corresponding to MERA0 and MERA± interacting with the pseudoknot (Fig. 8b,c). These two protonation states of the same inhibitor cause qualitatively different FELs for the same RNA topology. When MERA0 interacts with the threaded form of the pseudoknot, it creates an asymmetric landscape with a very low FS†1 →S2 barrier. In contrast, MERA± produces a more symmetric three-basin landscape with similar barriers and free-energy differences. For the unthreaded form, only the zwitterionic form MERA± significantly perturbs the FEL of RU. Thus, the physiological protonation state of MERA is essential for accurate prediction of RNA–drug interactions.
5
Conclusions
We have used spectral map (SM), a thermodynamics-driven ML method, to characterize the conformational dynamics of the SARS-CoV-2 −1 PRF pseudoknot in complex with merafloxacin and two related analogs, for two distinct RNA fold topologies and two ligand protonation states. By learning CVs that directly target slow dynamic modes of the RNA– ligand systems, we obtain free-energy landscapes (FELs) that link chemical structure of
25
ligands, their protonated states, and RNA topology to the rare conformational transitions that underlie pseudoknot function. We have shown that ligand-induced destabilization is topology-selective. In the threaded pseudoknot, the strongest distortion occurs in the S2 stem. In the unthreaded pseudoknot, distortion occurs in the S1 and S3 stems. This pattern matches the presence of a 5′ -threading cleft in the threaded topology, which is absent in the unthreaded one. It indicates that the same inhibitor can engage different RNA regions depending on how the pseudoknot is folded. The degree to which each ligand reshapes the learned FEL follows the order MERA > A1 > A2, reproducing the experimental ranking of anti-frameshifting and antiviral activity reported by Sun et al. 23 Importantly, this modulation is primarily kinetic rather than thermodynamic. Free-energy differences between metastable states span only a few kJ/mol. In contrast, barriers vary by more than an order of magnitude across the systems examined. Effective inhibitors are therefore those that most strongly modify transition rates, not those that stabilize a specific bound conformation. Furthermore, the neutral and zwitterionic forms of merafloxacin create qualitatively different landscapes for the same RNA topology. The zwitterion alone imposes a clear spectral gap on the otherwise featureless unthreaded pseudoknot. This demonstrates that the protonation state of the ligand should be considered in mechanistic models of RNA–ligand recognition. These findings support a topology- and pose-dependent view of ligand recognition by the SARS-CoV-2 pseudoknot and, more broadly, by structured RNA. In contrast to proteins, which typically present well-defined binding pockets, the RNA pseudoknot offers shallow, transient interaction sites: ligands initially placed in a docked pose frequently migrated during the simulation, and the same RNA–ligand pair did not always produce the same response across trajectories initiated from different poses. The perturbations we observe, therefore, should not be interpreted as a single reproducible bound state, but as a small number of transient binding poses that couple to structurally sensitive regions of the RNA. We view the simulations presented here as a mechanistic exploration of candidate microscopic
26
states rather than a one-to-one structural assignment of experimentally observed RNA– ligand complexes. From a methodological standpoint, our results demonstrate that SM can recover the slow modes governing pseudoknot conformational dynamics. The resulting FELs provide a quantitative, mechanistic bridge between the chemical structure of the ligand and rare conformational transitions of RNA. SM thus complements the growing set of ML methods developed specifically for RNA. 61
Acknowledgement M.I. acknowledges funding JSPS KAKENHI Grant Number JP24K23894 and computational resources provided by Research Institute for Information Technology at Kyushu University, including the HPCI System Research Project hp250160. J.R. acknowledges funding from the Ministry of Science and Higher Education in Poland and the National Science Center in Poland (Sonata 2021/43/D/ST4/00920, “Statistical Learning of Slow Collective Variables from Atomistic Simulations”).
References (1) Dinman,
J.
D.
Mechanisms
and
Implications
of
Programmed
Transla-
tional Frameshifting. Wiley Interdiscip. Rev. RNA 2012, 3, 661–673, DOI: https://doi.org/10.1002/wrna.1126. (2) Falese, J. P.; Donlic, A.; Hargrove, A. E. Targeting RNA with Small Molecules: From Fundamental Principles Towards the Clinic. Chem. Soc. Rev. 2021, 50, 2224–2243, DOI: https://doi.org/10.1039/D0CS01261K. (3) Childs-Disney, J. L.; Yang, X.; Gibaut, Q. M.; Tong, Y.; Batey, R. T.; Disney, M. D.
27
Targeting RNA Structures with Small Molecules. Nat. Rev. Drug Discov. 2022, 21, 736–762, DOI: https://doi.org/10.1038/s41573-022-00521-4. (4) Kovachka, S.; Panosetti, M.; Grimaldi, B.; Azoulay, S.; Di Giorgio, A.; Duca, M. Small Molecule Approaches to Targeting RNA. Nat. Rev. Chem. 2024, 8, 120–135, DOI: https://doi.org/10.1038/s41570-023-00569-9. (5) Almena Rodriguez, L.; Kallert, E.; Husmann, J.-Å.; Schaubruch, K.; Meisel, K. I. S.; Schwickert, M.; Hoba, S. N.; Heermann, R.; Kersten, C. Electrostatic Anchoring in RNA-Ligand Design—Dissecting the Effects of Positive Charges on Affinity, Selectivity, Binding Kinetics, and Thermodynamics. J. Med. Chem. 2025, 68, 8659–8678, DOI: https://doi.org/10.1021/acs.jmedchem.5c00339. (6) Finkel, Y.; Mizrahi, O.; Nachshon, A.; Weingarten-Gabbay, S.; Morgenstern, D.; Yahalom-Ronen, Y.;
Tamir, H.;
Achdout, H.;
Stein, D.;
Israeli, O.;
oth-
ers The Coding Capacity of SARS-CoV-2. Nature 2021, 589, 125–130, DOI: https://doi.org/10.1038/s41586-020-2739-1. (7) Jones, C. P.; Ferré-D’Amaré, A. R. Crystal Structure of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Frameshifting Pseudoknot. RNA 2022, 28, 239–249, DOI: https://doi.org/10.1261/rna.078825.121. (8) Roman, C.; Lewicka, A.; Koirala, D.; Li, N.-S.; Piccirilli, J. A. The SARS-CoV-2 Programmed −1 Ribosomal Frameshifting Element Crystal Structure Solved to 2.09 Å using Chaperone-Assisted RNA Crystallography. ACS Chem. Biol. 2021, 16, 1469– 1481, DOI: https://doi.org/10.1021/acschembio.1c00324. (9) Zhang, K.; Zheludev, I. N.; Hagey, R. J.; Haslecker, R.; Hou, Y. J.; Kretsch, R.; Pintilie, G. D.; Rangan, R.; Kladwang, W.; Li, S.; others Cryo-Em and Antisense Targeting of the 28-kDa Frameshift Stimulation Element From the SARS-
28
CoV-2 RNA Genome. Nat. Struct. & Mol. Biol. 2021, 28, 747–754, DOI: https://doi.org/10.1038/s41594-021-00653-y. (10) Bhatt, P. R.; Scaiola, A.; Loughran, G.; Leibundgut, M.; Kratzel, A.; Meurs, R.; Dreos, R.; O’Connor, K. M.; McMillan, A.; Bode, J. W.; others Structural Basis of Ribosomal Frameshifting during Translation of the SARS-CoV-2 RNA Genome. Science 2021, 372, 1306–1313, DOI: https://doi.org/10.1126/science.abf3546. (11) Peterson, J. M.;
Becker, S. T.;
O’Leary, C. A.;
Juneja, P.;
Yang, Y.;
Moss, W. N. Structure of the SARS-CoV-2 Frameshift Stimulatory Element With an Upstream Multibranch Loop. Biochemistry 2024, 63, 1287–1296, DOI: https://doi.org/10.1021/acs.biochem.3c00716. (12) Wacker, A.; Weigand, J. E.; Akabayov, S. R.; Altincekic, N.; Bains, J. K.; Banijamali, E.; Binas, O.; Castillo-Martinez, J.; Cetiner, E.; Ceylan, B.; others Secondary Structure Determination of Conserved SARS-CoV-2 RNA Elements by NMR Spectroscopy. Nucl. Acids Res. 2020, 48, 12415–12435, DOI: https://doi.org/10.1093/nar/gkaa1013. (13) Neupane, K.; Zhao, M.; Lyons, A.; Munshi, S.; Ileperuma, S. M.; Ritchie, D. B.; Hoffer, N. Q.; Narayan, A.; Woodside, M. T. Structural Dynamics of Single SARSCoV-2 Pseudoknot Molecules Reveal Topologically Distinct Conformers. Nat. Commun. 2021, 12, 4749, DOI: https://doi.org/10.1038/s41467-021-25085-6. (14) Pekarek, L.;
Zimmer, M. M.;
Gribling-Burrer, A.-S.;
Buck, S.;
Smyth, R.;
Caliskan, N. Cis-Mediated Interactions of the SARS-CoV-2 Frameshift RNA Alter its Conformations and Affect Function. Nucl. Acids Res. 2023, 51, 728–743, DOI: https://doi.org/10.1093/nar/gkac1184. (15) Omar, S. I.; Zhao, M.; Sekar, R. V.; Arbabimoghadam, S.; Tuszynski, J. A.; Woodside, M. T. Modeling the Structure of the Frameshift-Stimulatory Pseudoknot in SARS29
CoV-2 Reveals Multiple Possible Conformers. Biophys. J. 2021, 120, 313a–314a, DOI: https://doi.org/10.1016/j.bpj.2020.11.1986. (16) Schlick, T.; Zhu, Q.; Jain, S.; Yan, S. Structure-Altering Mutations of the SARSCoV-2 Frameshifting RNA Element. Biophys. J. 2021, 120, 1040–1053, DOI: https://doi.org/10.1016/j.bpj.2020.10.012. (17) Yan, S.; Zhu, Q.; Jain, S.; Schlick, T. Length-Dependent Motions of SARSCoV-2 Frameshifting RNA Pseudoknot and Alternative Conformations Suggest Avenues for Frameshifting Suppression. Nat. Commun. 2022, 13, 4284, DOI: https://doi.org/10.1038/s41467-022-31353-w. (18) He, W.;
San Emeterio, J.;
Woodside, M. T.;
Kirmizialtin, S.;
Pollack, L.
Atomistic Structure of the SARS-CoV-2 Pseudoknot in Solution from SAXSDriven Molecular Dynamics. Nucl. Acids Res. 2023, 51, 11332–11344, DOI: https://doi.org/10.1093/nar/gkad809. (19) Park, S.-J.; Kim, Y.-G.; Park, H.-J. Identification of RNA Pseudoknot-Binding Ligand That Inhibits the −1 Ribosomal Frameshifting of SARS-coronavirus by Structure-Based Virtual Screening. J. Am. Chem. Soc. 2011, 133, 10094–10100, DOI: https://doi.org/10.1021/ja1098325. (20) Kelly, J. A.; Olson, A. N.; Neupane, K.; Munshi, S.; San Emeterio, J.; Pollack, L.; Woodside, M. T.; Dinman, J. D. Structural and Functional Conservation of the Programmed −1 Ribosomal Frameshift Signal of SARS Coronavirus 2 (SARS-CoV-2). J. Biol. Chem. 2020,
295,
10741–10748,
DOI:
https://doi.org/10.1074/jbc.AC120.013449. (21) Neupane, K.; Munshi, S.; Zhao, M.; Ritchie, D. B.; Ileperuma, S. M.; Woodside, M. T. Anti-Frameshifting Ligand Active against SARS Coronavirus-2 is Resistant to Natural
30
Mutations of the Frameshift-Stimulatory Pseudoknot. J. Mol. Biol. 2020, 432, 5843– 5847, DOI: https://doi.org/10.1016/j.bpj.2020.11.1838. (22) Munshi, S.; Neupane, K.; Ileperuma, S. M.; Halma, M. T.; Kelly, J. A.; Halpern, C. F.; Dinman, J. D.; Loerch, S.; Woodside, M. T. Identifying Inhibitors of −1 Programmed Ribosomal Frameshifting in a Broad Spectrum of Coronaviruses. Viruses 2022, 14, 177, DOI: https://doi.org/10.3390/v14020177. (23) Sun, Y.; Abriola, L.; Niederer, R. O.; Pedersen, S. F.; Alfajaro, M. M.; Silva Monteiro, V.; Wilen, C. B.; Ho, Y.-C.; Gilbert, W. V.; Surovtseva, Y. V.; others Restriction of SARS-CoV-2 Replication by Targeting Programmed −1 Ribosomal Frameshifting. Proc. Natl. Acad. Sci. U.S.A. 2021, 118, e2023051118, DOI: https://doi.org/10.1073/pnas.202305111. (24) Zafferani, M.; Muralidharan, D.; Montalvan, N. I.; Hargrove, A. E. RT-qPCR as a Screening Platform for Mutational and Small Molecule Impacts on Structural Stability of RNA Tertiary Structures. RSC Chem. Biol. 2022, 3, 905–915, DOI: https://doi.org/10.1039/D2CB00015F. (25) Varricchio, C.; Mathez, G.; Pillonel, T.; Bertelli, C.; Kaiser, L.; Tapparel, C.; Brancale, A.; Cagno, V. Geneticin Shows Selective Antiviral Activity against SARS-CoV-2 by Interfering with Programmed −1 Ribosomal Frameshifting. Antiviral Res. 2022, 208, 105452 – 105452, DOI: https://doi.org/10.1016/j.antiviral.2022.105452. (26) Chen, Y.; Tao, H.; Shen, S.; Miao, Z.; Li, L.; Jia, Y.; Zhang, H.; Bai, X.; Fu, X. A Drug Screening Toolkit based on the −1 Ribosomal Frameshifting of SARS-CoV-2. Heliyon 2020, 6, DOI: https://doi.org/10.1016/j.heliyon.2020.e04793. (27) Ahn, D.-G.; Yoon, G. Y.; Lee, S.; Ku, K. B.; Kim, C.; Kim, K.-D.; Kwon, Y.C.; Kim, G.-W.; Kim, B.-T.; Kim, S.-J. A Novel Frameshifting Inhibitor Having
31
Antiviral Activity against Zoonotic Coronaviruses. Viruses 2021, 13, 1639, DOI: https://doi.org/10.3390/v13081639. (28) Mathez, G.; Brancale, A.; Cagno, V. Novel Inhibitors of SARS-CoV-2 RNA Identified through Virtual Screening. J. Chem. Inf. Model. 2024, 64, 6190–6196, DOI: https://doi.org/10.1021/acs.jcim.4c00758. (29) Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing Important Fluctuations: Rare Events and Metadynamics From a Conceptual Viewpoint. Annu. Rev. Phys. Chem. 2016, 67, 159–184, DOI: https://doi.org/10.1146/annurev-physchem-040215-112229. (30) Bussi, Energy
G.;
Laio,
Landscapes.
A.
Using Nat.
Metadynamics
Rev.
Phys.
to
Explore
2020,
2,
Complex
Free-
200–2012,
DOI:
https://doi.org/10.1038/s42254-020-0153-0. (31) Sega,
M.;
Landscapes.
Salvalaglio, Annu.
Rev.
M.
Molecular
Chem.
Understanding
Biomol.
Eng.
of
2026,
Free-Energy 17,
DOI:
https://doi.org/10.1146/annurev-chembioeng-100724-082451. (32) Van Doorslaer, X.; Dewulf, J.; Van Langenhove, H.; Demeestere, K. Fluoroquinolone Antibiotics: An Emerging Class of Environmental Micropollutants. Sci. Total Environ. 2014, 500, 250–269, DOI: https://doi.org/10.1016/j.scitotenv.2014.08.075. (33) Sponer, J.; Bussi, G.; Krepl, M.; Banas, P.; Bottaro, S.; Cunha, R. A.; Gil-Ley, A.; Pinamonti, G.; Poblete, S.; Jurecka, P.; others RNA Structural Dynamics as Captured by Molecular Simulations: A Comprehensive Overview. Chem. Rev. 2018, 118, 4177– 4338, DOI: https://doi.org/10.1021/acs.chemrev.7b00427. (34) Mlynsky, V.; Bussi, G. Exploring RNA Structure and Dynamics Through Enhanced Sampling Simulations. Curr. Opin. Struct. Biol. 2018, 49, 63–71, DOI: https://doi.org/10.1016/j.sbi.2018.01.004.
32
(35) Bernetti, M.; Bussi, G. Integrating Experimental Data with Molecular Simulations to Investigate RNA Structural Dynamics. Curr. Opin. Struct. Biol. 2023, 78, 102503, DOI: https://doi.org/10.1016/j.sbi.2022.102503. (36) Languin-Cattoën, O.;
Bussi, G. RNA Dynamics and Interactions Revealed
Through Atomistic Simulations. Annu. Rev. Phys. Chem. 2026,
77,
DOI:
https://doi.org/10.1146/annurev-physchem-082624-013453. (37) Peters, sis
Tests.
B.
Reaction Annu.
Rev.
Coordinates
and
Mechanistic
Phys.
2016,
67,
Chem.
Hypothe-
669–690,
DOI:
https://doi.org/10.1146/annurev-physchem-040215-112215. (38) Rogal, J. Reaction Coordinates in Complex Systems – A Perspective. Eur. Phys. J. B 2021, 94, 1–9, DOI: https://doi.org/10.1140/epjb/s10051-021-00233-5. (39) Rydzewski, J.; Chen, M.; Valsson, O. Manifold Learning in Atomistic Simulations: A Conceptual Review. Mach. Learn.: Sci. Technol. 2023, 4, 031001, DOI: https://doi.org/10.1088/2632-2153/ace81a. (40) Gokdemir, T.; Rydzewski, J. Machine Learning of Slow Collective Variables and Enhanced Sampling via Spatial Techniques. Chem. Phys. Rev. 2025, 6, 011304, DOI: https://doi.org/10.1063/5.0245177. (41) Zhu, K.; Trizio, E.; Zhang, J.; Hu, R.; Jiang, L.; Hou, T.; Bonati, L. Enhanced Sampling in the Age of Machine Learning: Algorithms and Applications. Chem. Rev. 2025, 126, 671–713, DOI: https://doi.org/10.1021/acs.chemrev.5c00700. (42) Rydzewski, tive
J.
Variables.
Spectral J.
Phys.
Map:
Embedding
Chem.
Lett.
2023,
https://doi.org/10.1021/acs.jpclett.3c01101.
33
Slow 14,
Kinetics
in
Collec-
5216–5220,
DOI:
(43) Rydzewski, J.; Gokdemir, T. Learning Markovian Dynamics With Spectral Maps. J. Chem. Phys. 2024, 160, 091102, DOI: https://doi.org/10.1063/5.0189241. (44) Rydzewski, J. Spectral Map for Slow Collective Variables, Markovian Dynamics, and Transition State Ensembles. J. Chem. Theory Comput. 2024, 20, 7775–7784, DOI: https://doi.org/10.1021/acs.jctc.4c00428. (45) Coifman, R. R.; Lafon, S. Diffusion Maps. Appl. Comput. Harmon. Anal. 2006, 21, 5–30, DOI: https://doi.org/10.1016/j.acha.2006.04.006. (46) Nadler, B.; Lafon, S.; Coifman, R. R.; Kevrekidis, I. G. Diffusion Maps, Spectral Clustering and Reaction Coordinates of Dynamical Systems. Appl. Comput. Harmon. Anal. 2006, 21, 113–127, DOI: https://doi.org/10.1016/j.acha.2005.07.004. (47) Lezcano
Casado,
M.
Trivializations
tion
Manifolds.
Adv.
Neural
on
Inf.
for Process.
Gradient-Based Syst.
2019,
Optimiza32,
DOI:
https://doi.org/10.48550/arXiv.1909.09501. (48) Yang, M.; Olatunji, F. P.; Rhodes, C.; Balaratnam, S.; Dunne-Dombrink, K.; Seshadri, S.; Liang, X.; Jones, C. P.; Le Grice, S. F.; Ferre-D’Amare, A. R.; others Discovery of Small Molecules Targeting the Frameshifting Element RNA in SARS-CoV-2 Viral Genome. ACS Med. Chem. Lett. 2023, 14, 757–765, DOI: https://doi.org/10.1021/acsmedchemlett.3c00051. (49) Trott, O.; Olson, A. J. AutoDock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J. Comput. Chem. 2010, 31, 455–461, DOI: https://doi.org/10.1002/jcc.21334. (50) Eberhardt, J.; Santos-Martins, D.; Tillack, A. F.; Forli, S. AutoDock Vina 1.2. 0: New Docking Methods, Expanded Force Field, and Python Bindings. J. Chem. Inf. Model. 2021, 61, 3891–3898, DOI: https://doi.org/10.1021/acs.jcim.1c00203.
34
(51) Quiroga, R.; Villarreal, M. A. Vinardo: A Scoring Function Based on Autodock Vina Improves Scoring, Docking, and Virtual Screening. PloS One 2016, 11, e0155183, DOI: https://doi.org/10.1371/journal.pone.0155183. (52) Denning, E. J.; Priyakumar, U. D.; Nilsson, L.; Mackerell Jr, A. D. Impact of 2‘hydroxyl Sampling on the Conformational Properties of RNA: Update of the CHARMM All-Atom Additive Force Field for RNA. J. Comput. Chem. 2011, 32, 1929–1943, DOI: https://doi.org/10.1002/jcc.21777. (53) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; others CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM AllAtom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671–690, DOI: https://doi.org/10.1002/jcc.21367. (54) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101, DOI: https://doi.org/10.1063/1.2408420. (55) Bernetti, M.; Bussi, G. Pressure Control Using Stochastic Cell Rescaling. J. Chem. Phys. 2020, 153, 114107, DOI: https://doi.org/10.1063/5.0020514. (56) McGibbon, R. T.; Beauchamp, K. A.; Harrigan, M. P.; Klein, C.; Swails, J. M.; Hernández, C. X.; Schwantes, C. R.; Wang, L.-P.; Lane, T. J.; Pande, V. S. MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophys. J. 2015, 109, 1528 – 1532, DOI: https://doi.org/10.1016/j.bpj.2015.08.015. (57) Best, R. B.; Hummer, G.; Eaton, W. A. Native Contacts Determine Protein Folding Mechanisms in Atomistic Simulations. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 17874– 17879, DOI: https://doi.org/10.1073/pnas.1311599110. (58) Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.; Chanan, G.; Killeen, T.;
35
Lin, Z.; Gimelshein, N.; Antiga, L.; others PyTorch: An Imperative Style, HighPerformance Deep Learning Library . Adv. Neural Inf. Process. Syst. 33, 8026–8037. (59) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces. J. Chem. Phys. 2006, 125, 024106, DOI: https://doi.org/10.1063/1.2212942. (60) Maragliano, for
Sampling
Events
L.;
Vanden-Eijnden,
Free
Energy
Simulations.
Chem.
and
E.
A
Temperature
Determining
Phys.
Lett.
Accelerated
Method
Pathways
in
Rare
168–175,
DOI:
Reaction
2006,
426,
https://doi.org/10.1016/j.cplett.2006.05.062. (61) Sacco, G.; Bussi, G.; Sanguinetti, G. Machine Learning for RNA Secondary Structure Prediction: A Review of Current Methods and Challenges. RNA 2026, 32, 443–456, DOI: https://doi.org/10.1261/rna.080840.125.
36