Multi-Site Targeting Potential of 3CLpro Inhibitors: Unveiling Diverse Mechanisms Through Blind Docking and Brownian Dynamics Simulations

Blind docking of 3CLpro inhibitors to 3CLpro monomer

To explore the full range of potential binding interactions to the 3CLpro monomer without prior bias, blind docking simulations^[42,43] were performed for each of the 73 identified inhibitors against the 3CLpro monomer using the AutoDock 4.1 suite.^[49] This approach allows for an unbiased mapping of possible inhibitor binding poses ^[43] independent of pre-defined binding site information. AutoDock 4.1 was chosen for this study due to its demonstrated ability to predict correct binding poses based on binding energy, even in the absence of prior knowledge about the binding site.^[42,43]

For each docking simulation, the 3CLpro monomer was treated as a rigid body held stationary in the simulation space. Conversely, each ligand (inhibitor) molecule was initiated from a random position and orientation. To thoroughly explore the conformational space, 150 independent docking runs were performed for each ligand using the hybrid Lamarckian Genetic Algorithm with Local Search (GA-LS) method.^[50] The parameters for these runs included a population size of 200, a maximum of 3,000,000 energy evaluations, a maximum of 27,000 generations, and 300 cycles of local search refinement. Energy grids encompassing the 3CLpro monomer were designed to extend beyond the protein’s spatial dimensions by 40% in each Cartesian direction, ensuring that a sufficient area was available for ligand exploration. The grid potentials were calculated with a spatial resolution of 0.5 Å.

The 150 docking simulations performed for each of the 73 inhibitors produced a total of 10,950 unique binding poses on the 3CLpro monomer surface. To characterize the local environment of these poses, we identified proximal 3CLpro surface residues based on a 7.0 Å distance cutoff between the residue’s alpha-carbon (Cα) atom and any non-hydrogen atom of the inhibitor. Subsequent hierarchical clustering of these 10,950 poses using Ward’s algorithm in R^[51] with a 2.0 Å distance threshold revealed seven distinct clusters of ligand poses whose centers are separated by at least 10.0 Å. The binding site for each cluster was then defined by the three proximal surface residues exhibiting the highest frequency of interaction based on the aforementioned distance threshold.

Modeling the diffusional association of 3CLpro monomers

The effect of inhibitors bound at the different identified sites on the diffusional association of 3CLpro monomers was investigated through Brownian dynamics (BD) simulations using the SDA software package.^[52,53] We compared the association of two unbound monomers (M+M → M:M) with the associations of an unbound monomer (M) and a monomer with an inhibitor bound at one of the seven identified binding sites, ML(n): n=1,2,3…7, as shown in the following reaction schemes:

$\text{M}+\text{M}\to \text{M}:\text{M}$

$\text{M}+\text{ML}\left(\text{1}\right)\to \text{M}:\text{ML}\left(\text{1}\right)$

$\text{M}+\text{ML}\left(\text{2}\right)\to \text{M}:\text{ML}\left(\text{2}\right)$

$\text{M}+\text{ML}\left(\text{3}\right)\to \text{M}:\text{ML}\left(\text{3}\right)$

$\text{M}+\text{ML}\left(\text{4}\right)\to \text{M}:\text{ML}\left(\text{4}\right)$

$\text{M}+\text{ML}\left(\text{5}\right)\to \text{M}:\text{ML}\left(\text{5}\right)$

$\text{M}+\text{ML}\left(\text{6}\right)\to \text{M}:\text{ML}\left(\text{6}\right)$

$\text{M}+\text{ML}\left(\text{7}\right)\to \text{M}:\text{ML}\left(\text{7}\right)$

The dynamic process of association for each reaction was simulated by propagating 50,000 Brownian dynamics (BD) trajectories for each of the reactions above. These trajectories were generated by solving the translational and rotational diffusion equations over time, employing the Ermak-McCammon^[54] algorithm as implemented in the SDA package (version 7.2). The resulting translational and rotational displacements and Δw are given by:

$\Delta r\left(t\right)=\frac{D\Delta t}{kBT}F\left(t\right)+R\left(t\right)$

$\Delta w\left(t\right)=\frac{D_r\Delta t}{kBT}T\left(t\right)+W\left(t\right),$

In these equations, D and D_r​ represent the relative translational and rotational diffusion coefficients, respectively, and Δt denotes the simulation time step. The terms F and T describe the potential force and torque acting on the interacting reactants, the terms R and W are the random translational and angular displacements due to reactant collisions with the solvent, and K_B and T denote the Boltzmann constant and absolute temperature, respectively. The translational and rotational diffusion coefficients (D and D_r​) were calculated using the HYDROPRO software.^[55]

For each BD simulation trajectory, the centers of mass (COMs) of the two reactants were initially separated by 200.0 Å, with the 3CLpro monomer (M) placed at the origin. This large starting distance was selected because the electrostatic potential around each monomer becomes nearly isotropic beyond 80 Å, and with an average radius of 35 Å, it guaranteed a minimum initial surface separation of 130 Å. The simulation time step was dynamically adjusted: 1.0 ps for COM-COM distances less than 130 Å, increasing linearly by 0.5 ps Å^-1 at larger separations. The trajectory was terminated once the COM-COM separation exceeded 300 Å.

The forces governing the interaction between reactants were derived from steric, desolvation, and electrostatic contributions. Steric exclusion was implicitly enforced by preventing overlap using a 1 Å spaced exclusion grid centered on each reactant. Electrostatic forces were calculated by determining the electrostatic potential generated by one monomer at each atom of the other and multiplying by the atom’s charge. This electrostatic potential was obtained by numerically solving the nonlinear Poisson-Boltzmann equation using the APBS program^[56] on a 201 Å × 201 Å × 201 Å grid with 1 Å spacing, centered on each reactant. The solvent and reactant interior relative dielectric constants were set to 78.5 and 4, respectively, with a salt concentration of 0.15 M. CHARMM27 atomic charges and radii were consistently used. The solute-solvent boundary was defined by the van der Waals surface, as molecular surface definition was found to result in significant underestimation of the association rates in some cases.^[57] Finally, electrostatic desolvation effects were empirically accounted for by calculating a desolvation penalty grid around each reactant using the SDA package’s desolvation grid module. To reduce the computational cost of the Brownian dynamics simulations, the full atomic charges of the reactants were approximated by a smaller set of effective charges. These effective charges were carefully derived using the ECM^[58] module of the SDA package to ensure an accurate reproduction of the calculated electrostatic potential. The fitting was performed at the accessible surface (defined by a 4 Å probe) within a 3 Å thick layer extending outward from the reactants’ surfaces.

Structural and functional characteristics of the identified binding sites

In binding site 1, the GLN110 and GLY109 residues are located in domain II, while THR292 is located in domain III. GLN110 resides in the S1 subsite, which is reported to form hydrogen bonds with the substrate, viz., the viral polyprotein recognized and bound by 3CLpro during proteolytic cleavage.^[59] In fact, the inhibitor AT7519 was reported to form a hydrogen bond with GLN110^[13] while the THR292 is known to interact with 3CLpro inhibitors like Remdesivir and Eltrombopag.^[15] Notably, THR292 and GLY109 are near the dimer interface, a region critical for maintaining 3CLpro’s dimer structure.^[60] Inhibitors binding to this site may, therefore, destabilize substrate recognition (via GLN110) or disrupt dimer stability (via THR292/GLY109), indirectly impairing catalysis. In binding site 2, ILE286 and LEU287 are both located in domain III and likely contribute to the domain’s hydrophobic core, which is essential for its structural stability and interactions with other parts of the enzyme.^[61] Interestingly, the super-active STI/A mutant of SARS-CoV 3CLpro involves a mutation at ILE286, suggesting that this region can greatly influence enzyme activity.^[62] As domain III is the primary driver of dimer formation,^[61] inhibitor binding to this site could directly affect the dimerization process. Similarly, in binding site 3, LYS5 is situated in the N-terminal finger region, which plays a critical role in dimerization and interaction with the active site of the opposing monomer.^[63] LYS5 stabilizes the N-terminal domain, which positions the catalytic HIS41.^[64] Within the same site, TYR126 and GLN127 are reported to contribute to the S2 subsite, where TYR126 engages in π-π interactions with substrate aromatic groups.^[64] Inhibitors at this site may, therefore, destabilize the active-site architecture or block substrate anchoring. In binding site 4, CYS145 is a key catalytic residue located in the active site of 3CLpro, directly responsible for the nucleophilic attack on the substrate.^[10] Its mutation to alanine results in complete inactivation of the enzyme.^[65] In fact, many potent inhibitors, particularly covalent inhibitors, target CYS145 to irreversibly deactivate the enzyme.^[65] Within the same binding site, the MET165, located in domain II, forms part of the S2 subsite of the active site and is involved in binding the hydrophobic residues of the substrate.^[66] Binding near MET165 could directly compete with the substrate for binding in the S2 subsite, thereby inhibiting the enzyme’s ability to recognize and cleave the viral polyprotein. Therefore, inhibitor binding to site 4 would directly inhibit the enzyme by physically blocking the active site and preventing the substrate from binding and undergoing catalysis. Binding site 5, however, is located near the interface of domains I and II and involves ALA70 and GLY71 residues in domain I, and SER121 in domain II. ALA70 and GLY71 form part of the S4 subsite, accommodating substrate side chains, while SER121 stabilizes the oxyanion hole during catalysis.^[67] Inhibitors at this site could, therefore, block substrate entry or destabilize transition-state interactions, thereby interfering with the enzyme’s ability to properly position or process the substrate. Interestingly, in binding site 6, THR225 is situated in the linker loop that connects domains II and III, providing flexibility to the enzyme structure. Notably, this loop is distant from both the active site and the dimer interface. Moreover, GLU270 is part of the S1 subsite that is involved in substrate binding via electrostatic interactions with its conserved P1 residues (e.g., glutamine).^[68] Blocking GLU270 would impair substrate recognition, a mechanism observed in peptidomimetic inhibitors.^[7] Finally, binding site 7 involves the ASP153, TYR154, and ILE152 residues that are all located in domain II. Notably, ASP153 and TYR154 are part of a flexible loop that adopts distinct conformations during substrate binding.^[69] Moreover, TYR154 is located at a tight turn in the protein,^[25] suggesting that inhibitor binding could alter the local conformation of domain II. Inhibitors restricting this loop’s mobility may prevent the transition from an open to a closed (active) state, akin to allosteric regulation.

In summary, our analysis of the identified binding sites reveals four distinct mechanisms by which inhibitors can disrupt 3CLpro function. Dimer interface disruption is achieved by inhibitors binding directly to or to critical components of the monomer-monomer interface, exemplified by binding sites 2, 3, and 7. Direct active site inhibition, on the other hand, involves binding site 4, which targets the catalytic CYS145 and the S2 subsite. Furthermore, substrate exclusion or subsite disruption can occur through interactions at binding sites 1, 3, 5, and 6, blocking key subsites or destabilizing substrate interactions. Lastly, allosteric modulation, primarily through binding sites 1 and 7, functions by restricting conformational flexibility of loops or regions near the interface, thereby indirectly impacting activity or dimerization.

It is noteworthy that our docking analysis revealed that a substantial proportion of the identified inhibitors (at least 77%) exhibited the potential to bind favorably (G ranges from -6.4 to-7.8 kcal∙mol^-1) to all seven identified sites, Figure 2(a). This suggests that the inhibitory activity of these compounds might not be attributable to a single mechanism of action. Instead, they could exert their effect through a combination of the proposed mechanisms, such as simultaneously hindering dimerization by binding near the dimer interface and directly interfering with the active site or substrate binding. This polypharmacological potential could contribute to their overall efficacy and might be an important consideration in the design of future 3CLpro inhibitors. Understanding the interplay of these multiple mechanisms for individual inhibitors warrants further investigation.

Dynamical characterization of the inhibitory mechanisms

Following the identification of seven distinct binding sites for potential 3CLpro inhibitors through blind docking, we sought to further explore the potential mechanisms of inhibition, particularly those related to dimer stability. As binding sites 1 (GLN110, GLY109, THR292), 2 (ILE286, LEU287), and 3 (LYS5) involved residues located at or near the dimer interface, we hypothesized that inhibitors binding to these regions could disrupt the crucial dimerization process. To investigate this, we performed Brownian Dynamics (BD) simulations using the SDA software package to assess the effect of inhibitor occupancy at each of the seven identified binding sites on the diffusional association of 3CLpro monomers. Our goal was to quantify how inhibitor binding at different locations might either hinder or facilitate the formation of the active dimeric enzyme. Given the computational cost of simulating all 73 inhibitors across the seven binding sites, we selected the inhibitor with the lowest IC50 value (CHEMBL222234), representing the highest potency, to be bound at each of the seven identified binding sites for the BD simulations. Throughout the BD simulations, formation of the diffusional encounter complex - a state prior to formation of the final bound complex^[53,70,71] was monitored using reaction criteria based on the X-ray structure of wild-type 3CLpro dimer. The rate constant of the formation of the diffusional encounter complex in all simulations was calculated using the Northrup algorithm^[72] within the SDA package. Compared with wild-type monomer-monomer interaction, binding at all seven binding sites results in a decrease of the rate constant of diffusional encounter complex formation by almost 50 % [Figure 3a], which corroborates the importance of ligand binding at these sites. This decrease is also reflected in the radial distribution functions (RDFs) of all seven reactions, where the RDFs retain the same shape of the wild-type dimer interaction, albeit with a decrease in their maxima [Figure 3b]. In a previous work^[36], we provided detailed analysis of the provided RDF for the wild-type monomer-monomer interaction, where it was shown to proceeds through three critical points, two maxima and a minimum, at 22.5, 28.4 31.4 Å, showing enhanced probability, compared with the bulk medium thereby indicating presence of three diffusional intermediate states, encounter complexes^[73] that could constitute precursors to reaching the monomer-monomer closest contact. Binding at all sites maintains the same RDF profile, which could be indicative of pursuing a similar mechanism of interaction. In order to understand the reason behind the modulation of the rate constants and the radial distribution functions upon ligand binding to the seven binding sites, we investigated the structures of the encounter complexes corresponding to the maximal peaks, at ∼ 23 Å, Figure 4. The encounter complex of the wild-type monomer-monomer interactions adopts a structure that is closely similar to that of the wild-type bound dimer, Figure 1(a). Binding at any of the seven sites, however, results in deformation of the encounter complex at maximal peaks, which explains the decrease in the rate constant upon ligand binding. Such disruption could be due to the interception of the trajectory pathways leading to the formation of the wild-type encounter complex, which is the subject of an ongoing study.

Close

Figure 3:

(a) Diffusional rate constants k (in M^-1∙s^-1) of wild-type 3CLpro monomers association (M+M) and wild-type monomer associations with the monomer-inhibitor complexes, ML(n), where the inhibitor L is bound at one of the seven identified binding sites: n=1,2,3…7. A 7.0 Å Cα-Cα distance threshold was used to define the reaction criteria in all association reactions. The average rate constant for all M+ ML(n) reactions is shown in red. (b) Corresponding radial distribution function (RDF) profiles along a reaction coordinate defined by the distance (r) between the Centers Of Mass (COM) of the reactants in each case.

Export to PPT

Close

Figure 4:

Configurations of the 3CLpro monomer-monomer encounter complex at 23 Å COM-COM distance for wild-type 3CLpro monomers association (M+M) and wild-type monomer associations with the monomer-inhibitor complexes, ML(n), where the inhibitor L is bound at one of the seven identified binding sites: n=1,2,3…7. Throughout, chain B of the 3CLpro protein is shown in surface representation with the color scheme used in Figure 1 adopted, and the bound ligand is shown in turquoise. Chain A is shown in green for (M+M) association and in yellow for M+ ML(n) associations.

Export to PPT