key: cord-0793508-6iwtr6tv
authors: Alvarado, Ysaias José; Olivarez, Yosmari; Lossada, Carla; Vera-Villalobos, Joan; Paz, José Luis; Vera, Eddy; Loroño, Marcos; Vivas, Alejandro; Torres, Fernando Javier; Jeffreys, Laura N.; Hurtado-León, María Laura; González-Paz, Lenin
title: Interaction of the New Inhibitor Paxlovid (PF-07321332) and Ivermectin With the Monomer of the Main Protease SARS-CoV-2: A Volumetric Study Based on Molecular Dynamics, Elastic Networks, Classical Thermodynamics and SPT
date: 2022-05-14
journal: Comput Biol Chem
DOI: 10.1016/j.compbiolchem.2022.107692
sha: 3fd16f6a717729797cdd48c2fa1fb39397177775
doc_id: 793508
cord_uid: 6iwtr6tv

The COVID-19 pandemic has accelerated the study of drugs, most notably ivermectin and more recently Paxlovid (PF-07321332) which is in phase III clinical trials with experimental data showing covalent binding to the viral protease M(pro). Theoretical developments of catalytic site-directed docking support thermodynamically feasible non-covalent binding to M(pro). Here we show that Paxlovid binds non-covalently at regions other than the catalytic sites with energies stronger than reported and at the same binding site as the ivermectin B1a homologue, all through theoretical methodologies, including blind docking. We volumetrically characterize the non-covalent interaction of the ivermectin homologues (avermectins B1a and B1b) and Paxlovid with the mM(pro) monomer, through molecular dynamics and scaled particle theory (SPT). Using the fluctuation-dissipation theorem (FDT), we estimated the electric dipole moment fluctuations at the surface of each of complex involved in this study, with similar trends to that observed in the interaction volume. Using fluctuations of the intrinsic volume and the number of flexible fragments of proteins using anisotropic and Gaussian elastic networks (ANM+GNM) suggests the complexes with ivermectin are more dynamic and flexible than the unbound monomer. In contrast, the binding of Paxlovid to mM(pro) shows that the mM(pro)-PF complex is the least structurally dynamic of all the species measured in this investigation. The results support a differential molecular mechanism of the ivermectin and PF homologues in the mM(pro) monomer. Finally, the results showed that Paxlovid despite beingbound in different sites through covalent or non-covalent forms behaves similarly in terms of its structural flexibility and volumetric behaviour.

During the COVID-19 pandemic researchers sought new therapeutic approaches to target SARS-CoV-2. A common approach was the repurposing of approved pharmaceutical drugs to reduce costs, risks and time to approval/adoption in clinical settings. Computational methods with a biophysical approach can be used to discover different interactions and relevant drug phenomena that are not considered and are not found during the clinical trial process, especially for existing drugs that could offer a broader medical scope such as ivermectin [1, 2] .

Ivermectin is composed of an approximately 80:20 mixtures of two homologues, avermectin B1a (B1a) and avermectin B1b (B1b), which differ in the presence of a secbutyl and an isopropyl group, at the C25 position, respectively (see Fig.1 ). A wide variety of computational biophysical studies have been carried out to determine the capacity of ivermectin to treat COVID-19 due to shown multitarget antiviral activity against SARS-CoV-2 in vitro [3, 4, 5, 6] . Ivermectin has shown a higher binding energy to different target proteins of SARS-CoV-2 compared to other antiviral drugs and new J o u r n a l P r e -p r o o f 4 candidates, such the alphaketoamide type inhibitor 13b [2] . Many proteins have been suggested to interact with ivermectin, such as the M pro protease which is responsible for most of the post-translational modifications of SARS-CoV-2 polypeptides [1] [2] [3] [4] .

Despite all the efforts made to demonstrate the structural and energetic changes induced in M pro by the binding of ivermectin [7] , the thermodynamic changes that accompany this process from a volumetric perspective still remain unknown. As some authors [8] [9] [10] have conclusively demonstrated, this type of study provides information on the role played by hydration in modulating protein stability and its role in ligand binding reactions. It is important to mention that it has been shown that protein-protein and protein-ligand interactions occur with the loss or gain of water molecules around the ligand and/or protein (biological water), a phenomenon with potential pharmacological impact because it is key to the biological activity of these proteins [11] [12] [13] [14] [15] [16] [17] [18] [19] . Despite the theoretical effort made and the in vitro studies of its antiviral activity to date, there are still no experimental reports that confirm the binding and perturbation of M pro by ivermectin. For this reason, we carried out a comparative study of these homologues of ivermectin with the new clinical phase III drug called Paxlovid or PF-07321332 (from here on we will call PF-07321332) of which there is theoretical and experimental data that show its binding to protease M pro in catalytic site [20] [21] [22] [23] [24] . Theoretically using directed Docking, a model has been proposed that this binding occurs through a noncovalent binding of PF at the catalytic site causing the formation of the covalently bound species at the site through the formation of a thiolate-imidazolium group [20] .

Interestingly this species is less thermodynamically favored than the non-covalent complex [24] . This study, in addition to being interesting, opens the possibility of evaluating by blind docking the possible formation of complexes with PF-07321332 by J o u r n a l P r e -p r o o f 5 thermodynamically favorable non-covalent bonding at sites other than the catalytic site that has been reported [25] .

In this study we focus on analyzing the volumetric changes of the structures using the minimum binding energy resulting from blind dockings. Blind docking represents one of the greatest challenges in the field of molecular docking [26] . Blind docking has been described as an unbiased molecular docking approach as it scans the protein structure to locate the ideal ligand binding site, representing a probabilistic approach.

However, blind docking represents one of the greatest challenges in the field of molecular docking [26] . In addition, it has been suggested for the study of potential therapeutic agents against coronaviruses [27] .

Despite the limitations of blind docking several tools have been produced such as DockThor which has been been optimized for ligand docking to SARS-CoV-2 targets [31] and has been shown to be relatively successful in situations where the coupling box spans the entire surface of the receiver [32] . Especially since, in a virtual screening circumstance, this method has a great advantage: it can determine the best binding pocket for each candidate ligand in a single docking run [32] .

It has been suggested that although some potent antivirals can bind strongly to the active site, exploring other regions of the protein with the blind docking approach should not be ruled out because a preference for binding to target sites has been demonstrated [33] . Especially since promising non-covalently bound inhibitors have been reported by applying the molecular docking strategies presented here [34] , and that both covalent and non-covalent binding free energy contributions are important in the binding affinity of an inhibitor towards its target [35] .

J o u r n a l P r e -p r o o f 6 In addition, it has been suggested that in M pro allosteric inhibition could take place due to the presence of an allosteric site in the protein structure that is not found in the active site. In fact, it is noteworthy that residues such as P122 have been implicated as keys to interactions at the dimer interface and contribute to M pro dimerization [28] .

Interestingly, these residues are part of those involved in the interactions with the homologue B1b and PF-07321332 (non-covalent bonding) described in our study (see supplementary material - Table 1S ). In light of these observations, we decided to apply the blind docking approach, so as not to exclude the possible inhibition strategies suggested against M pro , such as: (1) interaction with active/catalytic site regions, (2) interaction with allosteric site regions, and (3) interaction with interphase regions associated with enzyme dimer formation [28] , the latter strategy was the most relevant to our investigation, especially, since M pro is known to depend on homodimerization for its biological activity [29, 30] .

The volumetric results obtained provide further insight into the possible mechanism by which these drugs inhibit the reaction of formation of the homodimeric protease M pro . For this we build protein-ligand complexes (M pro + avermectin) and (M pro + PF-07321332) on which we have used the Voronoi model implemented in the 3Vee program [36, 37] , HullRad program [38] , the volumetric models proposed by Chalikian et. al [14, 15, 39] and Graziano-Lee models based in scaled particle theory (SPT) [40] [41] [42] [43] .

This approach can be translated to other protein-ligand complexes to more accurately determine interactions important for catalysis and the biological activity of these proteins.

In this manuscript we have, for the first time, reported the computational biophysics study of changes induced in the volumetric properties and the hydration of J o u r n a l P r e -p r o o f 7 the monomer of the M pro protease (mM pro ) and the respective non-covalent complexes formed by binding ivermectin homologues and PF-07321332. In the last section of this paper, a comparative study is also carried out between the non-covalent and covalent complexes thermodynamically more feasible between PF and mM pro detected using blind docking and molecular dynamics. The differential chemical group of each homologues is indicated (for B1a it is secbutyl and B1b is an isopropyl, respectively).

The crystal structure of the main protease of SARS-CoV-2 (PDB: 6LU7) was used as it is a key enzyme of coronaviruses and has a fundamental role in mediating viral replication and translation, making it an attractive target for drugs [44, 45] . All structures were obtained in PDB format from the RCSB Protein Data Bank J o u r n a l P r e -p r o o f 8 (https://www.rcsb.org/). Since the X-ray structure of the PF-07321332-M pro complex is currently available in the Protein Data Bank, and since a covalent bond has already been described, it was considered to study the non-covalent bond as it is a very likely thermodynamic bond. A comparative analysis was performed applying two sampling algorithms (DockThor and AutoDock) for the prediction of compound affinity by detection of covalent and non-covalent binding as suggested for multiple targets including those associated with SARS-CoV-2 [34, 46] . The reported covalent docking was predicted from the AutoDock algorithm's covalent docking approach, which has been validated with various biological systems [47, 48] . For this, the enzyme was loaded in the PDB format and the compound of interest was loaded separately as a mol2 file.

Subsequently, the binding center was assigned to the frame corresponding to residue

Cys145. In addition, the docking calculations were performed following the recommendations of the sampling region of the tests experimentally [24] . Predictions of the relative covalent and non-covalent binding energies were made from the sampling algorithm and scoring function offered by AutoDock, and using Molegro Molecular predictions, in order to calculate the contributions of the molecular interactions offered by the tool, such as hydrogen bonds, hydrophobic bonds, electrostatic bonds, and also covalent bonds at the binding site as suggested [49, 50] (see Supplementary material - Table 1S ). In addition, the non-covalent union was compared with the covalent one after the construction of the covalent model as suggested [20] . The X-ray frame of an alpha-ketoamide inhibitor covalently linked to M pro and strictly related (PDB: 7BPR) was used as suggested for the construction of the complex and the comparison between covalent and non-covalent binding of PF-07321332 [20] . .2, Fig.3, Fig.4, Fig.5 ). The DockThor-VS platform utilizes a topology file for the ligand and a specific input file for the protein containing the atom types and partial charges from the MMFF94S force field, and both are generated using the in-house tools MMFF Ligand and PdbThorBox. The file of the ligand is generated by the MMFF Ligand program, which utilizes the facilities of the Open Babel chemical toolbox for deriving partial charges and atom types with the MMFF94S force field, defining the rotatable bonds and the terminal hydroxyl groups, and calculating the properties necessary for computing the intramolecular interactions. In the DockThor program, both protein and ligand are treated with the same force field in the docking experiment. The complexes were built using the flexibility algorithm and blind docking.

The affinity prediction and ranking of distinct ligands are performed with the linear model and DockTScore GenLin scoring function. To increase accuracy 30 runs were made with 10 6 evaluations per run. As is usually done, all the water molecules were removed and the PDB files were separated into two different files, one containing the protein and the other containing the ligand structure. All the molecular force field parameterisations were performed automatically by the programs cited. The remaining settings, conditions and parameters offered by the program were used in the default mode [31] .

For the prediction of the covalent binding complex, the enzyme predetermined by DockThor was used for the study of M pro and the compound of interest was loaded separately as a PDB file. Subsequently, the binding center was assigned to the frame corresponding to residue Cys145, as has been suggested previously [46] . Similarly, a comparative analysis was carried out between the docking using the scoring functions included in the Molegro Molecular package (MolDock Score, Rank Score and Plants Score). We sought to increase the stringency of the predictions made in this study by using more than one scoring function to determine the best coupled confirmation. In this sense, if most of the scoring functions predict similar trends for the same pose and/or ligand, the energetic mean of various scoring functions could help to optimize the prediction at the probabilistic level and classify the coupled complexes correctly from a set of potential ligand binding conformations and to classify the different conformers based on binding mode (depending on energy score values). In fact, it has been suggested to perform molecular couplings and to apply more than one scoring function like the ones proposed in this study in efforts that consider covalent and non-covalent systems in order to further improve docking accuracy [47, 52] . Overall, predicting differences in covalent and non-covalent binding free energy contributions for inhibitors could be a plausible explanation for their in vitro differences in antiviral activity as observed in in vitro assays. This indicates that covalent and non-covalent binding free energy contributions are important in the binding affinity of an inhibitor towards its target [35] . Therefore, other authors have applied more than one scoring method for the selection of the most promising candidates during covalent and non-covalent bonding molecular dockings [53] (see Supplementary material - Table 1S ).

J o u r n a l P r e -p r o o f It is important to note that even though the interactions presented in the 2D diagram appear weaker than the reported relative coupling energies, non-covalent free energy contributions are important in the binding affinity of an inhibitor towards its target as has been reported [35] In addition, the 2D diagrams were made from BIOVIA Simulations were performed for docking to sample the minimum energy conformations and validate if there was a structural and conformational alteration of the mM pro complexes, as well as the stability of the ligand binding (see Fig.6 ). For a protein-ligand complex, the MD system went through three phases: the relaxation phase, the equilibrium phase (composed of two equilibrium runs), and the production phase for sampling the trajectories of interest as suggested previously [46, 50, [54] [55] [56] [57] [58] .

MD simulation of crystal structures was carried out in an explicit water system.

Specifically, the solvation of the system was performed in an 8.0 Å solvation box. Our In protein systems, TIP3P water is suitable for the PME. The cutoff distance of the Van der Waals interaction was 14.0 Å. After applying the steepest descent algorithm to meet the energy minimization of all systems, the canonical NVT set was applied for 100 ps balanced of all systems and thermalized at 300 K. This was followed by a second run of 100 ps on the isothermal-isobaric NPT set to equilibrate the system at 1 atm and 300 K and progressively drive the equilibrium of each system. The SHAKE algorithm (used to satisfy link geometry constraints) was applied to the system and the time step was set to 2 fs. Finally, the production step was performed with the output of the NPT set, which was used as the initial configuration of an MD production series at a constant temperature of 300 K for a total simulation time of 100 ns.

Minimum energy structures in PDB format were obtained every 10 ns as target structures extracted from a 100 ns trajectory to be used in the following analyses. All MD simulations and additional adjustments were performed using cosgene/myPresto. In this sense, we worked with each of the avermectin homologues (B1a and B1b) and showed these molecules establish a thermodynamically favorable (ΔG ≈ -10 kcal/mol) and stable docking with mM pro with an RMSD ≈ ≤ 4 Å, noting that a minor fluctuation of the mM pro + B1b complex was predicted. The relative binding energy of the PF-07321332 control was ΔG ≈ -8 kcal/mol showing stable binding with mMpro as expected (RMSD ≈ ≤ 4 Å) (see Fig.6 ).

J o u r n a l P r e -p r o o f 17 

The partial molar volume of protein by definition [12, 13, 40, 43, 61] is considered to be constituted of two volumetric contributions, a molecular or geometric contribution (volume of van der Waals and volume of internal voids ) and a non-intrinsic contribution〈 〉 (volumetric contribution from repulsive and attractive interactions ):

is the van der Waals volumes of all the protein constitutive atoms of protein and is volume of cavities within the protein from imperfect atomic packing, dependent of temperature, proportional to molar mass , and equal to the geometric volume of protein impenetrable to surrounding solvent molecules [13] [14] [15] [16] [17] 19, 39, 45] . Where is the cavity volume that hosts the solute and by definition is [39, 54, 63, 64] . In the SPT model proposed by Lee-Graziano [40] [41] [42] [43] , the cavity volume of solute in a binary mixture of hard spheres can be estimated as 

Where is solvent accessible surface area of protein and as defined in the equation 7.

It is important to mention that this model although descriptive, has limitations due to the difficulty of exactly separating the border between the surface of the protein and the solvent (interfacial region). This model has been very useful to understand important processes such as hydration/dehydration of proteins involved in protein-ligand binding, protein-protein binding and chemical-induced unfolding [12, 14, 15, 39, 62, 63, 65] .

In computational biophysics studies, the molecular volume in a fluid is determined using the Voronoi-Delaunay model and this amount is 15% larger than the molecular volume previously defined in equations 1 and 2 [66] . This thermodynamic quantity is equivalent to the volume of a cavity dry or anhydrous (know as Voronoi volume) [62] in solution where the protein molecule has been placed [39, 62, 66] . This anhydrous volume will be defined herein after as of protein, which strictly and by definition, is equal to J o u r n a l P r e -p r o o f 20 (9) where is the internal void volume of protein and is the part of boundary empty space assigned to the solute [62] . If we assume that the protein-water interface division is the same in both models, then the three first components of the volumetric equation 3, and geometric equation 6 based in molecular dynamics can be empirically approximated as (10) And the partial molar volume of a protein can then be written as [see equation (3) and (9)]

The volume of interaction is the only term of 〈 〉 that contains information about the redistribution of water molecules between the hydration shell and bulk phase [15, 67] . On the other hand, in molecular dynamics [62, 68] the contribution is defined as the difference between the volume of the hydration shell around of solute and the occupied volume in the bulk by the same number of water molecules [62, 68, 69] .

Or alternatively (13) J o u r n a l P r e -p r o o f 21 If in the Voronoi-Delaunay approach, the border is defined as the Voronoi surface of the solute molecule, the equations 12 and 13 are equivalent. Also, the fact that in volumetric studies the term has an analog definition to equation 12 is very interesting.

Where the number of water molecules around of protein is , is the partial molar volumes of hydration and bulk water, respectively [15, 39] . Therefore, the empirical eq. 11 is a good approximation. In fact, recent studies have shown that 〈 〉 and are equal [62] .

In the high dilution regime the thermodynamic theory leads to the following expression for the partial molecular volume of proteins [70] ( )

Where is the dry volume of protein (here it is approximated to the Voronoi volume ( )), is the dry molar mass of protein, is the specific volume of bound water to the protein molecule and is the specific volume of nearby water. Chalikian and co-workers [17] have proposed that , while is the amount of grams of water associated with the protein per gram of protein, it is a measure of the degree of hydration of protein and it does not depend on temperature and concentration of solute [70] [71] [72] . In hydrodynamic studies [70] [71] [72] ] the volume of a particle in solution is called hydrodynamic volume and contains two components, one corresponding to the partial molar volume of the dry protein (also defined here as to maintain the same nomenclature) and the other due to volume of the hydration layer [70, 71] . 

Then, evaluating the change induced in each of the contributions to the components and 〈 〉 of apparent molar volume in the protein-ligand binding reactions provides important information on the structural and hydration changes that the protein undergoes in biological processes.

The proteins are dynamic molecules [73, 74] , and at thermal equilibrium, the intrinsic 

Where 〈 〉 is the mean square intrinsic volume fluctuaction, is the intrinsic volume of protein, is the Boltzmann contant, the absolute temperature, and the isothermal intrinsic compressibility coefficient. The mean square intrinsic volume fluctuaction can be statistically estimated from values and mean value of obtained from dynamic molecular simulation as:

Then 〈 〉 is a quantitative measure of protein dynamics. However, in many studies the amount used for comparative studies is: As the hydrodynamic volume , Simha factor and anhydrous volume (Voronoi volume) can be estimated using different methodologies such as HullRad [38] and 3Vee programs [36, 37] , all parameters described here in this section can be 

J o u r n a l P r e -p r o o f 27 If the change in the van der Waals volume is considered negligible [100] , then the difference in is due to binding-induced changes in the volume of the internal cavities of the mM pro and the contribution of volume that comes from the repulsive protein-water interactions [see equations (6 -8) and (25) With this in mind, we can suggest that , and making use of the values obtained for from SPT model (see Table 1 and eq.7 or 8) we find that for the difference between the mM pro -B1b complex and the mM pro , a value of for the difference between the mM pro -B1a It should be mentioned that homologues of ivermectin interact to the mM pro at different sites as shown in Fig.2, Fig.3, Fig.4 and Fig.5 . While the ligand PF-07321332 bind in the same site that the homologue B1a in the M pro monomer (see Fig.5 ).

This previous analysis was possible due to the knowledge of the Voronoi volume and the thermal volume for each protein involved in this analysis. While in the case of the second term of equation (24), we must first determine the hydration volume and from the hydrodynamic volume and Voronoi volume of each species. Once is known, it is possible to determine and then estimate both and in each case.

The values obtained for the hydrodynamic volume are shown in Table 1 Additional information about the changes induced in the hydration of the proteins was obtained from the degree of hydration measured as grams of water per gram of monomer , which was calculated with equation 19 and the values of the hydrodynamic volume and Voronoi volume (see Table 1 ). The results obtained are shown in Table 1 Another relevant aspect to highlight from these results is that the hydration obtained for the mM pro ( 0.57) is equal to the value reported for lysozyme ( 0.57) and close to porcine elastase ( 0.53) using other proteases models [95] . In fact, this degree of hydration is high and any change in this can significantly affect the molecular recognition events for this protease, since it is known that for this type of biochemical reaction in proteins, hydration is an extremely important driving force [11] .

It is very important to note there are no reports of this type of data for proteases associated with SARS-CoV-2 and in particular for M pro .

J o u r n a l P r e -p r o o f 32

The interaction from the specific and non-specific van der Waals forces of water as a solvent with the polar, non-polar and charged surface residues of the protein is quantified in the attractive volumetric contribution (see eq. 14 and 16) to the partial or apparent molar volume . This volumetric quantity therefore is a measure of the contraction of the volume of water as a solvent in the proximity of the charged, and polar groups of the protein [15, 39, 61, 65] . Hence its importance to evaluate the induced changes in protein hydration due to ligand binding. Table 1 shows the values obtained for this volumetric property, demonstrating the following trend (mM pro -B1b) (mM pro ) (mM pro -PF) (mM pro -B1a).

The difference between mM pro -B1b and mM pro is 247.26 , + 52.66

between mM pro -B1a and mM pro and 15.78 between mM pro -PF and mM pro These results suggest that the binding of ivermectin B1b to mM pro induces a redistribution of the charged and polar surface groups, favoring the exposure and the population of these to the solvent. While in the case of the complex mM pro -B1a and PF-07321332, the exposure of these surface amino acid residues is disfavored and the exposure of the nonpolar ones could be favored, which possibly affects the solubility of the protein as occurs with hydrophobic solutes that tend to aggregate in water [106, 107] .

It is important to see the consistency of these results by remembering that B1a and PF-07321332 bind at the same site within the protein mM pro . In any case, the little water that surrounds the protein or the site where B1a or PF-07321332 are bound should be dispersed or poorly structured. Thus affecting the degree of local hydration could be important. It has been proposed that water plays an important role in the catalytic site in M pro [20] . These results are also important since depending on which region this J o u r n a l P r e -p r o o f 33 hydration or dehydration occurs in the protein, the flexibility, stability and mechanisms of recognition of the protein can be significantly affected [68, 74, 79, 82, 83] . The homodimerization of this monomer to form the biologically active species could also be affected [29] . In this direction, the results obtained are interesting as the negative value of observed for the formation of the mM pro -B1b complex is associated with a high hydration which suggests lower aggregation tendency, diminution in hydrophobicity and higher charge surface in proteins [65] . In contrast, in the case of homologue B1a

and PF-07321332, the formation of the mM pro -B1a complex leads to a positive value of which is associated with low hydration and as consequence a higher aggregation tendency or an increased tendency to form aggregates or disordered clusters (an increase in hydrophobicity and lower charge in surface of proteins) [65] . In any case, the results suggest the negative effect on the homodimerization reaction. It is very important to mention that these results support the reports of the inhibition of the dimerization reaction of mM pro by binding of the PF-07321332 drug [20] [21] [22] [23] [24] . Therefore, as PF-07321332 and B1a bind at the same site, it is possible to propose that this ivermectin homologue has a similar effect.

It is very interesting that the trend observed for was similar to obtained for the electric dipole moment fluctuations on the surface (see eq.23) of the protein mM pro -B1b the protein-water interface [108] . It should also be noted that the values obtained for are in the expected order for a protein with a mass of 33.4 KDa [95] .

The value of the mM pro -B1a complex is 0.18% lower than that obtained for mM pro and mM pro -PF complex while it is 2.87% lower when compared to the mM pro -B1b complex. Although these percentages appear to be very small, there are no other reports of comparative studies between protein-ligand complexes and free protein that allow an analysis of this situation. It is clear that more studies on this topic are necessary.

However, these small differences between makes sense if it is thought to be related to changes in the populations of exposed polar amino acid residues on the surface of the protein by binding to ivermectin homologues and PF-07321332, and it is expected that these polar residue populations change by only a small fraction during global conformational changes.

These results, together with those obtained for and , suggest that each protein- 

Graziano reported that the non-intrinsic volumetric contribution 〈 〉 at infinite dilution in water for mixture of hard spheres is always a small and positive quantity [40] . We have recently confirmed this for model globular proteins [12, 61] . This is because this contribution to is the result of a balance between the repulsive The difference in the sign of this parameter is due to the fact that in the two firsts cases where 〈 〉 the repulsive volumetric component water-protein dominates in the unbound monomer, while in the last case this repulsive component dominates in the complex, which reflects that the mechanism involved in the volume change for the formation of the complexes is different in each case. It is interesting to note that although data suggest the B1a homologue and PF-07321332 bind in the same site in mM pro (see Figure 5 ), the property 〈 〉 has a similar magnitude but a different sign.

Based on equation 11 and using the data for Voronoi volume and interaction volume of Table 1 , we have determined the partial molar volume in physiological medium of mM pro , the mM pro -B1a, mM pro -B1b and mM pro -PF complexes. In the same Whereas, for B1a and PF-07321332 binding the volume of cavities within the protein may have the greatest effect on the change in volume.

As previously discussed, the volumetric component contributing the most to in B1a/mM pro and PF/mM pro is the cavity volume or internal voids , whereas the interaction volume at 100 ns is the major contribution to in the case of B1b/mM pro . Moreover, the geometric component is within of intrinsic volume , while contributes in balancing .

By analyzing the data presented in Table 2 previously reported [79] . The mM pro -B1a complex and the mM pro -PF complex are 0.6 % and 7% less dynamic than the native monomer, respectively, while the mM pro -B1b complex is 4.5% more dynamic than the native monomer. Thus, the fraction of volume fluctuaction 〈 〉 in all cases was around 0.6%, which is close to that reported for J o u r n a l P r e -p r o o f 40 other smaller proteins than mM pro (34.5 kDa) [78, 79] . Nonetheless, it is twice the mean value reported for a larger set of proteins that includes proteins of greater mass than mM pro [89] . It is very interesting to see that the complexes with the smallest volume fluctuation (〈 〉 , e.g., the mMpro-PF and mMpro-B1a complexes, also have the lowest values.

Since is a measure of the shrinkage of the hydration layer due to the attractive interactions of the protein surface groups with water molecules, and since the major contribution comes from the presence of charged groups exposed on the protein surface to the solvent [14, 15] , then the fraction of charged groups in these complexes should be low compared to the monomer. Therefore, the interfacial tension between the protein and water, must be high relative to the native monomer. The opposite is true for the mM pro -B1b complex, which exhibits the highest structural dynamics and the highest value of all species, suggesting a high fraction of charged surface groups exposed to water. Accordingly, we can expect a lower protein-water interfacial tension of all values. It should be noted that a low value of is associated with a high degree of protein hydration and vice versa, being in agreement with the adhesive-cohesive model for protein compressibility, in which the attractive forces of water compete with the intraprotein interactions favoring folding proposed by Dadarlat and Post [75] .

To support this idea, we estimate the interfacial tension using the value of 〈 〉 and the average radius of the protein 〈 〉 (see Table 1 ) and the model proposed by Lee [110] .

Here,

and 〈 〉 √ , the results obtained showed the following trend for , mM pro -PF( ) mM pro -B1a( ) mM pro ( ) mM pro -B1b( ) , which supports what was previously discussed.

Therefore, the results obtained in this work showed that 〈 〉 , , and are correlated and agree with the Dadarlat-Post model [75] . This model suggests also that high compressibility, resulting from large fluctuations in the molecular volume of the protein, corresponds to an enthalpically less stable protein [75] . Another possible connection may come from the inverse relationship between protein flexibility and stability [82, 111] . Moreover, it has been reported that the hydration of a protein is favored as its rigidity increases [112] . Based on these arguments, we propose that the mM pro -B1a and mM pro -PF complexes are more enthalpically stable than the unbound monomer, whereas the mM pro -B1b complex, is enthalpically unstable with respect to the native monomer. We propose that the increase or decrease in monomer structural dynamics induced by PF-07321332 or ivermectin binding has some inhibitory effect on homodimer formation [20] [21] [22] [23] [24] . As demonstrated, complex formation involves a change in volume fluctuation. It is also observed that this change correlates with changes in hydration, internal cavity distribution, interfacial tension and enthalpy stability, since some of these factors play an important role in relevant biological processes such as biological recognition. It is interesting to consider that these volumetric and thermodynamic properties depend on the binding site in the monomer. For example, the similarity in mechanism between PF-07321332 and the B1a homologue that bind at the same site, PF-07321332 ligand is known to bind covalently at the catalytic site of M pro [24] , theoretical studies using directed docking suggest that this occurs by a two-stage process. The first involves the formation of a non-covalent complex at the same site and followed by an electrophilic coupling between the thiol group of a cysteine at the catalytic site and the nitrile group from the PF-07321332 a forming a thio-imidazolium bond giving stability to the covalent complex [20] . Although the covalently linked species is thermodynamically less stable than the non-covalent complex, this covalent complex is the only species detected so far [24] . This covalent complex indicates that the blocking of the catalytic site is its main disturbing action on the M pro protease [20] [21] [22] [23] [24] . However, in this study we have found through the use of blind docking that PF-07321332 can bind more strongly through non-covalent interactions in another distant region in the monomer mMpro, which in turn is the same site where the ivermectin B1a homologue also binds. An interesting result is that all the theoretical methods used suggest that this bonding is thermodynamically more feasible at that site in a noncovalent way than the covalent complex in the catalytic region (see supplementary material - Table 1S ). There is no report that has considered the effect of this on The PF-07321332 binding region different from those evaluated here is described (see Supplementary material - Table 1S ). It is important to highlight that the interactions of the PF-07321332 in the site covalently linked with the aminoacid residues coincide very well with that previously reported [20] [21] [22] [23] [24] .

Despite the predictions made in this study about a potential theoretical effect at the volumetric level of ivermectin similar to that of PF-07321332 on M pro , it is important to point out that the usefulness of ivermectin in the treatment of COVID-19 is a subject of controversy [113] . Additionally, data from the I-Tech randomized clinical trial do not support the use of ivermectin for COVID-19 [114] . While there are studies showing binding of ivermectin to M pro , most of them are in silico coupling-based or in vitro fluorescence-based assays that are prone to false-positive results (especially given the size of this drug molecule). A recent study also suggested that this drug does not show any significant activity in human airway-derived cell models [115] .

J o u r n a l P r e -p r o o f 46 On the other hand, the controversy increases because ivermectin has been included in clinical trials that have shown promising results [116] (https://clinicaltrials.gov/ct2/show/NCT04425863 -ClinicalTrials.gov Identifier:

NCT04425863) and in collective reviews of multiple efforts that have suggested that ivermectin may have a prophylactic effect and would be a strong candidate for clinical trials to treat SARS-CoV-2 [117] , as also observed in systematic reviews and metaanalyses of randomized clinical trial studies [118] . especially when seeking to evaluate the antiviral activity of drugs [123] . Our team is currently working in this direction.

In this sense, multiple authors have suggested continuing studies of ivermectin and its interactions with the various targets of SARS-CoV-2 in order to use it as a model to guide efforts towards the development of new compounds and treatment strategies [119, 120, 124] as has already been done [125] . It is important to point out that our study did not consider the competition between ligands, so it is recommended to carry out a competitive binding study taking into account the inhibition kinetics of each compound and evaluating the effect of the association and dissociation mechanisms of each ligand on the sites of interest.

However, this study focused on examining the complexes formed by the drugs once they were bound independently against the same target protein, with the aim of once the complexes were formed, to be able to evaluate by molecular dynamics the strength of the union predicted by the dockings and also take minimum energy structures every 10 ns and at 100 ns to perform an analysis of the changes induced in the volumetric and hydrodynamic properties of the M pro using different validated theoretical models such as HullRad [38] and 3vee.mol [37] .

Specifically, in the case of PF-07321332, the complex was reproduced as reported in the literature in which it is reported that its complex with M pro is mediated by a covalent bond and therefore its binding is covalent at the site [24] , at the same time non-covalent bonds were also predicted. While in the case of ivermectin, which is a mixture of two homologues (80:20 mixtures avermectin B1a and avermectin B1b, which differ in the presence of a secbutyl and an isopropyl group, at the C25 position, respectively), what is known in the literature, including what was reported by our work, is that its binding establishes only non-covalent bonds in the M pro [7, 38, 126] .

In this sense, we focus on further studying whether PF-07321332 and ivermectin homologues, once bound to the target protein using the same conditions and degrees of In terms of docking, the docking reported for PF-07321332 could be reproduced, and a differential docking for ivermectin homologs on M pro was predicted. It was found that compound B1a and PF-07321332, although they bind differently at the same site, both induce changes in volumetric properties in a similar but not identical way. While the complex formed by the monomer and compound B1b, bound to a different site, led to volumetric and hydrodynamic changes different from those of B1a and PF-07321332.

These observations are important because in the literature it is known that the volumetric and hydrodynamic study can provide information on possible mechanisms of action of various compounds [132] .

In this way, driven by the controversy that exists around ivermectin, it was sought to know if when comparing the volumetric and hydrodynamic changes induced by its counterparts, the resulting predictions had any difference or similarity with respect to the changes predicted for the complex formed by PF-07321332, which is an authorized and validated drug for the treatment of COVID-19 and is known to bind covalently and affect homodimerization [24, 133] . Additionally, this work provides the first report of the effect of PF-07321332 on M pro at a volumetric and hydrodynamic level.

The importance of our predictions is based on the fact that when an accurate thermodynamic analysis of the volumetric fluctuation is attempted to find its Finally, as we focus on analyzing the monomer as has been done in other theoretical studies [33] , in order to describe the perturbation that each of these ligands could induce on the main block of the main protease (the monomer) of SARS-CoV-2, and how these interactions can cause volumetric and hydrodynamic perturbation of homodimerization, especially, since M pro is known to depend on homodimerization for its biological activity [29, 30] . However, as the interest of various efforts is the inhibition of M pro by blocking or perturbing the active site [33] , we recommend, to perform the analyzes proposed here on the dimeric form of M pro , since the catalytic pocket of a monomer is capped at the N-terminus of the adjacent unit, as has also been done [136] .

The results obtained in this work show that the binding of Paxlovid (PF, PF- The volumetric results support our previous observations from a thermodynamic and structural perturbation point of view using various computational tools based on molecular dynamics that each homologue binds at different sites and disrupts the global conformation differently. Also, the present study supports the experimental and theoretical reports previously suggesting that PF and B1a perturb the structure of mM pro in a similar mechanism to each other but by a different mechanism to B1b. Under this same study strategy, it was found that PF-07321332 can also covalently bind to another region of the monomer (a drug multi-site for M pro ), but both complexes have similar structural flexibility and volumetric properties. Finally, we consider that this type of study can help to understand the mechanism by which a ligand can block the homodimerization of this important monomer to form the dimeric protease dM pro and in turn help in studies of activity-structure relationships for the design of new drugs.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sector. Writing-Reviewing and Editing, Conceptualization, Methodology, Investigation.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Using informative features in machine learning based method for COVID-19 drug repurposing

Comparative study of the interaction of ivermectin with proteins of interest associates with SARS-CoV-2: A computational and biophysical approach

Drug repurposing studies targeting SARS-CoV-2: an ensemble docking approach on drug target 3C-like protease (3CL pro )

An in-silico analysis of ivermectin interaction with potential SARS-CoV-2 targets and host nuclear importin α

Repurposing of FDA-approved drugs against active site and potential allosteric drug-binding sites of COVID-19 main protease

Exploring the binding efficacy of ivermectin against the key proteins of SARS-CoV-2 pathogenesis: and in silico approach

Structural deformability induced in proteins of potential interest associated with COVID-19 by binding of in ivermectin: Comparative study based in elastic networks models

Local conformational fluctuations can modulate the coupling between proton binding and global structural transitions in proteins

Excluded volume contribution to cosolvent-mediated modulation of macromolecular folding and binding reactions

Protein-solvent preferential interactions, protein hydration, and the modulation of biochemical reactions by solvent components

Does the release of hydration water come with a Gibbs energy contribution?

Surface behavior of BSA/water/carbohydrate systems from molecular polarizability measurements

How large are the volume changes accompanying protein transitions and binding?

Interactions of glycine betaine with proteins: Insights from volume and compressibility measurements

Volumetric characterization of tri-Nacetylglucosamine binding to lysozyme

Volume of Hsp90 protein-ligand binding determined by fluorescent pressure shift assay, densitometry, and NMR

Probing the binding ability of vitamin B 1 with bovine serum albumin: Calorimetric, light scattering, spectroscopic and volumetric studies

Hydration of proteins: Excess partial volumes of water and proteins

Conformational change of ovalbumin induced by Surface cavity binding of N-phthaloyl gamma-aminobutyric acid derivative: A study theoretical and experimental

Characterization of the non-covalent interaction between the PF-07321332 inhibitor and the SARS-CoV-2 main protease

Exploring the binding mechanism of PF-07321332 SARS-CoV-2 protease inhibitor through molecular dynamics and binding free energy simulations

Supervised Molecular Dynamics (SuMD) Insights into the mechanism of action of SARS-CoV-2 main protease inhibitor PF-07321332

Considerations for the discovery and development of 3-chymotrypsin-like cysteine protease inhibitors targeting SARS-CoV-2 infection

Crystal structure of SARS-CoV-2 main protease in complex with protease inhibitor PF-07321332

Covalent narlaprevir-and boceprevir-derived hybrid inhibitors of SARS-CoV-2 main protease: room-temperature X-ray and neutron crystallography, binding thermodynamics, and antiviral activity

Benchmarking of different molecular docking methods for protein-peptide docking

Potential therapeutic target identification in the novel 2019 coronavirus: insight from homology modeling and blind docking study

In silico evidence of beauvericin antiviral activity against SARS-CoV-2. Computers in biology and medicine

Targeting the dimerization of the main protease of coronaviruses: A potential broad-spectrum therapeutic strategy

Impact of dimerization and N3 binding on molecular dynamics of SARS-CoV and SARS-CoV-2 main proteases

Drug design and repurposing with DockThor-VS web server focusing on SARS-CoV-2 therapeutic targets and their non-synonym variants

Improving blind docking in DOCK6 through an automated preliminary fragment probing strategy

Site mapping and small molecule blind docking reveal a possible target site on the SARS-CoV-2 main protease dimer interface

Ligand and structure-based virtual screening applied to the SARS-CoV-2 main protease: an in silico repurposing study

Covalent and non-covalent binding free energy calculations for peptidomimetic inhibitors of SARS-CoV-2 main protease

MOLE: A Voronoi diagram-based explorer of molecular channels, pores, and tunnels

3V: Cavity, channel and cleft volume calculator and extractor

HullRad: Fast calculations of folded and disordered protein and nucleic acid hydrodynamic properties

On empirical decomposition of volumetric data

Non-intrinsic contribution to the partial molar volume of cavities in water

Partial molar volumen from the hard-sphere mixture model

Partial molar volumen of n-alcohols at infinite dilution in water calculated by means of scaled particle theory

Volume-related properties of thiophene and furan-2-carboxaldehyde phenylhydrazone derivatives in DMSO: A discussion about non-intrinsic contribution

Essential oils as an effective alternative for the treatment of COVID-19: Molecular interaction analysis of protease (M pro ) with pharmacokinetics and toxicological properties

Exploration of inhibitory action of Azo imidazole derivatives against COVID-19 main protease (M pro ): A computational study

New tyrosinases with putative action against contaminants of emerging concern

Identification of irreversible protein splicing inhibitors as potential anti-TB drugs: insight from hybrid non-covalent/covalent docking virtual screening and molecular dynamics simulations

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

Exploration of new druglike inhibitors for serine/threonine protein phosphatase 5 of Plasmodium falciparum: a docking and simulation study

Repurposing FDA-approved drugs to fight COVID-19 using in silico methods: targeting SARS-CoV-2 RdRp enzyme and host cell receptors (ACE2, CD147) through virtual screening and molecular dynamic simulations. Informatics in medicine unlocked

Open Babel: An open chemical toolbox

Molecular dynamics study on graphene-mediated pyrazinamide drug delivery onto the pncA protein

Repurposing known drugs as covalent and non-covalent inhibitors of the SARS-CoV-2 papain-like protease

Size dependence of cavity volume: A molecular dynamics study

Evaluation of drug repositioning by molecular docking of pharmaceutical resources available in the Brazilian healthcare system against SARS-CoV-2. Informatics in medicine unlocked

Loop dynamics behind the affinity of DARPins towards ERK2: Molecular dynamics simulations (MDs) and elastic network model (ENM)

Reorientational dynamics of molecules in liquid methane: A molecular dynamics simulation study

Electrostatic properties of water models evaluated by a long-range potential based solely on the Wolf charge-neutral condition

myPresto/omegagene 2020: a molecular dynamics simulation engine for virtual-system coupled sampling

Can non-steroidal anti-inflammatory drugs affect the interaction between receptor binding domain of SARS-COV-2 spike and the human ACE2 receptor? A computational biophysical study

Non-intrinsic contribution to the limiting partial molar volume of globular proteins in water: a study comparative between a new refractometric strategy and densitometric classical approach

Exploring volume, compressibility and hydration changes of folded proteins upon compression

Flexibility and hydration of amphiphilic hyperbranched arabinogalactan-protein from plant exudate: A volumetric perspective

Partial molar volumes of molecules of arbitrary shape and the effect of hydrogen bonding with water

Interaction volume is a measure of the aggregation propensity of amyloid-β

Binding of bovine pancreatic trypsin inhibitor to trypsinogen: Spectroscopic and volumetric studies

The hydration of globular proteins as derived from volume and compressibility measurements, cross correlating thermodynamic and structural data

Disentangling volumetric and hydrational properties of proteins

Volumetric properties of human islet amyloid polypeptide in liquid water

Viscosity and the shapes of macromolecules: A physical chemistry experiment using molecular-level models in the interpretation of macroscopic data obtained from simple measurements

A viscometric approach of pH effect on hydrodynamic properties of human serum albumin in the normal form

On the hydrodynamics and temperature dependence of the solution conformation of human serum albumin from viscometry approach

Exploring cavity dynamics in biomolecular systems

Compressibility of cavities and biological water from Voronoi volumes in hydrated proteins

Insights into protein compressibility from molecular dynamics simulations

Hydrational and intrinsic compressibilities of globular proteins

The influence of correlated protein-water volume fluctuations on the apparent compressibility of proteins determined by ultrasonic velocimetry

Evaluation of intrinsic compressibility of proteins by molecular dynamics simulation

Compressibility of the protein-water interface

Native state volume fluctuations in proteins as a mechanism for dynamic allostery

Change in vibrational entropy with change in protein volume estimated with mode Grüneisen parameters

Native protein fluctuations: The conformational-motion temperature and the inverse correlation of protein flexibility with protein stability

Long-range correlation in protein dynamics: Confirmation by structural data and normal mode analysis

Pressure perturbation: A prime tool to study conformational substates and volume fluctuations of biomolecular assemblies

Molecular science of fluctuations toward biological functions

Inhomogeneous molecular density: Reference packing densities and distribution of cavities within proteins

Dynamic void distribution in myoglobin and five mutants

Protein fluctuations and the thermodynamic uncertainty principle

Pressure-A gateway to fundamental insights into protein solvation, dynamics, and function

Compressibility-structure relationship of globular proteins

Ergometric studies of proteins: New insights into protein functionality in food systems

Fatty acid and retinolbinding protein: Unusual protein conformational and cavity changes dictated by ligand fluctuations

Protein fluctuations and cavity changes relationship

Quantifying the protein core flexibility through analysis of cavity formation

Protein binding pocket dynamics

Electrical fluctuations on the surfaces of proteins from hydrodynamic data

WAXSiS: a web server for the calculation of SAXS/WAXS curves based on explicit-solvent molecular dynamics

Validating solution ensembles from molecular dynamics simulation by wide-angle X-ray scattering data

Hinge Prot: automated prediction of hinges in protein structures

Volumetrically derived thermodynamic profile of interactions of urea with a native protein

Interactions of urea with native and unfolded proteins: A volumetric study

On volume changes accompanying conformational transitions of biopolymers

Investigations on the pH-dependent binding of sodium valproate with bovine serum albumin: A calorimetric, spectroscopic and volumetric approach

Water is an active matrix of life for cell and molecular biology

Water dynamics in the hydration shells of biomolecules

Hydrogen-bond dynamics and energetics of biological water

Temperature dependence of the pairwise association of hard spheres in water

Energetics of the contact minimum configuration of two hard spheres in water

Dipolar susceptibility of protein hydration shells

On the magnitude of border thickness in the partial molar volume of cavities in water

Calculations of fluctuations for globular models

The complex inter-relationships between protein flexibility and stability

Protein hydration thermodynamics: The influence of flexibility and salt on hydrophobin II hydration

Flawed ivermectin preprint highlights challenges of COVID drug studies

Efficacy of ivermectin treatment on disease progression among adults with mild to moderate COVID-19 and comorbidities: The I-TECH randomized clinical trial

Moxidectin and ivermectin inhibit SARS-CoV-2 replication in Vero E6 cells but not in human primary airway epithelium cells

Ivermectin, aspirin, dexametasone and enoxaparin as treatment for COVID 19

Repositioning Ivermectin for Covid-19 treatment: Molecular mechanisms of action against SARS-CoV-2 replication

Ivermectin and outcomes from Covid-19 pneumonia: a systematic review and meta-analysis of randomized clinical trial studies

Prospective mode of action of Ivermectin: SARS-CoV-2

Remdesivir-ivermectin combination displays synergistic interaction with improved in vitro activity against SARS-CoV-2

Evaluation of SARS-CoV-2 entry, inflammation and new therapeutics in human lung tissue cells

Entrectinib-A SARS-CoV-2 Inhibitor in Human Lung Tissue (HLT) Cells

Comparative study of SARS-CoV-2 infection in different cell types: Biophysical-computational approach to the role of potential receptors

Antiviral Activity of Repurposing Ivermectin against a Panel of 30 Clinical SARS-CoV-2 Strains Belonging to 14 Variants

Two antioxidant 2, 5-disubstituted-1, 3, 4-oxadiazoles (CoViTris2020 and ChloViD2020): successful repurposing against COVID-19 as the first potent multitarget anti-SARS-CoV-2 drugs

A bioinformatics study of structural perturbation of 3CL-protease and the J o u r n a l P r e -p r o o f HR2-domain of SARS-CoV-2 induced by synergistic interaction with ivermectins

Computationally driven discovery of SARS-CoV-2 M pro inhibitors: from design to experimental validation

Discovery of S-217622, a Noncovalent Oral SARS-CoV-2 3CL Protease Inhibitor Clinical Candidate for Treating COVID-19

In search of non-covalent inhibitors of SARS-CoV-2 main protease: Computer aided drug design using docking and quantum chemistry

Structure of papain-like protease from SARS-CoV-2 and its complexes with non-covalent inhibitors

Crystal structure of SARS-CoV-2 main protease in complex with the non-covalent inhibitor ML188

Intrinsically Disordered Proteins: Targets for the Future

Inhibition of the main protease of SARS-CoV-2 (Mpro) by repurposing/designing drug-like substances and bioactive compounds

Direct evaluation of polypeptide partial molar volumes in water using molecular dynamics simulations

Molecular dynamics studies of the nucleoprotein of influenza A virus: role of the protein flexibility in RNA binding

Targeting SARS-CoV-2 main protease by teicoplanin: A mechanistic insight by docking, MM/GBSA and molecular dynamics simulation