Exploring the conformation-dependent reactivity and dynamics of a dinucleotide inhibitor of ribonuclease A

Abstract

The dinucleotide molecule, 5′-phospho-2′-deoxyuridine-3′-pyrophosphate adenosine-3′-phosphate (pdUppA-3′-p) (PUA), is known to act as a highly potent inhibitor of bovine pancreatic ribonuclease A (RNase A) with K_i = 27 nM (D. D. Leonidas et al., Biochemistry, 1999, 38, 10287–10297). In the present work, PUA is investigated to study the effect of different conformations on its residence time and dynamical interactions at the active site of RNase A. Two conformers of PUA, designated as closed and open in terms of their end-to-end distances, have been screened in different environments based on a wide range of global and local reactivity parameters computed using conceptual density functional theory (CDFT). In all the cases studied, the open conformation is found to be more flexible with its reactive sites predominantly localized at the terminal uracil nucleobase. The closed conformer, on the other hand, is associated more with intramolecular hydrogen bonding and van der Waals interactions at both nucleobase termini. When subjected to classical molecular dynamics (MD) simulation, the open conformer is found to dissociate rapidly from the active site of RNase A within 24–300 ns, while the closed conformer remained bound throughout the 5 µs long MD trajectory. The transition from open to closed conformation is found to stabilize the PUA–RNase A complex, lowering its free energy by nearly 25 kcal mol⁻¹. The crystal structure of the PUA–RNase A complex reports the ligand molecule only in its open conformation at the active site of RNase A (N. Russo and R. Shapiro, J. Biol. Chem., 1999, 274, 14902–14908). The closed conformer has not been detected experimentally as yet. Our results therefore indicate a hereto unexplored role of the closed conformation in the inhibition of RNase A by PUA.

---

1 Introduction

Synthetic enzyme inhibitors are often designed based on the structure and interactive mechanism of the substrate with the active site residues of the enzyme.^1–3 To inhibit a target enzyme, new molecules may be designed or existing ones repurposed by exploring the available databases, including known drugs, pro-drugs and many other small, organic, drug-like molecules. The lead identification of a therapeutically suitable candidate from the databases still remains a challenge due to the rigorous requirements of potency, selectivity and reactivity within the complex biological environment. Additional complications in the search and/or design may arise, for instance, due to a non-trivial role of the dynamics of protein–ligand conjugates sampling multiple binding poses with varying residence times at the binding site(s). These dynamical factors compete with the average, time-independent principles usually employed to optimize the protein–ligand interaction while searching the databases. It is therefore imperative that one includes conformational diversity and the duration of binding of the inhibitor molecule in an effective design principle of an inhibitor.

In this article, we have investigated a synthetic di-nucleotide molecule, 5′-phospho-2′-deoxyuridine-3′-pyrophosphate adenosine-3′-phosphate (pdUppA-3′-p) (PUA),⁴ as shown in Fig. 1(a). It is a well-known inhibitor of the enzyme bovine pancreatic ribonuclease A (RNase A) with a K_i value of 27 nM.⁵ It is characterized by three important anchoring regions through the uracil, α-phosphate, and adenine moieties. The crystal structure of PUA has been resolved when bound to the active site of RNase A (PDB id: 1QHC)⁵ and is presented in Fig. 1(b). The adenine group of PUA has been shown to participate in significant hydrogen bonding and stacking interactions with Asn-67, Gln-69, Asn-71, and His-119 in the active site of RNase A while adopting a syn-conformation. The uridine 5′-phosphate contacts Thr-45 and Lys-66 while engaging other active site residues like Gln-11, His-12, His-119 and Phe-120 in its proximity. The inhibitor possesses a robust binding affinity not only for RNase A but also to RNase 4 and eosinophil-derived neurotoxin (EDN) (K_i = 180 and 260 nM, respectively).⁴

Fig. 1 Schematic illustration of the inhibitor PUA (a) and the crystal structure of the RNase A and PUA complex (b) (obtained from PDB id: 1QHC).⁵

The crystal structure of the PUA–RNase A complex, as shown in Fig. 1(b), reveals the conformation of PUA with its 3′ and 5′ ends located far apart from each other. No other stable conformation of PUA, when bound to the active site of RNase A, has been reported so far. The reliance of RNase A on the conformational selection of inhibitors has been noted in earlier structural and simulation studies.^6,7 2D-NMR and X-ray crystallographic data reveal that inhibitors, such as inosine 5′-phosphate (IMP, K_i = 4.6 ± 0.2 mM)⁸ when bound to the active site, are stabilized by a network of hydrogen bonds and electrostatic contacts tightly constraining the base orientation and ribose puckering. Similarly, 3′-CMP (K_i = 58 µM) bound at the active site of RNase A⁹ is reported to have a well-defined conformation enforced by an anti-glycosidic angle and a C3′-endo sugar pucker, even though the ligand samples multiple conformations in solution. The pyrophosphate-linked dinucleotides have been reported to bind when they take on compact, folded conformations while maintaining a syn orientation of the adenosine ring.^5,10

The results cited above do not reveal any direct correlation between conformational selection and high inhibitory action. In this article, we undertake a detailed computational study to explore different conformations of PUA and understand their difference in reactivity and dynamical interactions at the active site of RNase A. In principle, the flexible structure of PUA may adopt multiple stable conformations at the active site of RNase A under physiological conditions. An earlier classical molecular dynamics (MD) simulation study of the PUA–RNase A complex¹¹ until 4 ns revealed a conformation of PUA with a small and negligible deviation from that observed in its crystal structure. The stability of the complex has been attributed in this study primarily to electrostatic interactions involving the amino acid residues at and around the active site. However, this MD trajectory is too short to capture any transition to other metastable conformations of PUA within the active site. It is also possible that other conformations of PUA may not have been detected under the experimental conditions employed to determine the structure of the complex. In the absence of any dynamical studies, the implications of the conformational diversity of PUA for its putative inhibitory mechanism appear to be far from being completely understood.

To understand the interplay between the conformation and inhibitory action of any flexible molecule such as PUA, a complete set of functionally distinct, metastable states needs to be identified by searching the high-dimensional space spanning the conformations of the inhibitor molecule and the dynamic protein matrix hosting it. With the latest advancements in enhanced sampling combined with classical molecular dynamics simulations, it is indeed possible to identify a wide range of metastable conformational states of a flexible molecule along with the underlying free energy landscape and kinetics of inter-state transitions. However, a quantitative understanding of the function of these states with respect to the molecular mechanism of inhibition would warrant the validation of simulated atomistic models against the macroscopic, experimental values of K_i. This would involve prohibitively high computational costs associated with the required quantum mechanical (QM) or hybrid quantum mechanical–molecular mechanical (QM–MM) simulations interfaced with enhanced sampling methods. In this article, our major aim is to devise an alternative computational method whereby the search of the functionally relevant conformational space may be carried out in a reduced dimension by exploring a small number of collective variables (CVs, such as pair distances, dihedral angles, etc.), capable of delineating different metastable states of the inhibitor-protein complex. In the present article, we demonstrate how a reactivity-based screening of a single CV leads to the selection of two conformations of the ligand molecule PUA that eventually yield distinct residence times and dynamical interactions at the enzyme active site.

The large distance between the 3′ and 5′ ends of PUA in the crystal structure of the PUA–RNase A complex, as shown in Fig. 1(b), encouraged us to use it as the CV to characterize different conformations of PUA. As highlighted in Fig. 1(a), d_ee is defined as the distance between the atoms of adenine and of uracil moieties. We initiate the present study by empirically designating the PUA conformations as closed (when d_ee < 10 Å) and as open otherwise. Evidently, PUA adopts an open conformation in the crystal structure with d_ee = 13.59 Å. We generate an initial closed structure of PUA in silico by twisting both the termini of the open conformer and reducing d_ee to 3.18 Å. These initial structures are presented in Fig. 2(a) and (b).

Fig. 2 (a) and (b) Initial and (c) and (d) optimized geometries of the closed (a) and (c) and open (b) and (d) conformations of the inhibitor PUA in the gas phase. The atomic index numbering for the optimized structures is included in (c) and (d). Also shown are the end-to-end distance, d_ee (purple, dashed line) and the intramolecular H-bonds (orange, dotted line).

In view of an empirical choice of a single collective variable d_ee to generate different conformations of PUA, we next examine if the chosen conformations may be discriminated in terms of their reactivity in the gas phase and in polar and non-polar environments. We also investigate the effect of confinement on the reactivity of the chosen conformations of PUA in the presence of polar and non-polar amino acid residues, mimicking the active site of RNase A. The use of quantum chemical studies based on conventional density functional theory (DFT) and conceptual DFT on the chosen pair of conformations of PUA reveals distinct reactivity profiles and a preferential stabilization of the closed conformer in different media. We subsequently docked each conformation of PUA at the active site of RNase A using classical molecular mechanical force fields. The resulting distribution of molecular electrostatic potentials mapped onto the van der Waals surface of PUA qualitatively preserved the distinct distribution of electron-rich and electron-deficient terminal regions of the two conformers observed in our DFT-based calculations.

We have finally performed classical molecular dynamics (MD) simulations to probe the residence times of these two conformers at the active site of RNase A in water with the goal of identifying the key interactions stabilizing the protein–ligand conjugate in each case. Our results indicate that the closed conformation of PUA remains bound to the active site pocket until 5 µs while the open conformer unbinds from the active site within 300 ns. Therefore, the chosen models provide us with two states of the inhibitor that differ not only in terms of their conformations, but also with respect to their dynamical interactions at the binding site. Further analysis of the free energy profile of the unbinding of each conformation sheds light on a mechanism of inhibition of RNase A by PUA that has not been explored so far.

2 Methods

2.1 Geometry optimization of the initial structures

In this work, the initial molecular geometries of the inhibitors were built using the Gauss View 6.0 interface and the hydrogens were included using OpenBabel¹² and were kept fixed for all calculations.¹³ All quantum mechanical (QM) optimizations were performed first at the Hartree–Fock (HF) level and then re-optimized at the density functional theory (DFT) level using the B3LYP functional in conjunction with the 6-31++G(d,p) and def2-TZVP basis sets as implemented in the Gaussian suite to qualitatively assess the stability of the inhibitor. Frequency calculations at the same level of theory confirmed the absence of imaginary modes, indicating that the structures are true minima.

2.2 Reactivity descriptors using conceptual density functional theory (CDFT)

The wavefunction-based description of any molecule relies on representing the ground electronic state in terms of a 4N-dimensional complete wavefunction Ψ(x⃑^N) where the i-th electron (i = 1, N) is described in terms of the four-dimensional vector x⃑ comprised of three spatial and one spin coordinate. Conceptual density functional theory (CDFT), extended from DFT^14–16 simplifies the representation by using the single-particle electron density ρ(r⃑), a function of only three variables, r⃑: (x,y,z). In the absence of any other external field, CDFT further quantifies the well-known chemical concepts, such as electronegativity and hardness, in terms of the derivatives of the total energy (E) with respect to the total number of electrons, N and the potential, v(r⃑) felt by the electrons due to the nuclei. The foundation and applications of CDFT have been extensively reviewed in the literature.¹⁷ We summarize below the global¹⁵ and local¹⁸ reactivity descriptors utilized in the present work.

2.2.1 Global reactivity descriptors. The global (r⃑ independent) reactivity descriptors^15,17,19,20 have been investigated with the input of an accurate estimate of the HOMO and LUMO orbital energies, E. Among these descriptors are the global softness (S), global hardness (η), electronegativity (χ), chemical potential (μ), the HOMO–LUMO energy gap (E_LUMO–E_HOMO) and the global electrophilicity index (ω).^21–26 These are defined as follows.

The concept of electronegativity here represents the tendency of an atom or molecule to attract electrons in a chemical bond and is the average of ionization potential (I) and electron affinity (A).²³

Ionization potential (I) is the energy required to remove an electron from the molecule, which can be approximated as E_HOMO. Electron affinity (A) is the energy released when an electron is added to the molecule and approximated as E_LUMO, which is indicated by eqn (6).

As defined by Pauling, chemical hardness is the resistance of the chemical potential to a change in the number of electrons, which could be correlated with the stability and reactivity of a chemical system.²³ Moreover, the global hardness (η) and softness (S) are crucial while measuring molecular stability and reactivity. Low values of μ and ω indicate an excellent nucleophile, whereas high values of these parameters are indicative of good electrophiles. Depending on which molecule has a higher (lower) electrophilicity index, one will react as an electrophile or nucleophile. The electrophilicity index, according to Parr et al., was developed to quantify the energy reduction brought about by the electron transfer between the donor and the acceptor.^19,22

2.2.2 Local reactivity. Local reactivity descriptors are r⃑-dependent measures providing site-specific information about the electronic and reactive behavior of molecules^{14,17,19,27,28} such as electrophilicity or nucleophilicity.²⁹ Among the known local reactivity descriptors, the Fukui functions¹⁸ f_k measure the change in electron density at a specific atomic site due to the addition or removal of electrons. While the nucleophilic Fukui function, defined as

f⁺_k = ρ_N+1(r) − ρ_N(r) (7)

describes the response to the addition of an electron, the electrophilic Fukui function registers the effect of the removal of an electron as,

f⁻_k = ρ_N(r) − ρ_N−1(r) (8)

The radical Fukui function (average of f⁺_k and f⁻_k),

is also found to be useful in analyzing the local reactivity profiles.

The local softness associated with the corresponding Fukui function can be defined as,

s⁺_k = S·f⁺_k, s⁻_k = S·f⁻_k, s⁰_k = S·f⁰_k (10)

The local softness is provided with the reactivity of the atomic sites obtained from the product of the Fukui function and global softness (S). The local behavior of the molecule that is vulnerable to electrophilic (nucleophilic) attack is the source of its global tendency.^28,30 By addressing the identity related to the normalization of the Fukui function, a generalized treatment of nucleophilicity and both global and local electrophilicity has been discussed at length.²⁹

2.2.3 Non-covalent interactions. Non-covalent interactions (NCIs) are approximated from the reduced density gradient (RDG) method, providing a detailed map of weaker interactions in the molecule.³¹ These interactions that influence the binding efficiency of the inhibitor include hydrogen bonds, van der Waals forces, and π–π stacking, etc. The RDG function is essentially a dimensionless form of the electron density gradient norm function and is defined aswhere is the gradient of the electron density ρ(r⃑) at a point r⃑. Weak interaction regions will remain and be visible only after identifying the regions where ρ(r⃑) is very small. λ₂, the second eigenvalue of the Hessian matrix of the electron density is utilized to distinguish between repulsive and attractive interactions as sign(λ₂). Attractive forces are characterized by negative λ₂, for hydrogen-bonded regions. Repulsive forces are characterized by positive λ₂ values and van der Waals forces are identified when λ₂ ≈ 0. The RDG scatter plot is thus useful in predicting dense interaction networks and the sensitive non-covalent space.

2.2.4 Details of quantum chemical computation. The reactivity parameters, as listed above, have been computed not only in the gas phase, but also in polar (water) and non-polar (diethyl ether) solvent phases. In each system, the global reactivity descriptors are calculated from the total energy, E in gas and implicit solvent using the geometry optimized at the B3LYP/def2-TZVP level, including DFT-D3 dispersion corrections. The same computational level was maintained due to its reliable performance in capturing electron density and orbital energy characteristics. BJ-corrected DFT-D3 was excluded, as it yielded no significant difference in the total energy (Table S1). HOMO–LUMO gaps without the dispersion terms are provided in Fig. S1. The change in relative energy RE has been computed as ΔE_RE = [E_solvent − E_gas] to quantify the effect of the solvent. All optimizations have been conducted using the Gaussian 16 software suite¹³ and analyzed using Multiwfn^32,33 to study the conceptual DFT, weak interactions and electrostatic maps.

2.3 Classical molecular dynamics simulation

The inhibitor models (open and closed), both refined at the def2-TZVP level of theory, were subsequently docked into the active site of RNase A using AutoDock Vina,^34,35 which employs a gradient-based conformational search algorithm using a hybrid model of empirical scoring functions that accounts for steric, hydrophobic, hydrogen bonding, and repulsive interactions to estimate binding affinity. To proceed with docking, the protonation states of the ionizable residues were assigned using H++³⁶ at physiological pH and the residues at and near the active site were examined manually to ensure the chemically plausible protonation states validating the experimentally proposed catalytic mechanism.³⁷ For each model, the most favorable binding pose was selected and the docked enzyme–inhibitor complex subjected to classical molecular dynamics (MD) simulation in an explicit water environment using the AMBER24 software package.³⁸

The simulation system has been parameterized with the ff19SB protein force field for the enzyme and general Amber force field 2 (GAFF2)³⁹ for the inhibitor. Partial atomic charges on the inhibitor were calculated using the RESP fitting method, which is based on electrostatic potentials determined from quantum chemical computations. The GAFF2 parameters for the phosphate and pyrophosphate groups were allocated according to established techniques and evaluated against previously reported standards for phosphate-containing compounds.^40–42 The protonation states for the inhibitor were assigned using OpenBabel.¹² The complex is solvated in a cubic box of TIP3P water molecules⁴³ with a 15 Å buffer between the protein and the edges of the simulation box. Periodic boundary conditions are applied and the system is neutralized by adding 18 Na⁺ and 20 Cl⁻ ions to achieve a final concentration of 0.15 M. The use of water models like OPC has recently been reported to reproduce bulk solvent characteristics more accurately.⁴⁴ We have chosen to continue working with the TIP3P model as the system was found to be stable with a 1 fs time step without bond constraints. It also facilitated efficient sampling across prolonged simulation times.

We have performed several benchmarking MD simulations of length 1–5 µs to validate our choice of force field parameters for three systems, viz. protein, PUA and the PUA–RNase A complex solvated separately in TIP3P and OPC models of water. All three systems were found to remain stable in TIP3P water all along the simulated MD trajectories. However, despite following the same protocol for equilibration as in TIP3P water, unphysically large fluctuations in energy were observed while simulating the PUA–RNase A complex in OPC water at times larger than 1 ns. Needless to say, the electrostatic interactions play a pivotal role in stabilizing the complex and would require further benchmarking before we can adopt the OPC model to solvate it. We thus discuss only the results obtained using TIP3P water in the present study.

The solvated system for each model, comprised of 124 amino acid residues of RNase A and one PUA molecule in a box of ∼11 300 TIP3P water molecules, was then subjected to energy minimization, followed by gradual heating from 0 K to 300 K using a Langevin thermostat⁴⁵ with a collision frequency of 2 ps⁻¹. The system was further equilibrated under isothermal–isobaric (NPT) and isothermal–isochoric (NVT) conditions up to 200 ns. Finally, an NVT production run of 500 ns was carried out at 300 K. The trajectory was next extended up to 5 µs for both the closed conformers to track their location at and away from the active site. The particle mesh Ewald (PME)⁴⁶ method was employed for the long-range electrostatic interactions, whereas a cutoff of 10 Å was applied for short-range non-bonded interactions. The structural stability and flexibility of the simulated systems have been assessed through root mean square deviation (RMSD), root mean square fluctuation (RMSF) and radius of gyration (R_g) analyses of the equilibrated trajectory using the cpptraj⁴⁷ module of the AMBER24 software package.

2.4 Umbrella sampling and weighted histogram analysis method

The potential of mean force (PMF), projected along the end-to-end distance, d_ee of PUA in the PUA–RNase A complex, was calculated using the umbrella sampling technique.⁴⁸ To ensure sufficient configurational overlap between adjacent windows, a series of 18 overlapping umbrella windows were generated along d_ee. Each window was equilibrated and subsequently sampled for 5 ns at 300 K, employing a harmonic biasing potential of 10 kcal mol⁻¹ Å⁻² to restrain the system around the target distance. The resulting biased histograms were then unweighted and combined using the weighted histogram analysis method (WHAM),⁴⁹ applying a convergence tolerance of 10⁻⁴ to obtain the PMF profile. In the same manner, another PMF profile was generated where the major correlated active site residues were subjected to positional restraint with a biasing potential of 10 kcal mol⁻¹ Å⁻².

2.5 MM–GBSA free energy calculation

The binding free energy of a ligand–receptor complex in the presence of a solvent at a given temperature T is estimated as,

ΔG_bind,solv = ΔG_bind,vacuum + ΔG_solv,complex − (ΔG_solv,ligand + ΔG_{solv,receptor}) (12)

where ΔG_bind,vacuum is the binding free energy in the gas phase and the solvation free energy, ΔG_solv, consists of the solvation free energy due to the ligand (ΔG_solv,ligand), the receptor (ΔG_{solv,receptor}) and the complex (ΔG_solv,complex). ΔG_solv is expressed as the sum of its polar (ΔG_polar) and non-polar (ΔG_nonpolar) components.^50–52 ΔG_bind,vacuum is obtained by calculating the molecular mechanical energy, ΔE_MM and the entropic contributions. The overall molecular mechanical energy includes all the internal energies and interaction energies arising from electrostatic and van der Waals forces.^53,54

ΔE_MM = ΔE_int + ΔE_elec + ΔE_vdW (13)

The polar solvation free energy, ΔG_polar, was computed using the generalized Born (GB) model, which provides an efficient approximation to the linearized Poisson–Boltzmann equation.^55,56 The nonpolar solvation term, ΔG_nonpolar, was estimated from the solvent accessible surface area (SASA) of the solute.^57,58 All of the above terms were calculated as ensemble averages derived from the configurations obtained from the equilibrated classical MD trajectories.

3 Results and discussion

3.1 Geometry optimization

As mentioned earlier, the initial structure of the open conformation of PUA has been obtained from the crystal structure of the PUA–RNase A complex (PDB id: 1QHC).⁵ The initial closed conformation, on the other hand, was generated by twisting both the termini of the open conformer. These are presented in Fig. 2(a) and (b). The geometries of PUA for both conformations are first optimized in the gas phase at the HF/6-31G level and then re-optimized at the B3LYP/6-31G++(d,p) level. A significant reduction in total energy was observed due to the incorporation of both HF exchange and DFT correlation effects (Fig. S2 and Table S2). Further optimization with the larger def2-TZVP basis set led to a minimal change in the total energy within each model. The structures optimized with B3LYP/def2-TZVP have been used as input in the next steps of our work and presented in Fig. 2(c) and (d). Evidently, the closed and open conformations of PUA have been retained through the optimization process.

The total energy (in Ha) of the optimized structures are estimated to be −3995.2162 (closed) and −3995.1958 (open) (Table S2). In view of the energy difference of 53.56 kJ mol⁻¹, the open conformation appears to be more reactive and less kinetically stable than its closed counterpart. This observation is consistent with the HOMO–LUMO gap (ΔE_gap) values obtained from frontier molecular orbital (FMO) theory, as shown in Fig. 3.

Fig. 3 HOMO–LUMO gap for the (a), (c) and (e) closed and (b), (d) and (f) open conformers of PUA (a) and (b) in the gas phase, (c) and (d) water and (e) and (f) diethyl ether, including dispersion correction terms in the calculation of total energy, E.

3.2 Conformation-dependent global reactivity of PUA

The results presented in this section have been obtained starting with the optimized gas phase structures of PUA as shown in Fig. 2(c) and (d). The associated total energies of the HOMO and LUMO and the energy gap between them are highlighted in Fig. 3. The polarity of the environment is expected to be important in the stabilization and binding profile of PUA. Therefore, the results in the solvent phases have also been highlighted in this figure. Table 1 summarizes the computed values of changes in relative energy (ΔE_RE), energy gap (ΔE_gap), electronegativity (χ), chemical potential (μ), global hardness (η), global softness (S) and the electrophilicity index (ω) of the closed and open conformations of PUA in different phases. The implications of these results are discussed below.

Table 1 Calculated relative energy (RE), HOMO–LUMO energy gap (ΔE_gap), electronegativity (χ), chemical potential (μ), global hardness (η), softness (s), global softness (S) and the electrophilicity index (ω) for the gas phase, water and di-ethyl ether media. All the values except RE are given in eV and eV⁻¹ as applicable

PUA: closed	ΔERE (kJ mol^−1)	ΔEgap	χ	μ	η	S	s	ω
Gas	0	5.07	3.81	−3.81	2.53	0.20	0.39	2.86
Water	−135.21	5.15	3.95	−3.95	2.57	0.19	0.38	3.04
Diethyl ether	−81.13	5.14	3.87	−3.87	2.57	0.19	0.38	2.92

---

PUA: open	ΔERE (kJ mol^−1)	ΔEgap	χ	μ	η	S	s	ω
Gas	0	4.72	3.94	−3.94	2.36	0.21	0.42	3.29
Water	−141.51	4.89	3.95	−3.95	2.45	0.20	0.41	3.18
Diethyl ether	−87.95	4.82	3.89	−3.89	2.41	0.20	0.40	3.13

---

3.2.1 PUA in the gas phase. In the gas phase, the HOMO–LUMO gap is small for both conformations (compared to the solvated systems). This behavior is attributed to the strong electron-donating and accepting capabilities of the nucleobase components in PUA. The electron density of the HOMO is predominantly located on the pyrimidine base, whereas that of the LUMO is concentrated on the purine group (Fig. 3). The π-delocalization across the phenyl rings further stabilizes the planar structures of the purine and pyrimidine bases. In contrast, hydrogen bonding [–N–H⋯O] interaction between purine–N–H and pyrimidine–O (Fig. 2(c)) appears to stabilize the closed conformation. In the open conformer, this interaction is disrupted due to the wide separation of the aromatic rings, eliminating base-to-base interactions. The higher value of electrophilicity, ω (3.29 eV for open and 2.86 eV for closed) reflects a greater propensity of the open conformation to accept electronic charge, driven by its slightly more negative μ (−3.94 eV) and smaller η (2.36 eV).

3.2.2 PUA in water. When the system is transferred from the gas phase to a polar aqueous environment, both conformations show a significant increase in stability, as indicated by the substantial negative relative energy (ΔE_RE) values: −141.51 kJ mol⁻¹ for the open and −135.21 kJ mol⁻¹ for the closed conformers. The polar solvent efficiently screens charge fluctuations, resulting in a decreased polarizability and improved thermodynamic stability, as seen by an increase in energy gap (ΔE_gap) and chemical hardness (η) upon solvation. Remarkably, the open conformation has an equivalent electronegativity (χ = 3.95 eV) but a smaller HOMO–LUMO gap (ΔE_gap = 4.89 eV) than its gas-phase value, indicating that polar solvation promotes the reactivity of the extended structure with a lower number of intramolecular hydrogen bonds. This would also imply the formation of strong solute–solvent hydrogen-bonding interactions stabilizing the exposed polar functional groups of PUA that can compensate for the loss of intramolecular hydrogen bonding.

3.2.3 PUA in diethyl ether. As expected, a significantly lower stabilization of both conformations has been observed in a non-polar medium (diethyl ether) compared to water. As shown in Table 1, the ΔE_RE values are found to be −87.95 kJ mol⁻¹ for the open and −81.13 kJ mol⁻¹ for the closed conformations. The moderate rise in ΔE_gap shows a partial screening of charge delocalization by the non-polar environment in comparison to the gas phase. However, the global softness (S) and chemical potential (μ) vary for both conformations, suggesting that there is a little change in their overall reactivity. The open conformation has a reduction in its electrophilicity index (ω) compared to the gas phase. In contrast, a moderate increase in ω is observed for the closed conformer, suggesting that solvent stabilization improves the capacity of the closed conformer to take up electron density from a donor. A minimal variation in softness (s) between the gas phase and solvent phase for the closed conformer indicates the retention of reactivity of the conformer across variable polarity. The relative energy (RE) values vary across phases, as do other reactivity descriptors. Moreover, both the PUA conformations exhibit enhanced stabilization and structural flexibility upon solvation in both dielectric media.

3.3 Conformation-dependent local reactivity of PUA

Similar to Section 3.2, we have used the optimized gas phase structures of PUA, as shown in Fig. 2(c) and (d), to calculate the ground-state local reactivity of PUA conformers in terms of Fukui functions and intra-molecular non-covalent interactions.

3.3.1 Fukui functions. Table 2 summarizes the results of our estimation of the Fukui functions, as defined in Section 2.2.2, for each atom within the molecule. Specific atoms within PUA are found to exhibit high values of f⁺_k and f⁻_k, suggesting them as potential nucleophilic or electrophilic attack sites, respectively. To visualize these reactive regions, the iso-surfaces of Fukui functions f⁺_k and f⁻_k were generated using Multiwfn³² (Fig. 4(a)–(d)). In both the closed and open conformers, atomic regions in the nucleobases are found to be the most probable sites of attack, similar to the observation from the HOMO and LUMO orbitals shown in Fig. 3. The atoms C-2, C-10, and O-11, belonging to the nucleobase uracil, are more susceptible to nucleophilic attack in the closed conformer. On the other hand, opening of the dinucleotide exposes the N-6 and C-1 atoms in the uracil base widening the f⁺_k iso-surface. Atomic sites N-68, N-70, C-61, and C-62 in the adenine nucleobase are electron-deficient and rich centers as observed from the Fukui function values and are more accessible to nucleophiles and electrophiles, respectively, for the closed conformation. A negative f⁻_k has been observed at N-67 in the open conformer, indicating a low reactivity towards electrophiles. Both nucleobase regions in the closed conformer have been found to be highly susceptible to electron addition or removal, whereas for the open conformer, high susceptibility is displayed for the uracil base only.

Table 2 Fukui indices and local softness of most reactive sites in closed and open conformations in the gas phase

Atomic centres (PUA: closed)	f^+k	f^−k	s^+k	s^−k
O-11	0.047	0.026	0.01	0.006
C-10	0.044	0.012	0.008	0.002
C-2	0.075	0.035	0.015	0.007
C-61	0.041	0.054	0.008	0.01
N-68	0.028	0.063	0.006	0.012
N-70	0.021	0.097	0.004	0.019
C-62	0.006	0.046	0.001	0.009

---

Atomic centres (PUA: open)	f^+k	f^−k	s^+k	s^−k
C-2	0.08	0.1	0.016	0.02
N-67	0.029	−0.006	0.005	−0.001
N-6	0.017	0.15	0.003	0.03
C-1	0.05	0.2	0.01	0.04
O-11	0.07	0.11	0.014	0.22
C-10	0.055	0.037	0.011	0.037

---

Fig. 4 Iso-surfaces of the nucleophilic Fukui function f⁺_k (a) and (b) and the electrophilic Fukui function f⁻_k (c) and (d) for the (a) and (c) closed and (b) and (d) open conformations of PUA, respectively, in the gas phase. An iso-value of 0.005 was used for generating all surfaces. Green and blue regions denote areas of positive and negative Fukui function intensity, respectively, highlighting the spatial distribution of the most reactive atomic sites toward nucleophilic or electrophilic attack. NCI isosurfaces in the optimized (e) closed and (f) open conformations of PUA indicating stabilizing hydrogen-bonding interactions (blue), van der Waals contacts (green) and steric repulsion (red). (g) and (h) RDG scatter plots of closed and open PUA illustrating the relative strength and nature of non-covalent interactions through the distribution of sign(λ₂)ρ values.

3.3.2 Non-covalent interactions. The results of our calculations on both conformations, including the dispersion corrections, have been used to analyze the conformation-dependent intra-molecular non-covalent interaction (NCI) index in the gas phase using Multiwfn. The conformation dependence of three-dimensional NCI surfaces, shown in Fig. 4(e) and (f), has been generated using the RDG function and colored according to the sign of sign(λ₂)ρ as described in Fig. 4(g) and (h). The intra-molecular non-covalent interaction within the molecule is not that prominent, yet can be considered as an initial guide to the reactivity of the molecule. In the closed state, it has multiple predominant blue iso-surfaces (Fig. 4(e)), indicating the stabilization of the structure through intramolecular hydrogen bonds, accompanied by dispersive van der Waals interactions. In the open state, a strong hydrogen-bond interaction has been found in the 3′,5′-pyrophosphate group and a weaker interaction with the hydrogen donor from the 2′-OH group in the ribose ring to the P=O of 3′-phosphate of the adenosine half (Fig. 4(f)).

Several van der Waals interactions with minimal steric repulsions only within the rings have been observed in the molecule. The associated RDG scatter plot (Fig. 4(g) and (h)) shows strong peaks at the negative termini of the x-axis in the closed conformer. In contrast, the open conformation shows fewer hydrogen bonds and van der Waals interactions distributed across the backbone, while steric repulsion is extremely concentrated on the closed counterpart because of its locked conformation. These differences in non-covalent stabilization are consistent with the Fukui function results. In the closed conformer, electron accepting ability is localized near C-2 whereas in the open form, N-6 and C-1 coincide with altered non-covalent stabilization and enhanced nucleophilic character. These intramolecular interactions are expected to guide the pre-organization of the inhibitor, thereby affecting its conformational stability and flexibility within the protein binding pocket. This pre-organization may also indirectly influence interactions between PUA and RNase A by modifying the entropic and electrostatic factors during the complex formation.

3.4 Conformation dependent electrostatic potential

The molecular electrostatic potential (ESP) overlaid onto the van der Waals surface has been shown in Fig. 5(a) and (b) to highlight additional information about the conformational preferences of PUA. According to the distribution of ESP, the open conformation displays a narrower range, spanning −46.73 kcal mol⁻¹ to 81.64 kcal mol⁻¹ with a mean surface potential of 4.94 kcal mol⁻¹, and a diffused and spatially dispersed electrostatic profile. The closed conformation, on the other hand, shows highly localized pockets of electron density, with a broader range of ESP from −41.64 to 93.99 kcal mol⁻¹ and a mean of 4.82 kcal mol⁻¹, especially near the nucleobase moieties. Despite the fact that the closed conformation has both nucleobases prone to nucleophilic and electrophilic attacks obtained from its local electron densities (Section 2.2.2), their prominent potential localization produces distinct, well-defined interaction hot spots. The closed conformer will therefore be expected to sustain long-lasting non-covalent contacts with important residues due to these local propensities to bind more persistently with complementary sites in the active pocket. The open conformation, on the other hand, offers a more dynamic and less predictable interaction surface due to its fluctuating and less well-defined ESP properties. The open conformer is less advantageous for long-term binding because of these diffused potentials, which reduce the possibility of creating stable, directed connections.

Fig. 5 Molecular electrostatic potential (ESP) mapped onto the van der Waals surface for the (a) closed and (b) open conformations of PUA in the gas phase, highlighting the electron-rich sites with negative ESP (blue) and electron-deficient regions with positive ESP (red) relevant for electrophilic interactions.

Different features of ESP, as outlined above, qualitatively imply that the closed conformation may be a better candidate for sustained interaction with the active site, although the overall reactivity descriptors may still be comparable. This in turn suggests that conformational electrostatic characteristics alone may not be sufficient to determine binding affinity. Rather, the interaction between conformational variations of PUA and the dynamic behavior of the active site residues may be crucial in determining how strongly it binds and for how long.

The results presented above qualitatively justify the choice of d_ee as the sole discriminator of two conformations of PUA with different reactivity profiles associated with their optimized structures in the gas and solvent phases. To extend this study, we also observed analogous results when the –OH group orientation was changed at the C-2′ position of the furanose ring linked to the adenine base (an arabinose group instead of ribose) in the closed conformation of PUA. The overall energy remains unchanged, differing only by 3.6 kJ mol⁻¹, while preserving other reactivity parameters. Notably, the modified system exhibited a smaller HOMO–LUMO gap, accompanied by smaller electronegativity and electrophilicity indices (Table S3), implying that the inhibitor may have an enhanced preference towards electron donation and a reduced tendency for electron acceptance.

3.5 Effect of confinement on the reactivity of PUA at the active site of RNase A

In this section, we report yet another benchmarking study on the suitability of d_ee when the chosen conformations of PUA are embedded in the active site of RNase A. A kinetic pose of each conformer has been extracted from the classical MD trajectory at the end of the NVT equilibration step of 200 ns (Section 2.3) and subjected to single point energy calculation (B3LYP/def2-TZVP) followed by analysis of the associated reactivity parameters. The kinetic poses of PUA thus extracted are found to have d_ee = 2.92 Å (closed) and 15.2 Å (open). A significant kinetic rearrangement induced by confinement to the active site is evident from Fig. 6(a) and (b) obtained by superimposing these kinetic poses to the related optimized geometries shown in Fig. 2. We have also extracted the coordinates and point charges on the atoms belonging to active site residues surrounding each kinetic pose of PUA as shown in Fig. 6(c) and (d). Estimations of the single point energy and reactivity parameters of PUA are then repeated by including these charges at their extracted locations.

Fig. 6 Superposition of the optimized geometries (green) and the corresponding kinetic pose extracted from the equilibrated classical MD trajectory in the RNase A active site (pink) for the closed (a) and open (b) conformations of PUA. Closed (c) and open (d) conformations surrounded by the point charges of active site residues. Blue and red-colored vdW representations show the positive and negative partial charges, respectively.

The global reactivity parameters, obtained for the kinetic poses with and without the active site charges are presented in Table 3. Both conformations are found to demonstrate similar global softness, with the open form exhibiting lower softness (S = 0.20 eV⁻¹) compared to the closed conformer (S = 0.23 eV⁻¹), indicative of a more flexible and extended structure of the former. Interestingly, when the electrostatic field of the active site is explicitly considered, the closed conformation exhibits a significant increase in softness (S = 0.28 eV⁻¹), while that of the open conformation remains mostly unchanged (S = 0.20 eV⁻¹). Table 4(A) and (B) present the values of the Fukui functions for the highly reactive atomic regions in PUA. A significant reactivity concentrated at C-2 and O-11 in the uracil moiety and N-68,70 in the adenine moiety are observed from local reactivity analysis of the closed conformer in the first approach. The open conformation exhibits enhanced electrophilic properties as shown in Fig. 7 at several locations, notably N-68 and N-70, consistent with the reactivities obtained for optimized PUA in the gas phase (Fig. 4).

Table 3 Total ground state energy (E_total), HOMO–LUMO energy gap (ΔE_gap) and other global reactivity descriptors for kinetic poses of the closed and open conformations of PUA with and without being surrounded by active site charges. All the values except RE are given in eV and eV⁻¹ as applicable

PUA, without active site charges	Etotal (Hartree)	ΔEgap	χ	μ	η	S	s	ω
Closed	−3994.9535	4.16	3.32	−3.72	2.08	0.23	0.47	3.325
Open	−3994.8517	4.85	4.01	−4.01	2.42	0.20	0.41	3.323

---

PUA, with active site charges	Etotal (Hartree)	ΔEgap	χ	μ	η	S	s	ω
Closed	−3998.8207	3.53	5.35	−5.35	1.76	0.28	0.56	8.11
Open	−3997.9402	4.78	6.13	−6.13	2.39	0.20	0.41	7.85

---

Table 4 Fukui indices for the kinetic poses of the closed and open conformers of PUA with and without embedded active site electrostatics

(A) PUA, without active site charges
	Closed		Open
	f^+k	f^−k	f^+k	f^−k
C-1	0.044	0.006	0.046	0.055
C-2	0.067	0.016	0.068	0.026
O-11	0.051	0.007	0.061	0.054
C-61	0.016	0.076	0.038	0.055
C-62	0.0009	0.054	0.007	0.034
N-67	0.007	0.045	0.028	0.042
N-68	0.012	0.098	0.008	0.065
N-70	0.013	0.136	0.017	0.076

---

(B) PUA, with active site charges
	Closed		Open
	f^+k	f^−k	f^+k	f^−k
C-1	0.028	0.0002	0.036	0.067
C-2	0.027	0.013	0.05	0.031
O-11	0.0028	0.0033	0.046	0.071
C-61	0.007	0.08	0.084	0.024
C-62	0.0004	0.057	0.015	0.012
N-67	0.0007	0.046	0.048	0.028
N-68	0.002	0.106	0.052	0.018
N-70	0.006	0.144	0.043	0.027

---

Fig. 7 Iso-surfaces of the nucleophilic Fukui function f⁺_k (a) and (b) and the electrophilic Fukui function f⁻_k (c) and (d) for the (a) and (c) closed and (b) and (d) open conformations of PUA, respectively, for the kinetic pose of PUA at the active site of RNase A. An iso-value of 0.005 was used for generating all surfaces. Green and blue regions denote areas of positive and negative Fukui function intensity, respectively, highlighting the spatial distribution of the most reactive atomic sites toward nucleophilic or electrophilic attack. NCI isosurfaces in the optimized (e) closed and (f) open conformations of PUA indicating stabilizing hydrogen-bonding interactions (blue), van der Waals contacts (green) and steric repulsion (red). (g) and (h) RDG scatter plots illustrating the relative strength and nature of non-covalent interactions through the distribution of sign(λ₂)ρ values.

The incorporation of active site electrostatic charges results in the closed conformation maintaining its reactive characteristics while displaying a refinement of the values of the Fukui functions (Fig. 8). This indicates how protein electrostatics may tune the reactivity at particular atoms of the ligand. Electrophilic contributions at nitrogen atoms of the adenine moiety are found to be considerable, reflecting their role in hydrogen bonding with RNase A residues. In the open conformation, the charges in the environment prompt a redistribution of Fukui indices towards the uracil moiety, resulting in elevated f⁻_k values at C-1, C-2 and O-11, while the electrophilic character of the adenine nitrogen atoms is comparatively diminished. In both the closed and the open conformations with and without the active site charges, a reduction of stabilizing hydrogen-bonds and dispersive interactions is observed compared to that encountered for the optimized PUA in the gas phase. Thus, the electrostatic charges at the active site of RNase A appear to polarize the closed conformation preferentially, increasing its sensitivity to electronic perturbations and strengthening its bound state character.

Fig. 8 Iso-surfaces of the nucleophilic Fukui function f⁺_k (a) and (b) and the electrophilic Fukui function f⁻_k (c) and (d) for the (a) and (c) closed and (b) and (d) open conformations of PUA, respectively, for the kinetic pose embedded with active site charges. An iso-value of 0.005 was used for generating all surfaces. Green and blue regions denote areas of positive and negative Fukui function intensity, respectively, highlighting the spatial distribution of the most reactive atomic sites toward nucleophilic or electrophilic attack. NCI isosurfaces in the optimized (e) closed and (f) open conformations of PUA indicate stabilizing hydrogen-bonding interactions (blue), van der Waals contacts (green) and steric repulsion (red). (g) and (h) RDG scatter plots illustrating the relative strength and nature of non-covalent interactions through the distribution of sign(λ₂)ρ values.

3.6 Results from classical MD simulation

Two classical MD trajectories (of length 5000 ns each), starting from two different conformations of PUA docked to RNase A, have been analyzed by monitoring the root mean square deviation (RMSD), radius of gyration (R_g), and root mean square fluctuations (RMSF) of the enzyme–ligand assembly in water over time (Fig. S3). Both simulated systems are found to have stable, compact structures with the thermal fluctuations attributed primarily to dynamic turns and coils on the surface of the enzyme.

Fig. 9(a) and (b) present the most stable poses of PUA docked to RNase A with ΔG_binding (kcal mol⁻¹) = −8.872 (closed) and −9.783 (open) estimated from AutoDock Vina.³⁴ Although these values do not represent rigorous binding free energies, their magnitudes are consistent with the experimentally reported ΔG_binding of approximately −10.4 kcal mol⁻¹ for the PUA–RNase A complex. Upon docking, the closed conformation remained closed with d_ee = 3.03 Å, registering a slight increase compared to d_ee = 2.96 Å in the optimized structure in the gas phase (Fig. 2(c)). The docked pose of the open conformer is found to have d_ee = 17.68 Å, which is increased from the corresponding value of 13.59 Å observed in the crystal structure of the PUA–RNase A complex.⁵

Fig. 9 Binding of PUA (licorice) to the active site (yellow, van der Waals surface) of RNase A (grey) showing the (a) and (b) initial docked pose of the (a) closed and (b) open conformations of PUA. Also shown in (b) is the alignment of the crystallographic pose of PUA (PDB id: 1QHC⁵) (purple, CPK). (c) and (d) Time evolution of the end-to-end distance, d_ee along the classical MD trajectories for the (c) closed and (d) open conformations of PUA. The open conformation unbinds from the active site at 300 ns. (e) and (f) Representative snapshots extracted from the classical MD trajectories at 350 ns showing (e) the bound state of the closed conformation of PUA at the active site of RNase A and (f) the un-bound, open conformation of PUA in water near the active site cavity.

Along the MD trajectory, two important distances have been monitored to understand the dynamical fluctuations of the initial binding poses described above. First, several active site residues are identified with atoms located within 3 Å from any atom of PUA in the bound state. These residues include (i) Val-43, Thr-45, Asn-67, Arg-85,  His-119, Phe-120 and Asp-121 (closed) and (ii) His-12, Gln-11, Lys-41, Asn-67, His-119, Phe-120 and Asp-121 (open). We have next estimated the distance, d between the center of mass (COM) of PUA and that of the active site residues as listed above in each case. In the initial docked poses, as shown in Fig. 9(a) and (b), the distance, d = 4.92 Å (closed) and 3.70Å (open) characterizes the respective bound states. When sampled along the MD trajectory until 5 µs, d remained unchanged for the closed conformation. The open conformation, on the other hand, displayed an increased value of 12.58 Å at 300 ns followed by a discontinuous jump to d = 19.29 Å at 350 ns. A monotonic increase in the value of d was observed thereafter. Therefore, we conclude that the open conformation of PUA underwent the first unbinding event at t = τ_f = 300 ns and a complete unbinding at t = τ_u = 350 ns.

We next analyzed the variations of d_ee for each conformation of PUA up to 500 ns as shown in Fig. 9(c) and (d). When sampled along the MD trajectory until 500 ns, the closed conformation, in its bound form, was found to fluctuate with transient increments in d_ee up to 20 Å. However, it remained predominantly (63.85%) closed until 500 ns. A representative structure, extracted from the MD trajectory at 350 ns and shown in Fig. 9(e), is found to be associated with d_ee = 8.19 Å characteristic of a closed conformation. Evidently, despite the dynamical fluctuations, the closed conformation of PUA remains strongly anchored within the binding cleft. The atoms of PUA are found to spend 97.03% of their residence time at the active site within 4 Å of hydrophilic residues retaining critical hydrogen bonds with minimal dynamic interaction with the hydrophobic residues.

As shown in Fig. 9(d), the open conformation of PUA is retained in 86.39% of the structures sampled along the associated MD trajectory with infrequent transitions to d_ee < 10 Å. During its short residence until 300 ns at the active site, the open conformation of PUA was found to spend only 73.85% within 4 Å of hydrophilic residues. In contrast to its closed counterpart, the open conformation was found to establish hydrophobic contacts with the active site residues in 25.05% of the sampled structures. After 300 ns, the open conformer of PUA gradually drifts away at the active site of RNase A. Fig. 9(f) records its location at 350 ns in bulk water near the mouth of the active site cavity. Evidently, the release from the active site not only increases the entropy of the open conformation, but also augments stabilizing interactions with the polar medium.

We have further examined the simulated systems by extending each of the MD trajectories to 5 µs. The closed conformation was found to remain bound to the active site cleft all throughout the simulation. In contrast, the open conformation did not return to the active site cavity within 5 µs. Instead, it explored the surface of the enzyme before eventually binding to a site near the N-terminal of RNase A. These results re-affirm our choice of PUA as an ideal model of two distinct conformations with remarkably different binding profiles to RNase A related to different dynamical interactions with the hydrophobic and hydrophilic residues at the active site of RNase A.

The conformational dependence of dynamical interaction of PUA with the amino acid residues of RNase A has been presented in Fig. 10. A comprehensive list of residues interacting with PUA in both conformers has been enumerated and the associated elements of the mass weighted dynamical cross correlation matrix have been plotted against each of the listed residues. Evidently, near the N-terminal of RNase A, fluctuations of both conformers are correlated with the motion of Lys-7 and Gln-11 and the catalytically important His-12. The open conformer is found to be more correlated with the residues at the N-terminus, while the closed conformer exhibits significant correlated fluctuations with the C-terminus as well as the N-terminus of the enzyme. The highest dynamical stabilization is exerted on the closed conformation by the residues Val-43, Lys-66, Arg-85 and Lys-98. The fluctuations of the bound open conformation are found to be anti-correlated with the motion of Arg-85 and Lys-98 and relatively weakly correlated with Val-43 and Lys-66. These results emphasize the conformational dependence of dynamical coupling of the fluctuations of both the ligand and enzyme.

Fig. 10 Mass-weighted dynamical cross-correlations of PUA with amino acid residues of RNase A for the (a) closed and (b) open conformers.

In Fig. 11, we have presented the potential-of-mean-force (PMF) as a function of d_ee that highlights the transition of the open conformation (d_ee ∼ 20 Å) to the closed conformation (d_ee ∼ 5 Å). It clearly shows that within the active site of RNase A, the transition from open to closed conformation of PUA is mostly downhill in free energy. A small barrier around d_ee ∼ 13 Å is encountered during this transition. We also present in Fig. 11 the PMF as derived from umbrella sampling where the open conformation has been restricted all throughout in strong non-covalent contact with selected, highly correlated amino acid residues such as Ala-4, Lys-7, Gln-69, Thr-70, Val-118, and His-119. Interestingly, the free energy barrier still persists along the PMF profile.

Fig. 11 PMF profiles projected along the end-to-end distance, d_ee with (red) and without (black) restraints, highlighting the higher barrier in the unrestrained simulation between 10 and 15 Å. Representative interaction maps from these regions are shown, reflecting several solvent interactions and hydrogen-bond networks. Hydrophobic contacts as well as the residues with restraint are shown in magenta, residues having hydrogen-bond interactions are depicted as sticks and hydrogen-bonded O atoms of water molecules are represented as cyan spheres. The dotted lines connect the donor and acceptor atoms participating in hydrogen bonds.

Five more independent initial structures of the open conformer of PUA have been chosen, one with d_ee = 13.59 Å from the crystal structure (Fig. 1(b))⁵ and four others with d_ee = 15.19–19.98 Å, extracted from the equilibrated MD trajectory of the PUA–RNase A complex. Each of these systems was then used as the input to a classical MD simulation until 500 ns. The open conformer dissociates from the active site of RNase A along all these trajectories. The times taken for the first event of unbinding (τ_f) and complete unbinding (τ_u) from the active site are highlighted in Table 5. It is therefore concluded that the rapid unbinding of the open conformation of PUA is independent of the initial structure.

Table 5 The time taken for the initial unbinding event (τ_f) and that for complete unbinding (τ_u) of the open conformation of PUA from the active site of RNase A as observed along five independent classical MD trajectories starting from different end-to-end distances, d_ee (t = 0) (Å)

	Initial dee (Å)	Time of first unbinding event (ns)	Time of complete unbinding (ns)
1	13.59	160.3	200
2	15.19	24	33.8
3	16.74	47	125
4	18.40	44.8	50.6
5	19.98	41.9	75

---

A detailed comparison of the restrained and unrestrained MD simulations of length 600 and 300 ns, respectively, demonstrated that Gln-69, a flexible loop residue, plays a critical role in regulating the stability of PUA within the RNase A active site. In the restrained simulation, Gln-69 is limited to a single orientation, whereas in the unrestrained trajectory, it exhibits two different orientations with respect to the active site that can be defined in terms of the sidechain dihedral angles χ₁ ≡ N–C_α–C_β–C_γ and χ₂ ≡ C_α–C_β–C_γ–C_δ of Gln-69. In the in orientation, the Gln-69 sidechain is directed towards the active site with χ₁ < 0° and −50° ≤ χ₂ ≤ −180°. The out orientation, on the other hand, faces the solvent with χ₁ > 0° and −50° ≤ χ₂ ≤ 180°. The fluctuations of the loop on which Gln-69 is located results in a coupled movement of the open conformation towards the solvent. Within the limited sampling of 600 ns of the restricted MD trajectory, Gln-69 is found to establish stable hydrogen bonds with PUA in 15.49% and 11.55% of the structures when Gln-69 is in its in and out orientations, respectively. This observation is further supplemented by the high degree of dynamical correlation between Gln-69 and PUA as shown in Fig. 10(b). It will be interesting to perform further enhanced sampling simulation to probe how the loop dynamics through Gln-69 can facilitate the displacement of the open conformation of PUA towards the solvent by inducing an important loss of stabilizing H-bonding interactions.

3.7 Correlation with binding affinity

All the results presented above emphasize how the conformational diversity of PUA is intricately correlated with the reactivity, interactions with neighboring residues and residence time at the active site of RNase A. No attempt has been made in this work to obtain accurate estimates of the free energy of (un)binding of PUA with RNase A. This is mainly because any such estimate will be reliable in the limit of a very long time of simulation, preferably sampling all possible binding modes of the ligand at the active site and on the surface of the enzyme. Instead, as depicted in Fig. 12, we present here evaluation of the MM–GBSA binding free energy, ΔG_bind,solv for both the closed and open conformations of PUA for a series of MD trajectory segments of different lengths extracted from the 5 µs long trajectory. For the closed conformation, the persistence of the observed, large negative values of ΔG_bind,solv, as shown in Fig. 12(a), correlates well with its stable bound state with minor fluctuations arising due to local conformational arrangements of PUA or the active site residues within the binding pocket. In contrast, the open conformation, despite having an initial negative binding free energy, becomes increasingly less favorable with the loss of PUA–RNase A interactions with time even when it is present within the binding pocket. This is highlighted in Fig. 12(b). After it exits the active site, as shown in Fig. 12(c), a non-monotonic stabilization pattern is reflected in the values of ΔG_bind,solv, consistent with transient interactions with different regions of the enzyme surface before its localization in water near the mouth of the active site cavity. However, as mentioned above, these estimates are indicative only. A more rigorous computation of ΔG_bind,solv will be necessary for quantitative mapping onto the experimental binding affinity.

Fig. 12 Estimate of the MM–GBSA binding free energy (ΔG_bind,solv) of PUA for (a) the closed conformation assessed at different time intervals up to 5 µs, (b) the open conformation during the initial 300 ns before the first unbinding event and (c) the open conformation during a prolonged 5 µs simulation.

4 Conclusion

In this article, when bound to the active site of RNase A, PUA, a highly potent inhibitor of RNase A, can adopt at least one more stable conformation in addition to that observed in the crystal structure. We selected two conformations of PUA and then systematically examined their structural, electronic, and interaction profiles under different conditions. The closed conformer is found to be stabilized by a strong intramolecular hydrogen bond and π-delocalization across the nucleobases, according to gas-phase electronic structure analysis. In contrast, the open conformer has increased softness and electrophilicity, because it does not have these intramolecular contacts. Although the open form becomes slightly more preferred as a result of the solvation of exposed polar groups and the breaking of intramolecular connections, solvent-phase calculations showed that both conformers were significantly stabilized upon hydration. The variation of global reactivity descriptors across solvents showed that the two differences decrease in polar media for both, indicating that environmental polarity influences their inherent reactivity in a similar fashion. Local reactivity analysis using Fukui functions indicates the closed conformer has well-defined, localized sites for both nucleophilic and electrophilic attack on each nucleobase. In contrast, the flexible and solvent-exposed geometry of the open form causes these reactive centers to shift toward the uracil moiety. These findings were supported by NCI and RDG analyses, which showed that the closed form had numerous stabilizing hydrogen-bonded and dispersion contacts while the open form had fewer, more diffuse non-covalent interactions. This result was supported by the electrostatic potential maps, which show that the open form has a more diffuse, dynamically fluctuating electrostatic surface that is less favorable for long-lived contacts, while the closed structure exhibits sharply localized electrostatic pockets that can interact with the enzyme active site over an extended period of time. These results justify our ad hoc choice of two specific conformations of PUA without undertaking a computationally expensive exploration of its high-dimensional conformational space. The results from classical MD simulations emphasize the dynamical variance of these two conformations. Most importantly, it is predicted that the closed conformation may indeed be the one responsible for the high inhibitor potency of PUA.

The latest advancements in computer simulation studies on the dynamics of ligand–protein binding/unbinding are capable of putting forward deeper insight into how the fluctuations and association/dissociation of a ligand–protein complex alter the molecular mechanism of inhibition by the ligand. It may be noted that the results presented in this article are expected to serve as the much-needed starting point for further investigation of the role of conformational diversity in the mechanism of inhibition by PUA. Work is in progress in this direction and will be reported in due course.

Author contributions

SP and AD performed the quantum chemical calculations. SP carried out the classical MD simulations. AD and ST conceptualized the article. SP, AD and ST wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5cp04621a.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements

The research presented in this article is funded by the Science and Engineering Research Board (SERB), India (grant number CRG/2018/002005 and CRG/2022/007493). Additional computational support received from the HPC facility created under the DST-FIST scheme (SR/FST/CSII-011/2005 and SR/FST/CSII-026/2013) and that under the National Supercomputing Mission (NSM), Government of India, supported by the Centre for Development of Advanced Computing (CDAC), Pune, is gratefully acknowledged.

Notes and references

1. E. Procko, R. Hedman, K. Hamilton, J. Seetharaman, S. J. Fleishman, M. Su, J. Aramini, G. Kornhaber, J. F. Hunt, L. Tong, G. T. Montelione and D. Baker, J. Mol. Biol., 2013, 425, 3563–3575.
2. C. A. Evans, T. Liu, A. Lescarbeau, S. J. Nair, L. Grenier, J. A. Pradeilles, Q. Glenadel, T. Tibbitts, A. M. Rowley, J. P. DiNitto, E. E. Brophy, E. L. O’Hearn, J. A. Ali, D. G. Winkler, S. I. Goldstein, P. O’Hearn, C. M. Martin, J. G. Hoyt, J. R. Soglia, C. Cheung, M. M. Pink, J. L. Proctor, V. J. Palombella, M. R. Tremblay and A. C. Castro, ACS Med. Chem. Lett., 2016, 7, 862–867.
3. L. El Khoury, Z. Jing, A. Cuzzolin, A. Deplano, D. Loco, B. Sattarov, F. Hédin, S. Wendeborn, C. Ho, D. El Ahdab, T. Jaffrelot Inizan, M. Sturlese, A. Sosic, M. Volpiana, A. Lugato, M. Barone, B. Gatto, M. L. Macchia, M. Bellanda, R. Battistutta, C. Salata, I. Kondratov, R. Iminov, A. Khairulin, Y. Mykhalonok, A. Pochepko, V. Chashka-Ratushnyi, I. Kos, S. Moro, M. Montes, P. Ren, J. W. Ponder, L. Lagardère, J.-P. Piquemal and D. Sabbadin, Chem. Sci., 2022, 13, 3674–3687.
4. N. Russo and R. Shapiro, J. Biol. Chem., 1999, 274, 14902–14908.
5. D. D. Leonidas, R. Shapiro, L. I. Irons, N. Russo and K. R. Acharya, Biochemistry, 1999, 38, 10287–10297.
6. N. Doucet, T. B. Jayasundera, M. Simonović and J. P. Loria, Proteins: Struct., Funct., Bioinf., 2010, 78, 2459–2468.
7. E. L. Kovrigin, R. Cole and J. P. Loria, Biochemistry, 2003, 42, 5279–5291.
8. G. N. Hatzopoulos, D. D. Leonidas, R. Kardakaris, J. Kobe and N. G. Oikonomakos, FEBS J., 2005, 272, 3988–4001.
9. Y.-C. Lee, P. L. Jackson, M. J. Jablonsky and D. D. Muccio, Arch. Biochem. Biophys., 2007, 463, 37–46.
10. D. D. Leonidas, R. Shapiro, L. I. Irons, N. Russo and K. R. Acharya, Biochemistry, 1997, 36, 5578–5588.
11. S. Polydoridis, D. D. Leonidas, N. G. Oikonomakos and G. Archontis, Biophys. J., 2007, 92, 1659–1672.
12. N. M. O'Boyle, M. Banck, C. A. James, C. Morley, T. Vandermeersch and G. R. Hutchison, J. Cheminf., 2011, 3, 33.
13. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. M. Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16 Revision C.01, Gaussian, Inc., Wallingford, CT, 2016.
14. P. Geerlings, F. De Proft and W. Langenaeker, Chem. Rev., 2003, 103, 1793–1874.
15. R. G. Parr, Density functional theory of atoms and molecules, 1989, pp. 5–15.
16. S. K. Ghosh and P. K. Chattaraj, Concepts and methods in modern theoretical chemistry: electronic structure and reactivity, CRC Press, 2016.
17. P. K. Chattaraj, Chemical reactivity theory: a density functional view, CRC Press, 2009.
18. Y. Li and J. N. Evans, J. Am. Chem. Soc., 1995, 117, 7756–7759.
19. R. G. Parr and W. Yang, J. Am. Chem. Soc., 1984, 106, 4049–4050.
20. P. K. Chattaraj and D. Chakraborty, Chemical reactivity in confined systems: theory, modelling and applications, John Wiley & Sons, 2021.
21. C. Morell, A. Grand and A. Toro-Labbé, J. Phys. Chem. A, 2005, 109, 205–212.
22. U. Sarkar, J. Padmanabhan, R. Parthasarathi, V. Subramanian and P. Chattaraj, J. Mol. Struct. THEOCHEM, 2006, 758, 119–125.
23. L. Pauling, J. Am. Chem. Soc., 1932, 54, 3570–3582.
24. R. G. Parr and R. G. Pearson, J. Am. Chem. Soc., 1983, 105, 7512–7516.
25. R. Vijayaraj, V. Subramanian and P. Chattaraj, J. Chem. Theory Comput., 2009, 5, 2744–2753.
26. R. Pal and P. K. Chattaraj, J. Comput. Chem., 2023, 44, 278–297.
27. H. Chermette, P. Boulet and S. Portmann, Reviews of Modern Quantum Chemistry: A Celebration of the Contributions of Robert G Parr (In 2 Volumes), World Scientific, 2002, pp. 992–1012.
28. P. K. Chattaraj and D. Chakraborty, Electron Density: Concepts, Computation and DFT Applications, John Wiley & Sons, 2024.
29. W. Yang and W. J. Mortier, J. Am. Chem. Soc., 1986, 108, 5708–5711.
30. W. Langenaeker, F. De Proft and P. Geerlings, J. Phys. Chem., 1995, 99, 6424–6431.
31. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-Garca, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.
32. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
33. T. Lu, J. Chem. Phys., 2024, 161, 082503.
34. O. Trott and A. Olson, Effic. Optim. Multithreading, 2009, 31, 455–461.
35. J. Eberhardt, D. Santos-Martins, A. F. Tillack and S. Forli, J. Chem. Inf. Model., 2021, 61, 3891–3898.
36. R. Anandakrishnan, B. Aguilar and A. V. Onufriev, Nucleic Acids Res., 2012, 40, W537–W541.
37. R. T. Raines, Chem. Rev., 1998, 98, 1045–1066.
38. D. Case, H. Aktulga, K. Belfon, I. Ben-Shalom, J. Berryman, S. Brozell, D. Cerutti, T. Cheatham, G. Cisneros and V. Cruzeiro, et al., J. Am. Chem. Soc., 2020, 142, 3823–3835.
39. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
40. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
41. B. Koca Fndk, M. Jafari, L. F. Song, Z. Li, V. Aviyente and K. M. Merz Jr, J. Chem. Theory Comput., 2024, 20, 4298–4307.
42. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
43. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
44. S. Bosio, D. Gazzoni, C. Domene, M. Masetti and S. Furini, J. Chem. Theory Comput., 2026, 22, 1177–1186.
45. R. J. Loncharich, B. R. Brooks and R. W. Pastor, Biopolymers, 1992, 32, 523–535.
46. T. Darden, D. York and L. Pedersen, et al., J. Chem. Phys., 1993, 98, 10089.
47. D. R. Roe and T. E. Cheatham III, J. Chem. Theory Comput., 2013, 9, 3084–3095.
48. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
49. P. Tao, A. J. Sodt, Y. Shao, G. Konig and B. R. Brooks, J. Chem. Theory Comput., 2014, 10, 4198–4207.
50. A. Onufriev, D. Bashford and D. A. Case, Proteins: Struct., Funct., Bioinf., 2004, 55, 383–394.
51. H. Luo and K. Sharp, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 10399–10404.
52. P. A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. Huo, L. Chong, M. Lee, T. Lee, Y. Duan and W. Wang, et al., Acc. Chem. Res., 2000, 33, 889–897.
53. V. Tsui and D. A. Case, Biopolymers, 2000, 56, 275–291.
54. A. Onufriev, D. Bashford and D. A. Case, J. Phys. Chem. B, 2000, 104, 3712–3720.
55. W. C. Still, A. Tempczyk, R. C. Hawley and T. Hendrickson, J. Am. Chem. Soc., 1990, 112, 6127–6129.
56. J. Srinivasan, M. W. Trevathan, P. Beroza and D. A. Case, Theor. Chem. Acc., 1999, 101, 426–434.
57. L. R. Pratt and D. Chandler, J. Chem. Phys., 1980, 73, 3434–3441.
58. J. Dzubiella, J. Swanson and J. McCammon, J. Chem. Phys., 2006, 124, 084905.

---