Lysozyme allosteric interactions with β-blocker drugs

Lysozyme allosteric interactions with b-blocker drugs Michael González-Durruthy a , Ramón Rial a , Zhen Liu b , Juan M. Ruso a, ⇑ a Soft Matter and Molecular Biophysics Group, Department of Applied Physics,University of Santiago de Compostela, 15782 Santiago de Compostela, Spain b Department of Physics and Engineering, Frostburg State University, Frostburg, MD 21532, United States article info Article history: Received 26 July 2022 Revised 8 September 2022 Accepted 11 September 2022 Available online 16 September 2022 Keywords: Betablocker drugs Lysozyme Allosteric interactions Molecular docking abstract Effective and reliable prediction of allosteric molecular interactions involved in protein-ligand systems are essential to understand pharmacological modulation and toxicology processes that are driven by multiple factors covering from the atomistic to cellular level. Even though the interactions taking place within a defined biophysical environment are usually intricate and complex, having a preliminary knowledge of the structural determinant and biochemical function of target enzyme in the physiological or unbound state represent a step forward in the characterization of the forces involved these processes under interaction conditions as induced by drugs. In the present work, we tackle the study of relevant binding interactions between two well-recognized betablocker drugs and the lysozyme biological target from an experimental-computational perspective. In this way, molecular docking, machine learning and perturbation analysis combined with UV–vis and fluorescence measurements will allow us to determine the allosteric regulation and functional dynamics of lysozyme by binding propranolol and acebutolol. Ó2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction The wide variety of behaviors exhibited by proteins is a direct result of the type of bonding that brings atoms together. These molecules are neither fixed to space nor bound to time, i.e., they do not keep a shape for a long time and are constantly ready and willing to change it. Several authors have described proteins as complex systems and as deformable polymers [1,2]. Other descriptions characterized by small amplitude vibrations and low frequency modes do not fully cover the delicate balance between thermal stability and fluctuations necessary for their biological functionality. The concept of fractal sponge was recently introduced, after the evaluation of 911 protein sequences through principal component analysis, where proteins are represented as objects more like sponges with fractal dimension rather than densely packed spheres [3]. In any case, all these models have common patterns, are dynamic and are based on discrete mathematics [4]. The complexity of the three-dimensional structure of proteins, with interactions of different intensity and scope, makes any movement of a residue or even atoms highly cooperative. These conformational changes can modify the catalytic activity of an enzyme rapidly, reversibly and by non-covalent interaction, this mechanism, known as allosterism, plays an important role in the dynamic molecular control of cellular metabolism [5]. Lysozyme is a very popular protein that has been used as a model on many occasions. Beyond this, lysozyme exhibit significant antibacterial properties, and is widespread in nature. It occurs in both plants and animals and is involved in defense against bacterial invasion. Due to the ability to cleave the b-(1,4)-glycosidic bond between N-acetylmuramic acid and N-acetylglucosamine of bacterial cell wall peptidoglycans, it has been reported to have anti-inflammatory, anti-cancer, and analgesic properties reports [6,7]. It has also proven its enormous versatility in other fields. For example, the porous architecture of tissue-engineered hydrogels can be obtained by modulating hydrogel degradation via various mechanisms. Thus, the incorporation of lysozyme to chitosan hydrogels accelerates the degradation rate of the crosslinked hydrogels [8]. Also, lysozyme air-filled microbubbles can be easily prepared with a high degree of control over their size and polydispersity by applying high-intensity ultrasound. Lysozyme microbubbles are highly suited as theragnostic agents because of their stable and echogenic gas core encapsulated within a thick and elastic shell [9]. Propranolol is a widely used drug for hypertension, although it is also prescribed in cases of chronic stable angina. It is a badrenergic blocking agent. It is an optically inactive compound, and only isomer 1 has adrenergic blocking activity. Acebutolol, a beta-selective beta-blocker, is also very popular in the pharmaceutical industry. Its molecular structure features an aryloxypropranolamine chain as a non-selective blocker but attached to an aromatic ring containing a polar-4 substituent. This architecture https://doi.org/10.1016/j.molliq.2022.120370 0167-7322/Ó2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ⇑ Corresponding author. E-mail address: [email protected] (J.M. Ruso). Journal of Molecular Liquids 366 (2022) 120370 Contents lists available at ScienceDirect Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq allows the molecule to self-aggregate into different structures whose result depends both on the rigidity and planarity of the aromatic ring and on the alkyl chains present [10]. The potential of this molecule lies in the fact that many b-adrenoceptor blocking agents have pharmacological effects independent of their bblocking activity, such as modification of the cell membrane. These effects are known as membrane-stabilizing activities and include nonspecific cardiac depression, myocardial conduction velocity depression, and local anesthetic activity. Like propranolol, its pharmacological effects originate in the alteration of cell membrane properties [11,12]. The synergy between experimental methods coupled with nonclassic computational structure-based strategies to explore protein-ligand allosteric interactions have undergone significant progress to date. Allowing to efficiently predict and modeling multiple molecular mechanisms with potential biomedical relevance and pretty precision. These computational methods have simultaneously demonstrated the ability to evaluate the interaction mechanisms of potential therapeutic allosteric compounds, while demonstrating accuracy, specificity, and sensitivity, to reasonably reproduce the experimental insights [13–20]. We have recently reported several studies on the complexation of drugs and proteins. These studies were performed by a combination of computational and experimental methods. Thus, Molecular Dynamics, Density Functional Theory, Machine Learning, Ligand Transport Analysis, Molecular Docking, Fast Fourier Transform, Fractal Analysis and Local Perturbation Response Maps; on the one hand, combined with Isothermal Titration Calorimetry, Differential Scanning Calorimetry, Ultrasound and Density, Zeta potential, UV–vis, Raman, Fluorescence and Small Angle X-ray Scattering, on the other, enabled us to build a detailed and contrasted picture of the complex formed detailing at the atomic level its architecture, most relevant interactions and the dynamics of the perturbation that ligand adsorption exerts on the protein structure [21–23]. In view of the foregoing herein, we present a study oriented to effectively assess the allosteric molecular mechanism of interaction of two recognized betablocker drugs (propranolol and acebutolol) with the lysozyme as target receptor by using a combination of structure-based molecular docking and spectrofluorimetric techniques. 2. Materials and methods 2.1. Computational section The computational methodology workflow was carried out had the stages that follow: i) Selection and parameterization (i.e., optimization of the ligands betablocker drugs as propranolol and acebutolol) ii) performing the perturbation analysis-based allosteric signal propagation in the whole structure of lysozyme. iii) performing the prediction of the main or best-ranked allosteric binding site in lysozyme structure and setting up the flexibility profile based on the flexibility/rigidity index and the cross-correlation matrix; (iii) Performing structure-based molecular docking with the corresponding 2D-lig-plot diagrams of interaction for the cited betablocker drugs with lysozyme as target receptor. 2.1.1. Betablocker drug selection and lysozyme allosteric binding-sites prediction In this step, the query betablocker drugs propranolol and acebutolol were retrieved from the Pubchem Data Base Chemical Structure Search as propranolol (PubChem CID: 4946; MF: 259.339 g/mol) and acebutolol (PubChem CID: 1978; MF: 336.4 g/mol) [24]. Both betablocker drugs were submitted to ligand parameterization by using a fast and accurate force field (FF) based on neural network potentials (NNPs), which are previously trained to predict drugs QM energies from a quantum-level machine learning algorithm consisting the following three steps: i) the drugs are submitted to an initial FF formed by GAFF2 parameters and AM1-BCC atomic partial charges, ii) the ligand rotable bonds (i.e., parameterizable dihedral angles) are then identified and scanned by generating a rotamers set. Then, the ligand rotamers from both betablockers are minimized with the initial FF and their reference energies are evaluated using the NNP method. iii) lastly the ligand dihedral angle parameters were fitted to reproduce the reference energies, resulting into an improved FF for both betablockers [25]. On the other hand, the prediction of the main lysozyme (PDB ID: 1HER) allosteric binding-site was identified by using DeepSite online webserver followed by an allosteric network analysis [26], which accurately ranks all the potential lysozyme cavities by excluding convex Van der Waals volumetric regions and accepting concave ones as relevant binding sites which are more likely to bind the betablocker drugs into the Van der Waals surfaces. For this instance, DeepSite tool consider all the molecular descriptors such as main protuberances and collecting information-based crystallographic descriptors on the lysozyme pocket geometry and topology, associated with high ligand binding probability or druggability properties. The cavity detection method is based on an AI- algorithm coupled with 3D-deep convolutional neural networks (DCNN) which allows the unequivocal prediction of the bestranked catalytic binding site on the lysozyme [27]. Afterward, the allosteric network procedure can calculate the allosteric coupling intensity (ACI) parameter for each residue in the previously identified lysozyme active binding site. Then, the ACI hotspots in the neighborhood of the lysozyme structure (i.e., with distance cutoff at 3.4 Å was set up) were considered as potential allosteric sites [28]. The calculation of ACI parameter always needs to be associated to the alfa parameter which dictates the ACI magnitude of the lysozyme allosteric binding site (herein by default set up at 3). Finally, the corresponding Van der Waals volumetric maps is generated and used to define the search space (i.e., the cubic grid box) to perform the computational simulations. Herein, the allosteric binding site volumetric map of the lysozyme was obtained together with the Cartesian XYZ-coordinates which offer the 3D- discrete space for setting the lysozyme docking grid box with dimensions of X = 20 Å, Y = 20 Å, Z = 20 Å and lysozyme grid box center X = 3.26 Å, Y = 24.20 Å, Z = 27.20 Å, and volume equal to 431.26 Å 3 [29]. Afterward, a crystallographic structural validation of the catalytic binding residues belonging to the lysozyme was carried performing a Ramachandran analysis [30]. 2.1.2. Performing the allosteric coupled flexibility profile of lysozyme receptor In this stage, an allosteric coupled local flexibility perturbation (allosteric-LFP) models were proposed to characterize the local flexibility properties of the lysozyme receptor in both, the unbound and bound state (i.e., in the presence of the propranolol and acebutolol). This computational algorithm is bio-inspired on the traditional elastic network models and allows to evaluate concerted movements and the directions associated with elastic normal modes of fluctuations in the allosteric binding site of the lysozyme by representing the C-( a )-atoms belonging to the residues allosteric network of communication as a Hookean potential (or ‘‘springs” model) [18]. A more in-depth description of the applied ENM algorithm will be described later in the next section. The allosteric binding site flexibility profile was set up regarding the flexibility fluctuation patterns and the flexibility/rigidity index, whose values oscillate around its equilibrium conformation and collective motions in the residues network in the unbound state M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 2 of the allosteric binding site [31,32]. For more clarity, the computational results were presented as follows: i) unbound lysozyme allosteric binding site as the control simulation experiment as reference simulation for comparison purposes only; with the ii) bound lysozyme allosteric binding site interacting with the betablockers drugs (i.e., propranolol and acebutolol). 2.1.3. Performing molecular docking on betablockers plus lysozyme interaction Then, to explore the betablockers docking mechanisms with the predicted lysozyme allosteric binding site, we applied the wellrecognized Vina scoring function developed by Trott et al. to achieve the Gibbs free energy of binging (FEB, kcal/mol) in the bound state as the sum of force-field interactions terms (i.e., Van der Waals, hydrogen bounds and electrostatic interactions) [33– 35]. In this context, the accuracy of docking was set at 50 and for this instance, the best betablocker crystallographic binding pose was selected from both drugs (propranolol and acebutolol). Here it is important to note that allosteric Gibbs free energy of binding affinity ( D G) obtained from the docking complexes were categorized as conformationally unfavorable when the D Gofthe lysozyme/betablocker complexes 0 kcal/mol, thus pointing either extremely low (or complete absence of interaction); otherwise, the lysozyme/betablocker docking complexes were categorized like thermodynamically stable from the affinity point of view. Next, the 2D lig-plot interaction diagrams were generated for the bound state to identify the different types of docking interactions of the two betablockers within the predicted lysozyme cavity, using to this end the Discovery Studio software which allows to identify and visualize the relevant non-covalent inter-molecular interactions present in the docking complexes. It automatically generating a 2D-interaction diagram that includes hydrophobic van der Waals interactions, H-bonds, electrostatics, and p - p stacking interactions [36]. 2.2. Experimental procedures 2.2.1. Reagents Propranolol (1-[Isopropylamino]-3-[1-naphthyloxy]-2-propra nol) hydrochloride (No. P-0884) and acebutolol (N-[3-acetyl-4-(2 -hydroxy-3-[isopropylamino] propoxy) phenyl] butanamide) hydrochloride (No. A-3669) were purchased from Sigma Chemical and sufficiently well characterized as pure to be used as received. Lysozyme from chicken egg white, (lyophilized power, protein 90%, 40,000 units/mg) was purchased from Sigma-Aldrich and used without further purification. Samples were freshly prepared for each experiment within 1 h. prior to usage. Solutions were made using tripe-distilled and degassed water. 2.2.2. UV–vis absorption spectra A Cary 100 Bio UV–Vis Spectrophotometer was used to measure UV–vis absorption spectra. The spectral range examined was 240– 400 nm. A Lysozyme aqueous solution of 0.286 g/L, or 0.02 mM, was utilized in the UV measurements as reference. While keeping the Lysozyme concentration constant, solutions of the two betablocker agents were added at different rates. In the case of acebutolol, concentrations ranged from 0.0 mM to 0.4 mM with increasing rates of 0.05 mM. For propranolol, the same procedure was used but, as it presented a stronger response, concentration range chosen for this compound went from 0.0 mM to 0.04 mM with an increasing rate of 0.005 mM. UV measurements were repeated twice for each drug to avoid any possible mistake or misrepresentation. No remarkable variations were found between replicates, and the data plotted is the mean value of the two results. 2.2.3. Fluorescence measurements The spectra of fluorescence emission were acquired using a Cary Eclipse spectrofluorometer. The excitation and emission splits were 5 nm. The average duration and the data interval were set at 0.5 s and 1 nm, respectively. The excitation wavelength was fixed at 280 nm and the range used for the interaction study was 240 – 450 nm. With the intent of preventing erroneous results, inner filter effects were adjusted for the quenching experiments by using the expression: F corr ¼F obs exp ðA exc þA em Þ=2  ; where F obs and F corr are the observed and corrected fluorescence intensities, respectively; and A exc and A em are the absorptions of the systems at the excitation and the emission wavelength. UV–Vis-IR Spectral Software (FluorTools) [37,38] was used for data treatment. Following a similar procedure as in the UV–vis analysis, the fluorescence spectra of Lysozyme-beta blocker complexes were measured at different concentrations. Similarly, the ranges chosen were 0.0 mM to 0.5 mM with increasing rates of 0.05 mM and 0.0 mM to 0.05 mM with an increasing rate of 0.005 mM for acebutolol and propranolol, respectively. 3. Results and discussion 3.1. Computational section 3.1.1. Perturbation analysis-based allosteric signal propagation. First, it is important to keep in mind that specific allosteric perturbations (docking interactions) on lysozyme induced by betablockers drugs (i.e., propranolol and acebutolol) can directly impact on the general biochemical function of the enzyme by affecting its ligand binding (druggability), light absorption, flexibility-based conformational properties, and induce nonphysiological post-translational modifications which could be spread or propagate over long distances across the network of residues of the lysozyme, affecting the intrinsic conformational dynamics in the whole structure of the enzyme (lysozyme). As mentioned above, the study of allosteric molecular mechanisms is essential to understand not only enzymatic activity but also cell signaling, implement pharmacological modulation strategies as well as minimize potential toxicology processes. Herein, it will be discussed in depth a computational approach to decrypt allosteric coupling under the unbound and bound state of lysozyme interacting with the betablockers propranolol and acebutolol. The allosteric-based perturbation propagation algorithm consists in a repeated stochastic process on a lysozyme whole structure network of communication residues. This method allows to extract in the first instance all the residue contacts belonging to the three-dimensional structure of the lysozyme [19,39–43].In the next step, the algorithm calculates the number of contacts between residue pairs an after, it divides the number of contacts by the atom number in each lysozyme residue. This graphtheoretic information, is turn used to obtain a perturbation propagation probability matrix according to the equations (1) and (2): P ij ¼1p ij ¼1e  a N ij ð1Þ P ji ¼1p ji ¼1e  a N ji ð2Þ where the term P ij (or P ji ) represents the perturbation probability to propagate allosteric signals from the lysozyme residue (i) to residue j and vice-versa (i.e., bi-directional allosteric coupling); the term p ij (or p ji ) represents the probability in which the perturbation of residue (i) will not induce allosteric signals to residue j and vice-versa, here Nrepresents the number of residues in the lysozyme structure and the a parameter was set to 3.0 by default to amplify the probability of occurrence of signal propagation by increasing the ACI M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 3 absolute values of lysozyme allosteric residues. Besides, an optimal value of # rounds of perturbation propagation equal to 10,000 was set by default to improve the simulations [28]. For this instance, a perturbation analysis-based allosteric signal propagation was carried out in the lysozyme whole structure by determining the allosteric coupling intensity (ACI) values in the 3D-crystallographic structure of lysozyme that represent the frequency in which several blocks or clusters of residues can be affected from the native state (unbound state) to the bound state (i.e., under betablocker interaction). See Fig. 1. Considering the obtained results, it is evident that under binding conditions both betablockers could theoretically affect the physiological signal propagation of the lysozyme in a different way when we compare the anisotropic topology of the allosteric network of the bound condition (Fig. 1,Eand H) with the lysozyme reference control (or unbound condition represented in Fig. 1,B). Particularly, it is important to note that the allosteric interaction between acebutolol and lysozyme (Fig. 1,H) exhibits a more complex pattern of interaction than its counterpart propranolol plus lysozyme (Fig. 1,E). Besides, they were detected different hotspots allosteric residues as visibly prominent hotspot-peaks (ACI values 0.90) in the three conditions evaluated as i) lysozyme plus acebutolol with ACI values (ACI_Tyr53 = 0.95, ACI_Ile58 = 1) (Fig. 1,I), ii) lysozyme plus propranolol with ACI values (ACI_- Tyr53 = 1, ACI_Ile58 = 0.90, ACI_Leu83 = 0.93) (Fig. 1,F) and lysozyme unbound state with ACI hotspot-peaks values (ACI_Tyr53 = 1, ACI_Leu83 = 0.93). Here, it is important to note that the allosteric residue Tyr53 in all the cases exhibits ACI values 0.9, strongly suggesting that it plays a critical role in the allosteric signal propagation in both, the native and under interaction conditions. In this regard, it is well-known that the tyrosine residues actively participate in phosphorylation/dephosphorylation mechanisms which are directly involved in allosteric mechanisms and posttranslational modifications for most enzymes (e.g., lysozyme). 3.1.2. Allosteric coupled flexibility analysis and structure-based molecular docking To deepen on the allosteric binding interactions induced by the betablocker drugs interacting with the lysozyme, an allosteric coupled local flexibility perturbation (allosteric-LFP) approach was performed. This theoretical procedure allows to distinguish between allosteric and binding active sites of the lysozyme [28]. The allosteric-FLP approach provides a biophysical and structural view on how the allosteric inter-residue network of the lysozyme Fig. 1. Results of the allosteric analysis of lysozyme under unbound panels (A-C) and bound state represented by the panels (D-F) and (G-I) as lysozyme plus propranolol and lysozyme plus acebutolol; respectively. In the top, panels A),D), and G) represent the 3D-crystallograhic structure of lysozyme rendered according to the allosteric coupling intensity (ACI) parameter. Here, labeled-blue regions in the structure depict lysozyme residues with low-allosteric influence while labeled-yellow tored regions represent lysozyme residues with high-allosteric influence, in the whole structure. In the middle, panels B),E), and H) represent the anisotropic network topology of allosteric signal propagation composed by critical allosteric residues and evaluated regarding their relative importance in the allosteric pathways from low importance (0) to high importance (1) in the color intensity bar. Herein, the starting and ending point of the anisotropic signal propagation is represented by the labeled-dark red circles in the allosteric network of the lysozyme for the unbound panel (B) and bound state panels (Eand H; propranolol and acebutolol, respectively). At the bottom, bar chart representation of the ACI values (y-axis) for each lysozyme residue position index (x-axis) for the unbound panel (C) and bound state panels (Fand I; propranolol and acebutolol, respectively). Herein, the most prominent peaks are labeled-red and represent ACI hotspots denoting the maximum influence in the allosteric signal propagation in the lysozyme structure (e.g., allosteric regulatory residue Tyr53). M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 4 could be theoretically affected from its local flexibility properties under the unbound and bound state. In fact, the collective conformational dynamics of the lysozyme is strictly dependent on the flexibility/rigidity index (FRI) parameter that means the qualitative and quantitative balance of the flexibility/rigidity properties within a given binding site (e.g., allosteric binding site) evaluated. Fig. 2 and Fig. 3. The allosteric maps performed for the lysozyme bound state (Fig. 2 C,Dand Fig. 3 C,D), allow the visual representation of the strength of local fluctuation-based on perturbations between backbone–backbone atom contacts belonging to the two sequenceadjacent residues of the allosteric site by generating the matrix of local perturbation signals from the binding states and binding conditions (unbound and bound state of lysozyme) [43]. In this regard, we can observe that both betablockers theoretically affect the flexibility properties in the allosteric binding site when compared the allosteric LFP of the bond states in the presence of propranolol and acebutolol (Fig. 2Dand Fig. 3D; respectively) with their corresponding theoretical or reference control in the lysozyme unbound state (Fig. 2Cand Fig. 3C). Herein, we can observe significant flexibility perturbations in large blocks or coupled sensors (i) and effectors (j) residue clusters, which in turn involve large movements or allosteric fluctuations in a correlated way and with an accentuated propagation of allosteric signals over great distances across the lysozyme sensors (i) and effectors (j) inter-residue network. Please, note the concerted and correlated allosteric signal fluctuations in the 2D-matrixes represented by the labelled-red regions spreading in the allosteric LFP maps of the lysozyme plus propranolol (Fig. 2D) and lysozyme plus acebutolol (Fig. 3D) compared with the LFP map of the reference state (i.e., native lysozyme). Following this idea, it is important to note the existence of subtle differences between the allosteric LFP maps of both propranolol and acebutolol. Particularly, the cited differences appear associated to a different pattern of allosteric signalbased local flexibility perturbations that cover the region of sensor-effector residues (i, j) in the interval from the lysozyme regulatory residue SER40 to the binding target residue Ile98. The results obtained from the LFP maps fit excellent with the previous analysis on the allosteric coupling intensity (ACI) hotspots (refer to Fig. 1, panels C, F and I) because they represent the cluster of allosteric residues with the maximum influence in the allosteric signal propagation in the lysozyme structure (e.g., allosteric regulatory residue Tyr53) with the with the highest ACI value in the three binding conditions evaluated. Likewise, we can perform the allosteric analysis from a rigidity perspective, due to the allosteric conformational dynamics of the lysozyme also depend on the rigidity properties of the lysozyme target (lysozyme catalytic binding site) and how it could be affected under allosteric interactions with the betablocker drugs (i.e., propranolol and acebutolol) in the vicinity of the cited active site. In this regard, we can apply the continuum elasticity approach with atomic rigidity which represents a robust multiscale formalism for describing the main transition conformational states based only in a rigidity assessment of the LFP maps [32,44]. In this sense, the rigidity index can be estimated by considering an N-atom representation of the lysozyme structure under Fig. 2. On the left, 3D-Van der Waals surface representation of the main lysozyme allosteric binding site rendered according to the flexibility/rigidity index (FRI) obtained by applying an allosteric local flexibility perturbation analysis (LFP model) for both states A) unbound lysozyme (i.e., unperturbed allosteric lysozyme binding site in the native state) and B) bound lysozyme (or locally perturbed allosteric binding site by docking interaction with propranolol). Herein, the color intensity bar depicts the flexibility properties-based FRI index for the cited states from blue (conformational rigidity properties) to red (conformational flexibility properties). In the middle, 2D-matrix representation of allosteric local flexibility perturbation map (i.e., LFP maps based i,j-residue motion correlation C ij ) showing the local flexibility properties between sensor (i) and effector (j) residues pairs (i,j) for the allosteric binding site for both states as C) unperturbed allosteric lysozyme binding site in the native state and D) perturbed allosteric lysozyme binding site state induced by the interaction with propranolol. Herein, the pairs subject to fully correlated allosteric motions (0.2 C ij 1) are labeled-red indicating the same direction for the flexibility of lysozyme residues i and j, while those with anticorrelated allosteric motions (i.e.,C ij 0) are labeled-blue (i.e., opposite direction for the allosteric motions of lysozyme residues i and j). The moderately correlated and uncorrelated (C ij 0) map regions are labeled light red and blue, respectively. On the right, E) 2D-Lig-plot interaction diagram of the lysozyme plus propranolol docking complex. M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 5 unbound and bound state. Being the coordinates of lysozyme active site residues (or C( a )-atoms) as rigid bodies. This consideration has a paramount relevance from the biochemical point of view because it is well-known that the active binding site of a given enzyme (lysozyme) should necessarily maintain its complementary rigidity attributes (i.e., to maintain the catalytic residue orientation) from the three-dimensional point of view given by [32,44]. This fact allows the correct adjustment of physiological substrates (or ligands as propranolol and acebutolol). In this regard, we have the ||d i d j || that represents the Euclidean space distance between a given lysozyme i-th sensor catalytic residue and the corresponding effector j-th residue. Herein, a general correlation kernel, U (||d d j ||; g j ), represent a real-valued monotonically based on a decreasing function following conditions (Eqs. (7) and (8)). Ukdd j k; g j  ¼1likekdd j k!0ð3Þ Ukdd j k; g j  ¼0likekdd j k!1 ð4Þ where, g j is to denotes an atomic scale parameter. While the correlation transition allosteric matrix (C ij ) between the i-th and j-th lysozyme catalytic residues (C( a )-atoms) is given by the eq. (9). C ij ¼Ukdd j k; g j  ð5Þ the C ij represent a correlation transition matrix (refer to Fig. 2 C, D, and Fig. 3 Cand Dabove) which contain relevant information regarding rigidity properties as labelled-blue regions in the allosteric LFP maps. Then, we can define a position-dependent rigidity function RdðÞas eq. (6): RdðÞ¼ X N j¼1 a j Ukdd j k; g j  ð6Þ herein, a j represents an atom type dependent weight as: C, N, S, and O atoms which have different weights. Next, to model the position (d) dependent rigidity function RdðÞwe can apply two generalized equations as an exponential function (eq. (7)) and the Lorentz functions (eq. (8)); to explore the lysozyme rigidity properties in the context of local flexibility perturbations maps. RdðÞ¼ X N j¼1 Ukd i d j k; g ij  ¼e  kdidjk g j  k ;k>0ð7Þ RdðÞ¼ X N j¼1 Ukd i d j k; g ij  ¼1 1þ kd i d j k g j  v ; v >0ð8Þ It is expected that the correlation between any two allosteric target residues of the lysozyme should decay according to their distance when these are close to rigid regions which is expressed in the allosteric LFP maps anticorrelated allosteric motions (i.e.,C ij - Fig. 3. On the left, 3D-Van der Waals surface representation of the main lysozyme allosteric binding site rendered according to the flexibility/rigidity index (FRI) obtained by applying an allosteric local flexibility perturbation analysis (LFP model) for both states A) unbound lysozyme (i.e., unperturbed allosteric lysozyme binding site in the native state) and B) bound lysozyme (or locally perturbed allosteric binding site by docking interaction with acebutolol). Herein, the color intensity bar depicts the flexibility properties-based FRI index for the cited states from blue (conformational rigidity properties) to red (conformational flexibility properties). In the middle, 2D-matrix representation of allosteric local flexibility perturbation map (i.e., LFP maps based i,j-residue motion correlation C ij ) showing the local flexibility properties between sensor (i) and effector (j) residues pairs (i,j) for the allosteric binding site for both states as C) unperturbed allosteric lysozyme binding site in the native state and D) perturbed allosteric lysozyme binding site state induced by the interaction with acebutolol. Herein, the pairs subject to fully correlated allosteric motions (0.2 C ij 1) are labeled-red indicating the same direction for the flexibility of lysozyme residues i and j, while those with anticorrelated allosteric motions (i.e.,C ij 0) are labeled-blue (i.e., opposite direction for the allosteric motions of lysozyme residues i and j). The moderately correlated and uncorrelated (C ij 0) map regions are labeled light red and blue, respectively. On the right, E) 2D-Lig-plot interaction diagram of the lysozyme plus propranolol docking complex. M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 6 0). They are labeled-blue (i.e., opposite direction for the allosteric motions of lysozyme residues i and j) as observed in the bound state simulation conditions for lysozyme plus propranolol, and lysozyme plus acebutolol. On the other hand, it is important to highlight that under propranolol and acebutolol interaction (i.e., bound, or perturbed lysozyme state) a more complete analysis on the allosteric-LFP should necessarily consider a quantitative thermodynamics description of the causal effect of allosteric ligand binding by including the determination of the allosteric Gibbs free energy of binding affinity ( D G). See Fig. 2Eand Fig. 3E. Following the theoretical results, we strongly suggest that, in general terms, both betablockers can allosterically interact with the lysozyme allosteric binding site (denoted by F) following a spontaneous thermodynamic process. Particularly, we identify a prevalence of hydrophobic binding events based on a significant number and diversity of aromatic p - p stacking allosteric interactions (eg., p -sigma, p -alkyl, p -donor) in both betablocker drugs. In this sense, in the case of propranolol, its relative allosteric Gibbs free energy of binding affinity D G (or D G F )=6.5 kcal/mol) slightly greater respect to the acebutolol with allosteric affinity of D G (or D G F ) = –5.7 kcal/mol theoretically suggests that in the propranolol this fact could be directly attributed to the presence of an additional thermodynamic contributions from benzene-ring belonging to naphthalene ring which is absent in the case of acebutolol. Herein, the aim is to evaluate the quantitative causal effect of the betablockers given by the allosteric free energy gain/loss as D G(L E ?BL E ) upon the transition binding states between the unbound lysozyme (indicated as L E ) and the betablockers denoted as B (acebutolol or propranolol) in the docking complex. This allosteric D G(L E ?BL E ) represent a scoring function that include Vina interaction terms as hydrophobic Van der Waals interactions, H- bonds, electrostatics, and p - p stacking interactions [41], and consider the summatory of independent energetic contributions of sensor (i) and effector (j) residues belonging to the predicted allosteric site denoted as F in the lysozyme (L E ) according to the equation (9): D G F L E !BL E ðÞ¼ 1 n F X i2F D G i L E !BL E ðÞ ð9Þ where the term D G i L E !BL E ðÞ represent free energy gain/loss release per allosteric binding residues from L E , and n F represents the total number of residues in the allosteric site F. Besides, considering that the signal propagation in the whole lysozyme structure was evaluated through the strength of the allosteric effects in the whole lysozyme structure should be equally considered as in the equation (10): D G L E L E !BL E ðÞ¼ 1 n L E X i2L E D G i L E !BL E ðÞ ð10Þ herein, n L E represent the total number of lysozyme residues (i.e., catalytic, and allosteric residues). Therefore, the terms D G i L E !BL E ðÞ , D G F L E !BL E ðÞ and D G L E L E !BL E ðÞ allows to evaluate the causality and the energetic free energy contribution at different levels as i) per residue, ii) per allosteric site in the lysozyme, and the whole lysozyme structure; respectively. Then, we can also model the allosteric free energy for the case of multiple allosteric sites formed by large residue clusters (c i ) in the lysozyme following the equation (11): D G c i B n1ðÞ !B n L E  ¼ D G c i L E !B n L E ðÞ D G c i L E !B n1 L E  ð11Þ Which provides information on the thermodynamics contribution under sequential binding of allosteric ligand modulators (i.e., propranolol and acebutolol). Here the term B n L E indicates the docking interaction system with n ligand molecules (propranolol and/or acebutolol) bound to the lysozyme L E . The same reasoning could be efficiently adapted in the previous equations (9) and (10) as represented by the new ones as equations (12) and (13): D G F B n1ðÞ !B n L E  ¼1 n F X i2F D G i B n1 L E !B n L E ð12Þ D G L E B n1ðÞ !B n L E  ¼1 n P X i2L E D G i B n1 L E !B n L E  ð13Þ On the other hand, it is important to note that despite of the close values of allosteric free energy obtained for the D G F (propranolol) ffi D G F (acebutolol) different types of allosteric interactions with noticeable differences according to the crystallographic orientation of these ligand within the allosteric site were found (Please refer to the 2D-lig-plot diagrams Fig. 2Eand Fig. 3E). In this regard, in the case of the propranolol we identify a prevalence of three types of hydrophobic based p -interactions involving the naphthalene aromatic moiety with allosteric lysozyme residues as: i) mixed p -donor hydrogen bond plus p -alkyl interactions with the Ala107 which has the highest influence in the stabilization of the propranolol/lysozyme docking complex, following by ii) an interaction based p - p -T-shaped between the naphthalene aromatic moiety with the TRP108 where the edge-to-face geometry orientation as a sidewise electron cloud of 1 of the ring belonging to the naphthalene aromatic moiety and the head on electron cloud of the ring from the TRP108 was identify, iii) an interaction-based p -alkyl between the naphthalene moiety of propranolol and the target residue Ile98 that also plays a relevant role in the stabilization of the docking complex by involving process of charge transfer which help to an efficient drug accommodation in the allosteric binding site of the lysozyme receptor. Besides, we have identified iv) an isolated p -alkyl interaction of the 1-(isopropylamino) with the lysozyme residue Val109. In the same sense, the p -donor interaction of the propranolol naphthalene moiety contributes to the docking complex stabilization by involving p -donor hydrogen interaction with the allosteric residue Trp63. Herein, we suggest that the unfavorable donor-donor interaction identified for the 1- (isopropylamino) moiety interacting with the lysozyme residue Asn59 does not appear to significantly affect the strength of the interaction in the allosteric site. (Refer to Fig. 2E). On the contrary, the 2D-lig-plot interaction diagram for the acebutolol plus lysozyme revealed a different pattern. In this case, a high prevalence of stabilizing hydrogen bond interactions was identified which involve several groups of the 2-acetyl-4- (butanoylamino)-phenyl ether of acebutolol with the lysozyme allosteric residues Asp52, Gln57 with classic (or conventional) type of hydrogen bonds. While the aromatic benzene ring of acebutolol was able to strongly interact with the Ala107 through multiple interactions as i) carbon-hydrogen bond interaction, ii) p -donor interaction and iii) a single p -alkyl interaction combining for this instance hydrophobic forces plus charge transfer processes between the ligand and the allosteric site. Lastly, relevant interactions between the 5 0 -butyramido-2 0 -(2-hydroxy-3-isopropylamino propoxy) acetophenone fragment of acebutolol with the neighboring allosteric residues Trp62 and Trp63 were identified. So, two stabilizing interactions triangular-shaped formed by p -sigma interaction plus a hybrid p -donor and carbon-hydrogen bond interactions specifically with the Trp62, and a single but not less important with the p -alkyl interaction with the allosteric residue Trp63. (Refer to Fig. 2E). Finally, a crystallographic structural validation of the allosteric binding residues belonging to the lysozyme was carried by performing a Ramachandran analysis. Our work incorporates the Ramachandran crystallographic validation as provided in the Supplementary Materials(Figure SM1)[30], in order to highlight the importance only for the allosteric residues M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 7 involved in the interactions with the betablocker drugs (i.e., propranolol and acebutolol). Then, taking into account that allosteric effects are generally triggered by local perturbations in the protein (i.e., from lysozyme allosteric binding residues) we do not consider the Ramachandran crystallographic validation in the whole lysozyme structure to improve the quality of our theoretical observations. However, it this important to note that the same Ramachandran principles were followed being an important in silico validation procedure, that corresponds to a 2D-projection on the plane from the 3D-crystallographic structure of lysozyme (. pdb model) where all the possible allosteric residue conformations (flexibilities) are defined into the plot based on the torsion dihedral angles /vs. w around the lysozyme peptide-bonds. In this methodology, the allowed torsion values of the dihedral angles were found within the Ramachandran colored light-green contour for individual allosteric residues being considered as conformational-favored. As such, theoretical observations allow to avoid the presence of false positives docking results from the structural and conformational point of view. 3.1.3. Experimental validation through spectroscopic evaluation. To validate the computational results, a well-proven and generally utilized method for assessing conformational changes in protein-ligand complexes is spectroscopic analysis [45].By inspecting the spectral curves under various circumstances and concentrations, it is possible to determine the development of a complex or the structural alterations that may result from the organization of new structures. In this case, lysozyme absorbs in two primary bands. The p - p * electronic transitions of the peptide backbone C=O cause the stronger one from 200 to 230 nm and it represents the lysozyme framework conformation. The weaker one from 260 to 290 nm is caused by the absorption of aromatic amino acids and shows the alteration of the chromophore microenvironment [46]. Tyrosine, phenylalanine, and tryptophan are three aromatic amino acid residues that are specifically present in lysozyme and are the primary causes of this absorption peak. Tryptophan is the one of the three stated whose contribution is strongest, resulting in a peak at 278 nm [47]. Fig. 4 shows how the absorption spectra of lysozyme get affected by the presence of drugs. At 298 K, the peak at 278 nm changed in a similar way for both beta-blocker agents. Upon complexation, the 278 nm band of lysozyme slightly increase its maximum value, suggesting a decrease in hydrophobicity and an increase in polarity [48]. This fact reinforces the idea of existing multiple interactions between the aromatic benzene ring of acebutolol with lysozyme through p relations, combining hydrophobic forces with charge transfer. In the case of propranolol, the hyperchromic effect also justify the hydrophobic based p -interactions between the naphthalene aromatic moiety of the propanolol and the allosteric lysozyme residues previously predicted by molecular docking. Apart from the UV–vis technique, one of the most sensitive and trustworthy techniques for examining supramolecular host-guest interactions is fluorescence spectroscopy, which uses the guest’s noticeable change in fluorescence spectra to both qualitatively and quantitatively evaluate the interaction. As shown in Fig. 5,the presence of acebutolol promotes the quenching of lysozyme fluorescence and the red shift in the highest emission wavelength suggests again a change in the Trp and Tyr residues’ surroundings, suiting the solvent to a more polar, less hydrophobic, fluorophore environment. For the propranolol, the same bathochromic effect is observed, but, in this case, the fluorescence intensity, instead of diminishing, it increased regularly on interaction. This enhancing fluorescence tendency was already probed for this drug when combined with other blood plasma proteins, like human serum albumin [49,50]. The increase in the absorption and in fluorescence as well, implies that binding propranolol to lysozyme caused microenvironmental changes in the protein and suggest the formation of lysozyme–propranolol complexes. For its part, the presence of an isosbestic point at 318 nm indicates the existence of an equilibrium between energy transfer from lysozyme to the ligands is highly probable for free and bond propranolol [51]. To provide more information about the molecular interactions present in the studied biochemical systems, the fluorescence quenching was investigated. It is the process through which the intensity of a fluorophore’s fluorescence decreases, and it can be caused due to a variety of chemical interactions with quencher molecules. There are three basic types of quenching processes that can occur when a fluorophore is paired with a ligand: collisional, static, and mixed processes. Collisional quenching refers to the neutralization of an excited state fluorophore by a quencher. Instead of relying on diffusion and molecular collision, static quenching requires the growth of a ground state non-fluorescent complexation. Finally, mixed quenching is brought on by both col- Fig. 4. UV absorption spectra of lysozyme in the absence and presence of a) Acebutolol and c) Propranolol. Insets: b) From yellow to blue: UV absorption spectra of increasing concentrations of Acebutolol, from black to grey: UV absorption spectra of lysozyme (0.02 mM) and increasing concentrations of Cloxacillin combined. d) From yellow to blue: UV absorption spectra of increasing concentrations of Propranolol, from black to grey: UV absorption spectra of lysozyme and increasing concentrations of Propranolol combined. C lysozyme = 0.02 mM. M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 8 lision and complex formation with the same quencher [52]. Typically, to dilucidated the corresponding mechanism, the Stern- Volmer equation is used [53]: F 0 =F¼1þK sv ½Q¼1þk q s 0 ½Qð14Þ where F 0 and F are the fluorescence intensities in the absence and presence of the quencher; K sv ,k q , and s 0 are the Stern–Volmer quenching constant, the biomolecular quenching constant, the excited state lifetime of the biomolecule in the absence of the quencher ( s 0 = 5.9 10 -9 )[54], respectively and [Q], the concentration of the quencher. Figure SM2 a) depicts the results from the Stern Volmer plot for acebutolol. The slope of the plot, K SV , can be used to determine the quenching constant. In this case K SV = (3.92 ± 0.13) 10 3 L mol 1 , which corresponds to a k q = (6.64 ± 0.03) 10 11 M 1 s 1 .Itis important to notice that even at the higher values of the concentration of acebutolol, the Stern-Volmer plot does not display any deviation from linearity toward the y-axis, which suggests that there is no coexistence between static and mixed quenching. Specifically, the type of quenching may be inferred from the analysis of this condition and the measured values of k q . This analysis is typically regarded as a good first observation for identifying the mechanism in the interaction between proteins and ligands [55]. When the interaction is mostly controlled by diffusion, the values of k q for dynamic quenching are in the range of 1 10 10 M 1 s 1 . The quenching type is static for values greater than the diffusioncontrolled limit. For acebutolol and lysozyme the quenching rate constant at 25 °C is in the order of 10 times the maximum diffusion rate limit, which is a strong indication that the quenching involved is a static one. To obtain the binding parameters (i.e. the number of binding sites (n) and the association constant (k A )) a non-fluorescent complex was considered between the protein and drugs. So, for the binding interactions, an equilibrium should be assumed which is given by the expression [56]: PþnDPD n ð15Þ where nis the number of binding sites, Pand Dare the protein and drug concentrations. So, the equilibrium constant can be obtained from: k A ¼PD n ½ P½D½ n ð16Þ PD n ½is the equilibrium concentration. Considering that the complex and drug are non-fluorescent, it can be assumed: P½ 0 ¼kF 0 ;ð17Þ P½¼kFð18Þ D½¼D½ 0 nPD n ½ ð19Þ P½ 0 is the total concentration of protein (0.02 mM) and D½ 0 the concentration of the drug. Then: k A ¼P½ 0 P½ P½fD½ 0 nðP½ 0 P½ ð20Þ log P½ 0 -P½ P½ ¼logk A þnlog D½ 0 nðP½ 0 P½Þ  ð21Þ logF 0 -F F¼logk A þnlog D½ 0 nP½ 0 ðF 0 FÞ F 0  ð22Þ where F 0 and F are the same as in Eq. (14). From Eq. (22), the equilibrium constant (k A ) and the number of binding sites (n) can be obtained. Data from the fitting for acebutolol is presented in Figure SM2 b). In the case of propranolol, as it produced enhanced fluorescence, the binding constant was calculated through the Benesi-Hildebrand equation, assuming a 1:1 stoichiometric interaction [57 58]: 1 FF 0 ðÞ ¼1 F 1 F 0 ðÞ þ1 F 1 F 0 ðÞK eq D½ ð23Þ F 1 corresponds to the maximum possible fluorescence intensity when all molecules are complexed. Figure SM2 c) represents a plot of 1/(F – F 0 ) as a function of 1/[propranolol], and the binding constant value K is determined using the intercept and value is shown in Table 1. Fig. 5. Fluorescence spectra of lysozyme in the absence and presence of a) Acebutolol and b) Propranolol at k ex = 280 nm, T = 298 K. Black lines correspond to lysozyme alone. C lysozyme = 0.02 mM; C Aceb (x 10 -2 ) (i-xi) = (0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50) mM; C Propr (10 -3 ) (I-XI) = (0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50) mM. Table 1 Constants and binding parameters for the interaction of lysozyme with Acebutolol and Propranolol at 298 K. Binding parameters n10 -2 K A (mol 1 )R 2 Acebutolol 0.67 ± 0.01 2.58 ± 0.04 0.986 Propranolol 1.00 ± 0.00 56.66 ± 0.03 0.981 M. González-Durruthy, Ramón Rial, Z. Liu et al. Journal of Molecular Liquids 366 (2022) 120370 9