Background

Proteins are dynamic molecules that adopt multiple conformations. Protein structure is closely linked to how it responds to internal factors and interactions with the external environment, which in turn determines function in diverse biological processes. The energy landscape of a protein describes the relative free energies of its possible conformations; the probability of occupying a given conformation is reflected by the depth of that state on the landscape. Factors such as ligand binding [1], pH [2], ionic strength [3], or mutations [4] can reshape the landscape and drive structural changes [5]. The physicochemical properties of amino acids also influence these responses. Charged, polar, or hydrophobic residues can affect local interactions, packing, and collective motions [6,7]. In this context, our research focuses on identifying allosteric regions. Allostery refers to a perturbation in one region of a protein that leads to structural and dynamic changes in distant regions, translating a local signal into a global functional response.
Molecular dynamics (MD) simulations are widely used to study protein dynamics and allostery; however, observing large and rare conformational transitions often requires long simulation times and substantial computational resources. In some cases, directly observing the effects of environmental conditions is not feasible, which limits the applicability of purely MD-based approaches [8]. Alternative approaches have been developed to overcome these limitations. For example, elastic network models (ENMs) and methods based on linear response theory have been introduced. In the ENM framework, proteins are represented as coarse-grained networks (typically at the C_α level), enabling the simple yet efficient computation of residue-level fluctuations and their modal decompositions [9–11].
Perturbation Response Scanning (PRS) is a widely used computational method for identifying allosteric residues, grounded in linear response theory [12]. This framework assumes that structural responses near equilibrium are linearly related to applied forces [13]. Accordingly, in PRS, proteins are coarse-grained at the C_α level and subjected to systematic single-residue perturbations to identify residues whose perturbations best reproduce experimentally observed conformational changes [14]. PRS has been successfully applied to various proteins to identify functional hotspots, including enzymes (Caspase-6, ADK) [15,16], molecular motors (Kinesin, Myosin) [16], signaling proteins (CheY, PDZ) [17,18], viral proteins (HIV-1 RT) [16], and many others. Interestingly, these studies have revealed that the residues that maximize the overlap with the identified conformational change are not limited to surface-exposed positions; depending on the system, they may include exposed residues, partially buried interface residues, or buried residues near functional cores. However, PRS is practically limited to single-residue perturbations, and it may not detect conformational changes caused by the simultaneous and collective effects of multiple residues. When it comes to domain or hinge-movements, or rearrangements dependent on environmental conditions, multiple coordinated perturbations may be required to trigger the observed structural rearrangements [16].
To overcome this limitation, we introduce the Multiply Perturbed Response (MPR) method [19]. This approach extends PRS by applying forces simultaneously to multiple residues, capturing conformational transitions driven by collective residue effects. When only static initial and final structures are available, the protocol can be executed using the structure-based formulation alone. When MD trajectories are provided, the trajectory-based method captures equilibrium conformational fluctuations, improving the precision of the predicted conformational response. We note that for proteins without experimentally determined structures, MPR can be applied using predicted structural models. In this case, alternative conformational states may be generated using AlphaFold/ColabFold with different MSA depths, RoseTTAFold, ESMFold, or other validated structure-prediction workflows. However, because the input conformations are computational predictions, the resulting MPR residues should be interpreted as hypotheses for future experimental testing rather than experimentally validated allosteric sites.
MPR integrates enumeration- and optimization-based strategies to explore solutions across different numbers of perturbed residues (k ≥ 1). Enumeration enables exhaustive identification of optimal residue combinations for small k values, while optimization-based formulations allow efficient exploration of higher-order perturbations without prohibitive computational costs. As a result, the method substantially reduces computational expense compared with fully atomistic approaches while retaining sensitivity to collective allosteric effects. Using open-source software and widely used bioinformatics libraries increases the accessibility and applicability of the method.