Toward the Design of Allosteric Effectors: Gaining Comprehensive Control of Drug Properties and Actions

Generic Computational Framework for the Identification and Design of Allosteric Sites and Effectors

We developed and implemented here a computational framework, which is based on the directed design protocol that we recently proposed for allosteric drug development.⁵ The framework identifies allosteric sites on a target protein and yields a set of allosteric ligands that bind to a designated site. Figure 1 illustrates key components of the framework.

Structure-Based Statistical Mechanical Model of Allostery

The computational framework developed here is based on the SBSMMA^18,19—a physical model that describes causality of allosteric communication and quantifies it via free energy of the allosteric signaling (in real energy units, kcal/mol). The allosteric signal can originate from different types of perturbations, including but not limited to ligand binding and mutations (Figure 1a). The SBSMMA consists of three major steps. First, it takes a protein structure as an input and constructs a C_α elastic network representing it, the dynamics of which is simulated on the basis of the C_α harmonic model (see Methods for the energy function in eq 1). Ligand binding is modeled as a perturbation of the energy function through an inclusion of harmonic restraining terms among binding residues (shown as yellow spheres), which simulates a stabilization effect due to interactions with a ligand. Here, the protein dynamics in native and perturbed states are characterized by respective sets of low-frequency normal modes. In other implementations, protein dynamics may also be explored using principal components of the covariance matrices of trajectories obtained in molecular dynamics simulations. Second, the normal modes are used as an input to calculate the per-residue allosteric potential in free and perturbed states (see eq 2 in the Methods), allowing us to evaluate elastic work exerted on a residue as a result of the dynamics of neighboring residues in both states. Lastly, the allosteric potential is used for calculation of the partition function (eq 3 and Methods) to estimate the per-residue free energy of the allosteric signal Δg on targeted residues originated by the ligand binding (eq 4 and Methods). To estimate the work caused by purely allosteric signaling, the per-residue allosteric modulation Δh (kcal/mol) is calculated as the difference of Δg from the average value over all residues of the protein chain (eq 5 and Methods). A positive modulation value (Δh > 0) indicates work exerted on involved residues, which might lead to conformational changes, whereas a negative value (Δh < 0) points to decrease in dynamics of the affected residues, suggesting overstabilization of corresponding parts of the protein. To compare the effects at different sites, allosteric modulation can be evaluated by averaging the per-residue Δh over all residues forming the site of interest.

Reverse Perturbation and APM

The framework utilizes the reverse perturbation^27,38 and allosteric probing map (APM²⁸) approaches previously proposed and developed by us. A generic reverse perturbation concept³⁸ was implemented for a large set of proteins²⁷ to detect latent/known allosteric sites using information on functional sites and for characterizing their regulatory functions and drug design potential.^8,22,24 Noteworthy, an independent experimental work¹⁶ demonstrated a presence of bidirectional allosteric communication between the ATP-binding and PIF pocket in kinases. Alterations of protein dynamics associated with the reversible signaling across both sites were also characterized using HDX-MS and SBSMMA.¹⁷ Recently, the reverse perturbation concept has been adapted and utilized by other groups for various tasks.^39,40 Here, we use the reverse perturbation approach for identifying distant protein residues that are allosterically coupled to a site of interest, including the functional or active site of the protein (Figure 1b, left). First, a binding perturbation to a known functional site is performed using a protein structure and a set of known functional site residues as inputs. Then, distant residues (with C_α–C_α distance at least 11 Å to all functional site residues, the cutoff distance used in the elastic network model of protein) exhibiting the strongest allosteric response (i.e., high-magnitude modulation values) are shortlisted. The exhaustive APM is computed by simulating binding of a small probe ligand to every three-residue segment/triplet in a protein chain and evaluating the resulting modulation at each residue (Figure 1b, right). It allows identification of positions in a protein structure that originate a desired signaling to any single residue or residue group in the targeted site. The computation of APM requires only a protein structure as an input and information on a functional site is not necessary.

Depending on the tasks, one may shortlist segments/residues that elicit required modulation to a distant functional site known a priori or specific residue(s) that may play a functional role. Importantly, due to the extensivity of the free energy, the exhaustive APM may be utilized to explore and estimate the allosteric signal as a sum of contributions from various combinations of probed three-residue segments. The APM and reverse perturbation are typically used as complementary approaches. The former can be used to evaluate allosteric signaling that originated from binding to individual segments, including residues strongly linked to the functional site determined by the latter approach. Analyzing allosteric communication linking regulated and regulatory sites in opposite directions, both APM and reverse perturbation methods can be used independently or in conjunction. A combination of top-ranking distant residues/segments shortlisted by the reverse perturbation and/or APM method yields a list of residues strongly coupled to the remote functional site, delineating protein regions/locations consisting of latent allosteric sites with varied residual composition.

Generation of Site–Ligand Repertoire

In this component, residues in allosterically linked locations identified by the reverse perturbation and APM approaches form a basis for the screening or design of the corresponding ligands. Structure-based molecular docking of compound libraries and de novo molecular design can be used to sample interactions of various ligands to (non)overlapping sets of binding residues with sufficient estimated binding energy (Figure 1c). Virtual screening is demonstrated here with AutoDock Vina⁴¹ (see also Methods for details), and docking results with very low estimated binding free energy were discarded. Since the reverse perturbation and APM approaches are based on the analysis of protein dynamics on a per-residue resolution, they reveal and quantify signaling to and from single residues in sites/pockets apparent in static crystal structures as well as residues that may form relatively transient binding sites. Noteworthy, commonly used definitions and prediction for the so-called cryptic sites,⁴²⁻⁴⁵ emerging in protein dynamics as transient binding pockets, generally utilize properties of known binding pockets of functional sites without explicit characterization on the allosteric signaling and not considering many other potential allosteric sites with properties distinct to those of functional pockets.^5,8 Recognizing the signaling as a prerequisite for defining the site as allosteric and in a full agreement with the original work⁴⁶ by Bowman and Geissler, which shows a multitude of distant regulatory sites capable of originating allosteric signaling to functional sites, including the “allosteric cryptic sites”, we evaluate the allosteric modulation and binding strength characterizing each ligand and corresponding site (see below).

Quantification of Signaling and Binding Properties

The simultaneous evaluation of corresponding pairs of binding sites and ligands is an important part of the framework (Figure 1d). It estimates the energetics of allosteric signaling induced by an alteration of the protein dynamics upon binding, which should be optimized in parallel with binding affinity and other relevant properties. Shortlisted pairs of binding sites and ligands with required specifics of allosteric modulation and affinity, allosteric site and effector pairs, serve as starting points for further derivation and design of the chemical structures of effectors and residual composition of binding sites. It is followed by evaluation and fine-tuning of the signaling, binding, and other properties. Importantly, statistics and descriptors on characteristics of known allosteric drugs/drug candidates are utilized to complement the screening/design protocols and for shortlisting of site–effector pairs within the framework.

In the current implementation of the generic framework outlined above, we used the advantages of well-characterized drugs/drug candidates to efficiently utilize the information and benefits of characterized “drug-likeness”, pharmacokinetics, and effects of available therapeutics for new modes/mechanisms of actions or targets. Here, molecular docking of the Drug Repurposing Hub (DRH)⁴⁷ revealed approved and clinical trial drugs binding to latent allosteric sites on protein(s) that are not original targets for these drugs. Noteworthily, while the pharmacological and therapeutic profiles of repositioned orthosteric drugs are typically well-characterized, allosteric and orthosteric compounds/sites can exhibit distinct physicochemical and molecular properties,^8,48 which should be taken into account. On the basis of known drugs providing a source of drug-like molecular fragments for fragment-based design, new “allosteric drug-like” effectors were derived for various (non)overlapping allosteric sites with diverse features and requirements.

We have previously argued for the simultaneous design and adjustment of site–effector pairs by evaluation and optimization of allosteric signaling and other properties, termed as the directed design protocol.⁵ Here, the framework based on the directed design protocol yields allosteric site–effector pairs after optimizing them toward desired characteristics, including but not limited to (i) sufficient “drug-likeness” with respect to orally administered medicines and (ii) “allosteric drug-likeness” based on stereotypical properties of known allosteric drugs/effectors;⁸ (iii) strong binding affinity to regulatory sites; (iv) required allosteric modulation of the protein function; and (v) specific long-range signaling to precisely target protein residues/sites. The framework is demonstrated using M^Pro and K-Ras^G12D as case studies, each showing important aspects in the discovery and design of allosteric medicines. For example, the M^Pro case study illustrates how one can derive new compounds with molecular properties typical for allosteric effectors and the required level of drug-likeness, based on fragments of previously known orthosteric drugs instrumental for originating targeted allosteric responses. In the latter case study, a unique approach on the design of new effectors highly specific in the distinct regulation of switch regions I and II was presented.

Revealing Potential Allosteric Sites Including the Domain I/II Hinge Region

The M^Pro harbors a catalytic dyad formed by His41 and Cys145 in the catalytic cleft situated between 2 six-stranded β barrels known as domains I (residues 10–99) and II (residues 100–182). Domain III (residues 198–306) with five α helices is involved in the dimerization of the enzyme, which is functional only as homodimers. In an X-ray crystallographic screen, clinical trial drugs pelitinib and AT7519 were discovered to bind outside of the catalytic cleft and showed moderate antiviral activity in viral reduction assays.⁴⁹ However, the mechanism underlying the antiviral activity, specifically the involvement of allosteric communication, was not elucidated from the crystallographic structures. Information on the existing allosteric inhibitors, including their binding sites, allowed us to first evaluate to what extent the SBSMMA, a key component of the framework, accounts for the allosteric modulation. Binding residues within 4 Å from any heavy atom of pelitinib and AT7519 were obtained from the respective ligand-bound crystal structures (PDBs: 7axm and 7aga). Then, SBSMMA was used to measure the allosteric modulation caused by perturbations of the binding sites (see Methods for detail) in a substrate and inhibitor-free M^Pro dimer structure (PDB: 7jp1). Both pelitinib and AT7519 cause negative allosteric modulation in the helical portion of domain III, suggesting stabilization near their binding sites (Figure S1). We also observed positive modulation at the distant catalytic cleft and domain III dimer interface (0.42 and 1.2 kcal/mol), suggesting conformational changes upon AT7519 binding. For pelitinib, the interface is also positively modulated but not the catalytic cleft. Overall, we found that pelitinib and AT7519, effectors identified from previous experimental screening, allosterically destabilize the distant functional sites via alterations in the structure and dynamics, affecting the enzyme activity and/or dimer formation.

The computational framework was applied to demonstrate a rational approach (Figure 1) for the identification and design of latent allosteric sites and prospective effector drugs. Reverse perturbation by simulating substrate binding at catalytic clefts showed positive modulation exerted predominantly to residues in the domain I/II hinge and, to a lesser extent, the domain III interface (Figure 2a). The allosterically coupled regions identified by the method also overlap with the AT7519 effector binding site, corroborating the allosteric mode of action. The APM for the M^Pro dimer (chains A and B) was computed to obtain an exhaustive per-residue description of allosteric signaling upon simulated probe binding to each residue triplet (Figure 2b). It revealed major patterns and specifics of allosteric communication across regions, domains, and monomers, for example, (i) probe binding positions in domain I often result in positive modulation in domain II of the same monomer; (ii) probe binding to domain III positively modulates domain I/II of both monomers; and (iii) stabilizing the flexible connecting loop between domain II/III (residues 190–200) of each monomer (by interactions with probes) may result in decreased dynamics of residues in domain III (negative modulation); at the same time, it sends allosteric signals to domains I/II of chains A and B that may lead to conformational changes. Using the APM, we identified residue triplets that can induce the strongest allosteric modulation to distant catalytic cleft residues upon perturbation by a small probe. Two clusters of triplet positions, which might form allosteric sites, were found to exert positive modulation in the catalytic cleft, indicating conformational changes (Figure S1). The first location L1_MPro is centered at residues 14–16 and 99–102 near the domain I/II hinge (opposite face of the monomer relative to the catalytic cleft, indicated in Figures 2a and S1). The second area L2_MPro located in the vicinity of residues 197–203 in domain III causes weaker modulation in the catalytic cleft (Figure S1). Together, the results from APM and reverse perturbation simultaneously point to the L1_MPro allosteric location at the domain I/II hinge, consistent with previous experimental and computational findings.^50,51

Toward the Repurposing of Drugs for Allosteric Regulation of the M^Pro Catalytic Cleft

Molecular docking was performed on L1_MPro to obtain binding poses of existing drugs from the DRH⁴⁷ (with about 6800 clinical and launched drugs). For each ligand, several docked poses were considered to elucidate how compositional differences of binding sites (all residues within 4 Å from any heavy atom of the ligand) may fine-tune the signaling, providing the information for designing new allosteric effectors and sites. Prioritizing those with sufficient interaction strength (within the top 20% of all evaluated binding free energies), we obtained a repertoire of site–effector pairs with varied combinations of residues bound by different effectors. Subsequently, respective allosteric modulations at the distant catalytic cleft (averaged per-residue Δh) originated by corresponding site–effector pairs were calculated using SBSMMA.^18,19 The results show that site–effector pairs at L1_MPro induce a range of positive modulation values (Figure 2c). The majority of site–effector pairs lead to modulation values of about 0.3 kcal/mol at each catalytic cleft residue, and some pairs cause even stronger modulation than the known effector AT7519 described above (Figure S1). This result illustrates that while ligands may bind strongly to a location, the perturbation caused by the ligand binding may not necessarily induce a large modulation at targeted distant sites. The two-step action of allosteric effectors—binding and sending a modulatory allosteric signal—is an important distinction in comparison to traditional orthosteric drugs, which works via competitive binding to a functional site. The strengths of binding affinity and of allosteric modulation depend individually on characteristics of the targeted protein—the former is largely determined by the intermolecular interactions and the latter resulted from an alteration of protein dynamics. Both binding and signaling aspects are chiefly independent of each other, being determined by a number of characteristics, such as structure, size, and oligomerization state of the protein, structures and compositions of its allosteric and functional sites, and others. The binding and signaling values may, thus, be individually controlled and adjusted in designing of prospective allosteric drugs. The site–effector repertoire can be further analyzed to identify those with optimal modulation and binding capability, as well as other properties. For example, itacitinib, SR-1664, and entrectinib (compounds 1–3) showed large modulations and binding energies. Compound 2 is an antagonist for the PPARγ nuclear receptor, whereas compounds 1 and 3 are kinase inhibitors. Details on interacting residues in respective binding sites of compounds 1–3 and elicited per-residue modulations can be found in Figure S1 and Table S1. Figure 2d illustrates that long-range allosteric signaling induced by binding of compound 2 (SR-1664) at L1_MPro of chain A (which also interacts with chain B residues) causes positive modulations at both the catalytic cleft and dimer interface, indicating potential conformational changes that disrupt the protein function and dimerization. The parts A and D of Figure 2 also illustrate that allosteric communication is reversible: perturbation of a functional site originates a signaling to an allosteric site and vice versa.

Elucidating Key Residues in Allosteric Communication for the Design of New Effectors

Prospective effectors should be able to initiate long-range signaling to a targeted site upon binding to an allosteric site with a sufficient affinity. Figure 2c shows that while some ligands may bind tightly to a protein site, the perturbation caused by their binding induces rather weak modulation. Therefore, methods for obtaining and tuning of site–effector pairs at the molecular level should consider two important aspects—affinity and allosteric signaling,^5,8 which are chiefly determined by specifics of perturbed site residues and intermolecular interactions. We focus here on elucidating key determinants of the strength of allosteric signaling exerted to the catalytic cleft. The catalytic cleft modulation is largely uncorrelated with the number of perturbed residues n at the distant L1_MPro (Figure S2). On the other hand, the heatmaps in Figure 2e show that some site–effector pairs with specific interactions/contacts can induce a large allosteric response in the catalytic cleft. Comparing the frequencies of contact residues in 100 site–effector pairs ranked top in modulation versus the lowest ones (using n = 10 and 16, for instance) shows that perturbing residues Lys12.B, Lys97.B, and Gln306.B (in chain B) is crucial for inducing stronger allosteric signals. Full contact analysis with n from 10 to 16 also revealed other key residues including chain A residues Ala70.A and Gly71.A (Figure S2). Noteworthy, the above residues are also involved in the interactions with compounds 1–3 obtained from the DRH.⁴⁷ Specific residues obtained in the analysis of the binding sites further emphasize the demand for mutual design of site–effector pairs, in which adjustment and tuning of both, the site and the effector components, are of similar importance.

To obtain effectors for diverse allosteric sites, which, in principle, can exist at any location on a protein structure, one may explore beyond the chemical space populated by existing drugs and drug candidates that work orthosterically. We previously showed in a statistical analysis that allosteric and orthosteric drug candidates exhibit distinct properties. The statistics and knowledge of the salient properties of known allosteric compounds can be used to prioritize promising candidates. For example, allosteric effectors generally contain more aromatic rings and are typically flatter.⁸ Using compound 2 (SR-1664), we show how new effectors with “allosteric drug-like” properties can be derived from existing drugs. We obtained molecular fragments of compound 2 using the retrosynthetic combinatorial analysis procedure (RECAP) method⁵² (Figure S3) and selected the 2-phenylbenzoic acid part, which binds near key residues Ala70 and Gly71, for subsequent substructure-based analogue search (Figure 3a). The fragment is queried against the large-scale ZINC20 database⁵³ and yielded 2531 analogues for further investigation. Pairwise Tanimoto similarity shows the majority of analogue pairs with a coefficient at about 0.2, suggesting diversity in chemical structures (Figure S3). The quantitative estimate of the drug-likeness (QED), providing a continuous measurement of compound quality and potential as oral drugs,⁵⁴ allows elimination of those with very low likelihood of success in development. The criteria on flatness and aromatic ring contents were employed here for a preferential selection of “allosteric drug-like” effectors (Figure 3b). To assess molecular shapes, we first calculated the principal moments of inertia (PMI) of all obtained analogues (Figure S3) and used the second normalized PMI ratio (NPR2) as a measure for flatness.⁵⁵ The number of aromatic rings normalized by the molecular weight was also evaluated. Additional data on the distributions from the normalized PMI analysis and on the molecular properties are available in Figure S3. After discarding potentially promiscuous analogues using the pan assay interference compounds (PAINS) filter⁵⁶ and those containing chemical moieties with undesirable toxicity and pharmacokinetics using the Brenk filter,⁵⁷ followed by a selection of those satisfying the above criteria (QED ≥ 0.5, AROM/MW ≥ 6 kDa^–1, NPR2 ≤ 0.9; Figure 3b), 708 analogues of compound 2 were collected.

We next performed virtual screening (i.e., molecular docking) of the analogues on L1_MPro followed by a quantification of modulations and binding energies. Stratifying site–effector pairs by structural scaffolds of analogues allows one to explore how modifications within same scaffold types can lead to different binding free energies and modulations. The Bemis–Murcko method was used to characterize scaffolds/frameworks defined as the union of ring systems and linker atoms.⁵⁸ The modulation and estimated binding free energy of every site–effector pair belonging to scaffolds 1, 2, and 6 are shown in Figure 3c–e. Results on 6 selected scaffolds out of a total of 170 are shown in Figure S4. Analogues featuring a biphenyl scaffold (scaffold 1), which is the simplest and most populated one, display wide spectra of modulations, binding energies, and QED values, reflecting a potential diversity of structural analogues that can be derived. Analogues belonging to scaffold 2 show uniformly lower QED values, but many of the site–effector pairs result in modulation and affinity stronger than those from scaffold 1. Comparison of analogues 2a and 2b indicates that minor modifications in effector structures can lead to different site–effector pairs with distinct binding modes (see Figure S4) and respective large changes in allosteric response (Figure 3d). Conversely, scaffold-6 analogues generally feature higher QED (more oral drug-like) and show a different distribution relative to scaffold-2—typically weaker modulation and binding strength. While analogue 2d binds solely to residues in chain A, 2c binds to both monomers and forms a salt bridge with Lys12 in chain B (Figure S4), which may induce strong positive modulation (Figure 2e).

The major result of the design approach, starting from targeting the identified allosteric sites with orthosteric drugs to deriving new effectors, is that strong binders and those showing “drug-like” characteristics established based on orthosteric ligands are not necessarily optimal as allosteric ligands. The latter should also satisfy requirements on specifics of allosteric signaling and perturbations (site–effector pairs) required for inducing the desired signaling. Depending on the tasks, therefore, the molecular characteristics of allosteric effectors should guide explorations of the chemical space, aimed at obtaining new structural species with both necessary binding affinity and potential for exerting an allosteric signal.

Identification of Allosteric Sites Specific for the Switch-I and Switch-II Regions

K-Ras is a monomeric GTPase with a catalytic domain that is highly conserved in the Ras family. Conformational changes in the flexible switch I (SW-I, residues Asp30-Asp38) and switch II (SW-II, residues Gly60-Glu76) regions of the catalytic domain play instrumental roles in various aspects of K-Ras regulation and function.³² SW-I/II and the P loop constitute the catalytic site, and the former (mainly SW-I residues) interacts and activates different downstream effector proteins. Moreover, K-Ras activity is tightly regulated by guanine nucleotide exchange factors (GEFs) and GTPase-activating proteins (GAPs) which interact with SW-I/II to facilitate GDP dissociation and GTP hydrolysis, respectively.⁵⁹ We focus here on the identification of allosteric sites capable of emanating specific and strong signaling to the SW-I and/or SW-II regions of the oncogenic K-Ras^G12D (4B splice isoform), relative to the wildtype. To this end, comprehensive APMs were computed for the catalytic domains of K-Ras^G12D (Figure 4a) and K-Ras^WT (Figure S5). The results show probing residue triplets in the C-terminal lobe (residues 87–166, Figure 4b) generally causes large positive allosteric modulation in SW-I/II of both structures, suggesting conformational changes at these functional regions. From the APMs, we also obtained distant probed positions that originate much stronger modulation in the SW-I/II region of the K-Ras^G12D variant, compared to the modulation in K-Ras^WT (tabulated in Figure S5). Locations L1_K-Ras and L2_K-Ras were delineated based on these positions (Figure 4c), potentially forming allosteric sites to remotely modulate the protein dynamics of distant SW-I and SW-II of K-Ras^G12D, respectively. L1_K-Ras (colored in violet and magenta) and L2_K-Ras (red and magenta) are overlapping; perturbing residues 136–138 (magenta) is instrumental for sending allosteric signals to both switch regions. To our knowledge, both locations have received little attention for therapeutic developments. Emerging experimental data have shown that a small molecule termed KAL-21404358 (KAL)⁶⁰ and a monobody known as NS1⁶¹ bind locations overlapping with L1_K-Ras and L2_K-Ras, respectively, and allosterically inhibited K-Ras function. Here, SBSMMA shows that the binding of KAL and NS1 can remotely modulate the dynamics of SW-I/II and other distal regions (Figure S5).

Toward Designing New Effectors for Allosteric Regulation of the Switch-I Region

We next performed molecular docking of known drugs and drug candidates on L1_K-Ras, which is allosterically coupled with the SW-I region. Using the same protocol as that for the M^Pro, site–effector pairs were evaluated for modulation values at the remote SW-I region and predicted binding free energies at L1_K-Ras (Figure 4d). Compounds 4–6 (midostaurin, IKK-16, and WAY-207024, respectively) exhibiting high modulations and binding free energies are shown here. The binding of 5 and 6 to L1_K-Ras induces stronger allosteric modulation to SW-I/II than to KAL (Figure S5). For example, compound 6, an antagonist for the gonadotropin releasing hormone receptor, binds to and stabilizes L1_K-Ras. At the same time, it causes extensive positive modulation in the N-terminal lobe containing functional sites, especially the SW-I/II regions (Figure 4e). Residues important for signaling were elucidated by comparing their contact frequencies in site–effector pairs ranked as the top and the lowest in terms of modulation (Figure 4f). Binding sites with n = 8 and 13, which result in modulation values ranging from 0.2 to 0.6 kcal/mol and from 0.45 to 0.7 kcal/mol, respectively, are shown here as examples (see Figure S2 for complete data). The heatmaps show that in cases of lower numbers of contact residues (indicative of smaller ligands, such as n = 8), intermolecular interactions in top-ranking site–effector pairs are rather confined to specific preferential residues, for example, Pro110 and Arg162. The critical role of Pro110 and Arg162 in allosteric signaling is corroborated by several experiments. Both residues form the binding site, dubbed as the P110 site by authors, for KAL.⁶⁰ Using deep mutational scanning, it was also shown that mutations at Pro110 and Arg162 disrupted the binding of Raf-1 effector protein from a distance.⁶² In another example (n = 13, Figure 4f), perturbing residues His94 and Glu98 increases the modulation at SW-I at the expense of interactions with Lys165, showing the fine-tuning of signaling by mutual adjustments of the sites and effectors.

Subsequently, we explored the design of new site–effector pairs at L1_K-Ras on the basis of compound 6. Analogue search yielded about 5000 compounds containing a 6-(1-piperazinylmethyl)quinoxaline group (red, Figure 5a). Similar to the case of M^Pro, we used several criteria to prioritize flat analogues with a higher aromatic ring content and discarding those with low likelihood as orally administered drugs (AROM/MW ≥ 7 kDa^–1, NPR2 ≤ 0.95, QED ≥ 0.65; Figures 5b and S6). Molecular docking followed by SBSMMA-based evaluation of modulation was performed for all site–effector pairs obtained from 841 shortlisted analogues. The results from site–effector pairs belonging to scaffolds 1 and 6 revealed repertoires of analogues causing a range of positive modulations at SW-I upon interacting with different combinations of residues at L1_K-Ras (Figure 5c). Comparison among analogues 6a–d shows that interactions of 6a and 6c with both allosterically coupled residues Pro110 and Arg162 described above (Figure S7 and Table S2) likely contribute to larger modulation values and binding free energies compared to 6b (Figure 5c) and 6d (Figure 5d). While all analogues show slightly weaker modulations (with the highest being 0.6 kcal/mol) compared to that of precursor 6 (0.8 kcal/mol), other aspects, such as the estimate of drug-likeness (QED of 6 = 0.36), were improved considerably in these analogues. Similar to the case of M^Pro, we showed that modifications of chemical structures and interacting residues of the site–effector pairs allow optimization of the molecular properties and functional outcomes, such as affinity to the binding site and induced allosteric modulation. One can conclude that the directed design approach implemented here is universal, and it can be used in future protocols for designing new allosteric drugs/effectors for proteins with different functions and intrinsic mechanisms of their regulation.

Tuning the Signal Specificity of K-Ras Effectors for SW-I/II Regions

As one of the critical hubs in cellular signal transduction, activated Ras proteins regulate and, at the same time, are regulated by other proteins, such as GEFs and GAPs, via different interaction interfaces formed by various residues in SW-I/II. In principle, precise modulation of conformational dynamics at interface sites/residues would allow one to specifically target dysregulated interaction(s) mediated by a signaling protein with minimal interferences in the other pathways. Therefore, the specificity of intra- and intermolecular signaling, a quintessential advantage of prospective allosteric effectors, should be considered as one of the major goals in their design. Depending on a required modulation and targeted residue(s), effectors inducing broad modulation or those modulating only specific sites/residues may be favored one over another (Figure 6a). To characterize the signal specificity of compound 6 and its derived analogues, we evaluated the modulation at both switches (Figure 6b). Binding of compound 6 (magenta) results in a modulation of 0.8 kcal/mol at switch I and switch II, which may affect interactions with different proteins. Noteworthy, compound 6 interacts with Ser136-Gly138, which can originate signaling to both switches as revealed by the APM (Figure 4c). For the analogues, while the majority of site–effector pairs at L1_K–Ra causes larger modulation in SW-I than SW-II, varied modes and extents of switch specificity can be observed. Analogues 6e and 6f (orange and blue) are exemplified here. They differentially modulate SW-I and SW-II upon binding to subsets of residues bound by compound 6. Binding of 6f induces large modulation at switch I but not switch II, whereas the opposite was observed for 6e. The former binds three additional residues Met111, Ser136, and Glu162 compared to the latter (Figure 6c). As described above, Glu162 and, to a lesser extent, Met111 are both instrumental for signaling to SW-I (Figures 4f and S2), indicating that perturbing these residues generally increases the dynamics of SW-I residues and likely exerts the opposite effects in SW-II. Information about the binding sites is tabulated in Table S2. Altogether, the results here demonstrate that fine-tuning of the effectors and the composition of perturbed residues allow one to change the signal specificity of site–effector pairs toward distinct functional sites.

Structure-Based Statistical Mechanical Model of Allostery

The SBSMMA^18,19 was used here to quantify the energetics of allosteric communication at the per-residue resolution. The physical model was thoroughly benchmarked in multiple studies^25,27 accounting for the causality of allosteric signalling, and it describes the allosteric mechanisms caused by general perturbations P, including ligand/probe binding,^{17,19,22−24,28} mutation(s),^22,25 order–disorder transitions,²⁶ post-translational modifications,²² and their combinations. The model consists of three major steps. First, it takes a protein structure as an input to represent it via the Elastic Network Model (ENM) and to construct the C_α atom-based harmonic models for both native 0 and perturbed P states of a protein. The energy function for the perturbed P state is

1

The first term corresponds to the native state with d_ij and d⁰_ij as interatomic distances between C_α atoms in respective generic r and reference r⁰ structures. The k_ij is the distance-dependent force constant that decays as with a global distance cutoff at 25 Å, with a summation over pairs of neighboring C_α atoms i and j. The second term is an additional perturbation harmonic term for the ligand binding, where S contains a list of C_α atoms representing a binding site. Ligand binding is indirectly modeled with a stiffening α_P factor at 100, affecting residue pairs forming S and resulting in a local stabilization effect. Orthonormal modes e_μ characterizing the configurational ensemble of a native/perturbed state are obtained from the Hessian matrix K = ∂²E/∂r_i∂r_j. For each state, allosteric potential which evaluates the effect of a perturbation on a residue i as elastic work exerted due to changes in the displacements of proximal residues (with C_α–C_α distance not more than 11 Å) described by the first ten normal modes e^(P)_μ and e⁽⁰⁾_μ reads

2

where and σ = (σ₁, ..., σ_μ, ...) is a set of Gaussian distributed amplitudes with variance 1/ε_μ,i. As the generic displacement of a residue i is r_i(σ) - r⁰_i = ∑_μσ_μe_μ,i, vector σ can be considered as a configurational state of residue i. In the last part, integrating over the ensemble of all possible configurations σ of a residue, the partition function on the level of a single residue is obtained as

3

The per-residue partition function z_i allows one to estimate the free energy associated with an allosteric signal g_i = −k_BTlnz_i. The free energy difference (kcal/mol) reads

4

which quantifies the work exerted on residue i as a result of perturbation P. A positive Δg^(P)_i value indicates an increase of configurational work exerted on residue i which may result in alteration of the protein structures and dynamics; a negative value suggests overstabilization of the residue. Lastly, allosteric modulation Δh^(P)_i which evaluates exerted work due to purely allosteric signaling is obtained as a difference Δg^(P)_i from its mean over all residues of the protein chain containing i

5

A positive modulation value points to potential conformational changes, whereas a negative one indicates that conformational changes are precluded. Importantly, per-residue Δh^(P)_i values with magnitude above thermal fluctuations k_BT (0.6 kcal/mol) reflect significant allosteric communication. Generally, the k_BT value may be used to discriminate between an efficient signal and signaling noise due to thermal fluctuations. Nevertheless, in some cases, Δh^(P)_i values close to k_BT may still be considered in a homogeneously affected protein region due to the extensivity of the free energy that may result in strengthening of the signal by small additional contributions. While large magnitudes of modulation values may suggest strong allosteric effects on a protein function, the functional outcomes caused by binding of effectors should be validated in experiments or compared with known experimental data from other works. To compare allosteric modulation at different sites, a modulation on a site of interest S is given by where residue i∈S. We have previously shown that perturbations simulated at ligand-free input structures initially obtained from separate crystallographic experiments in apo and holo states typically lead to comparable results between them due to the overlapping conformational ensembles sampled in both states. Moreover, the SBSMMA may be applied to protein structures with or without apparent binding sites/pockets, as the model quantifies the per-residue energetics of signaling to and from single residues regardless of specific structural properties or dynamics of protein surfaces. The SBSMMA has been implemented in the AlloSigMA Web server²⁸ and AlloMAPS database,²⁹ widely used by experimental and theoretical researchers as explorative and analytical tools. Users should expect a computing time ranging from minutes to several hours to obtain downloadable data from AlloSigMA, which depends on the number of residues in the target protein and the specification of the tasks. More details on the computing time can be found in the Benchmarking section in the Tutorial page.

Identification of Allosterically Linked Locations

Crystal structures of M^Pro, K-Ras^G12D, and K-Ras^WT (PDB: 7JP1, 6GOF, and 6GOD) were used by the framework to identify potential allosteric locations. Residues constituting the SW-I/II regions in K-Ras and the catalytic cleft and dimer interface of the M^Pro were collected from the literature.^66,67 In reverse perturbation, substrate binding was simulated, and Δh^(Bind)_i was measured for all residues. The APM is obtained by simulating binding of a small probe ligand to every consecutive three-residue segment along a protein chain and evaluating the modulation Δh^(Probe)_i at each residue i. A combination of top-ranking residues/segments, with the largest |Δh| and with C_α–C_α distance at least 11 Å (the operational cutoff distance for the elastic network model of the protein) to all targeted site residues, is compiled by the reverse perturbation and/or APM methods. It presents a list of residues strongly coupled to the targeted site, allowing delineation of protein regions/locations containing allosteric sites with different compositions of residues.

Molecular Docking and Fragment-Based Design of Effectors

DRH⁴⁷ containing approximately 6800 approved and clinical trial drugs with molecular weight not more than 600 Da was used for molecular docking. Structures of all compounds in PDBQT format were converted from SMILES strings using OpenBabel⁶⁸ and protonated at pH 7. The protein structures were prepared using AutoDockTools⁶⁹ prior to docking. Molecular docking was performed at the geometric center of shortlisted strongly modulating residues/triplets on the crystal structures with a cubic box of lengths of 18 and 20 Å for K-Ras and M^Pro, respectively. AutoDock Vina⁴¹ was used with default scoring functions and parameters, providing an estimation of the free energy of binding. For each ligand, up to 10 docked poses were examined. Allosteric site residues were delineated based on a distance cutoff of 4 Å from any heavy atom of a bound ligand, which were then used for computing the allosteric modulation Δh^(Bind)_S at a distant regulated site using the SBSMMA formalism. Substructures/fragments of shortlisted molecules from DRH were obtained with the RECAP method, which performs fragmentation based on chemical knowledge. Subsequently, the SMILES strings of selected fragments were used to carry out an analogue search on the ZINC20 database.⁵³ The RDKit package⁷⁰ was used to flag and remove analogues using implementations of the PAINS and Brenk filters, to obtain Bemis–Murcko scaffolds⁵⁸ of derived analogues and to evaluate properties of compounds. In particular, the pairwise Tanimoto similarity was calculated using the Morgan fingerprints implemented in RDKit. The default average descriptor weights were employed in the QED calculation. The calculation of NPMI values was performed by the Descriptors3D module in RDKit after the embedding of a 3D conformation, which was set to prioritize known experimental torsional angles. The SBSMMA was used to evaluate the modulation caused by derived analogues. Details on compounds, binding sites, and modulations are tabulated in Tables S1 and S2.