pubs.acs.org/IC Article Getting Deeper into the Molecular Events of Heme Binding Mechanisms: A Comparative Multi-level Computational Study of HasAsm and HasAyp Hemophores Laura Tiessler-Sala, Giuseppe Sciortino, Lur Alonso-Cotchico, Laura Masgrau, Agustí Lledós, and Jean-Didier Maréchal* Cite This: Inorg. Chem. 2022, 61, 17068−17079 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: Many biological systems obtain their activity by the inclusion of metalloporphyrins into one or several binding pockets. However, decoding the molecular mechanism under which these compounds bind to their receptors is something that has not been widely explored and is a field with open questions. In the present work, we apply computational techniques to unravel and compare the mechanisms of two heme-binding systems, concretely the HasA hemophores from Gram negative bacteria Serratiamarcescens (HasAsm) and Yersinia pestis (HasAyp). Despite the high sequence identity between both systems, the comparison between the X-ray structures of their apo and holo forms suggests different heme-binding mechanisms. HasAyp has extremely similar structures for heme-free and heme-bound forms, while HasAsm presents a very large displacement of a loop that ultimately leads to an additional coordination to the metal with respect to HasAyp. We combined Gaussian accelerated molecular dynamics simulations (GaMDs) in explicit solvent and protein−ligand docking optimized for metalloligands. GaMDs were first carried out on heme-free forms of both hemophores. Then, protein−ligand dockings of the heme were performed on cluster representatives of these simulations and the best poses were then subjected to a new series of GaMDs. A series of analyses reveal the following: (1) HasAyp has a conformational landscape extremely similar between heme-bound and unbound states with no to limited impact on the binding of the cofactor, (2) HasAsm presents as a slightly broader conformational landscape in its apo state but can only visit conformations similar to the X-ray of the holo form when the heme has been bound. Such behavior results from a complex cascade of changes in interactions that spread from the heme-binding pocket to the flexible loop previously mentioned. This study sheds light on the diversity of molecular mechanisms of heme-binding and discusses the weight between the pre-organization of the receptor as well as the induced motions resulting in association. ■ INTRODUCTION The inclusion of inorganic moieties in proteins leads to unique to contribute to crucial biological functions, like transport/ storage of O2, electron transfer, or redox catalysis.6 More structural and catalytic features that overcome the limitations recently, heme has also been found to participate as a signaling of a purely organic living world. More than 40% of all proteins molecule in cellular processes like transcription regulation, in the cells exploit one or several metals to perform their protein complex assembly, microRNA processing, or cell functions.1 Decoding the interactions between inorganic growth and differentiation.7 Therefore, heme is a prototypical moieties and biological partners is therefore a fundamental system to study the binding of metallic cofactors to proteins, question in understanding the origin of life as well as in although there is no clear molecular description of heme opening new biotechnological routes, like in the design of uptake. artificial metalloenzymes obtained from the insertion of One could expect several scenarios in heme-binding inorganic homogeneous catalysts into proteic hosts.2−4 mechanisms that include conformational selection and induced However, metal-mediated binding processes are one of the most complex questions molecular biology could address. The biological−inorganic interplay overcomes standard knowledge Received: June 23, 2022 of both chemical and biological sciences, and assessing the Published: October 17, 2022 relative contributions of both partners in the recognition process is extremely challenging.5 Iron protoporphyrin IX or heme b is one of the most ubiquitous metal-containing ligands in nature. It is well known © 2022 The Authors. Published by American Chemical Society https://doi.org/10.1021/acs.inorgchem.2c02193 17068 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article Figure 1. Overall structure of hemophore HasA with names of secondary structure regions and heme coordinators indicated. Sequence alignment between HasAsm (up) and HasAyp (down) in ClustalX colors. Figure 2. Structural overlap between apo (1YBJ) and holo (1DKH) forms of hemophore HasA from S. marcescens (left) and apo (4JES) and holo (4JET) forms from Y. pestis (right). Shifts between apo and holo forms in loop L1 are indicated. fit. Conformational selection sustains that several unbound preferably to one or several well pre-organized ones.8 Induced states of the protein exist in equilibrium and the ligand binds fit implies that the conformational change of the receptor is a 17069 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article product of the entrance of the ligand.9 Both phenomena are HasAsm.24 To date, the unique computational study on likely to participate in ligand binding with different weights HasAsm also involves TMD simulations for the heme transfer depending on the system, and in some cases, one may prevail to HasR, and no simulations have been reported on HasAyp.25 over the other.10,11 In spite of extensive research of substrate It is important to notice that TMD,26 highly informative for binding to hemoproteins as P450 or cytochromes,12 the molecular mechanisms, implies that simulations are forced molecular mechanism under which the heme binds to its toward a defined final structure by using steering forces. Such receptor has as yet been rarely explored.13 restraints do not allow us to assess the conformational space Spectroscopic and crystallographic data on heme proteins that the apo-protein explores in standard conditions and how have shown that the tertiary structure of the apo and holo this is related to the heme-binding mechanism. forms are in general similar.13 Two main mechanisms seem to This study herein pretends to give further insights into the irrupt, though. In some systems like myoglobine,14 cytochrome heme-binding process in hemophores considering that it b5,15,16 or cytochrome b562,17,18 there are some secondary represents a ground material for similar systems. The structure rearrangements that occur because of the binding and underlying methodology consists in providing a wide generally in the proximal side of the heme. In others, only conformational exploration of the apo-hemophore structures, subtle rearrangements of the secondary structure are intermediate heme-bound complexes, and the final complex to observed.19 The latter systems are generally referred to as evaluate the pre-organization degree and identifies possible transient heme-binding proteins because they allow fast induced effects. To do so, we bridge together all-atom MD association−dissociation mechanisms. Among those systems with protein−ligand docking. In particular, we apply Gaussian are heme-chaperones, like HemS systems,19 or hemophores accelerated MD (GaMD), which allows simulations of long- from the HasA family. The crystal structures of HasA members range motions and explores the conformational space of the present, however, some variability that could give further system without imposing specific geometric restraints. More- insights into more complex heme-binding mechanisms. over, in both approaches, we use updated parameters for HasAs are extracellular heme-binding proteins that Gram- dealing with the metallic moieties, overcoming the challenges negative bacteria use for their heme uptake. They are able to still present for metal ion description in force field-based acquire free or hemoprotein-bound heme and to deliver it to approaches. The results give interesting insights into the heme- a specific receptor at the cell surface (HasR), whereby the binding mechanisms in heme binding proteins that could heme is internalized and used as an iron source.20 Interestingly, better envision protein engineering processes for heme- two hemophores from this family seem to have very distinctive containing enzymes and for artificial enzymes. heme-binding mechanisms: HasA from Yersenia pestis (HasAyp) and HasA from Serratia marcescens (HasAsm). HasA structure contains a α + β fold structure in which the ■ METHODOLOGY Overall Pipeline. The computational framework for this heme is found between loops L1 (27−43) and L2 (74−84) at work combines molecular dockings and GaMD simulations. It the interface of the α and β domains.20,21 Both have very starts with a crystallographic heme-free form available at the similar structures and a sequence identity of 31% (Figure 1), Protein Data Bank (PDB) to run GaMD simulations. The with the less-conserved part being located at loop L1.22 analysis of the trajectories, including cluster analysis, is Interestingly, experimental structures of the apo and holo followed by dockings of heme. From the docking results, forms of HasAyp and HasAsm show striking differences GaMD simulations on heme-bound intermediates are (Figure 2): (1) L1 position is almost invariant in HasAyp in performed to study the effect of heme binding to the protein apo and holo forms with a closed conformation upon the conformation (Scheme 1). heme-binding pocket, while a very large conformational change is observed in HasAsm, with an opened form in the apo Scheme 1. Multi-level Computational Protocol Followed structure and a closed form in the holo one with a motion of about 15 Å, and (2) this is also related to a difference in axial coordination of the iron in the two holo forms; the iron is bound to Tyr75 for the L2 loop in HasAyp, whereas it binds to both Tyr75 (L2) and His32 (L1) in HasAsm. This shows that very distinct binding mechanisms could occur although spectroscopic data and the comparison between the apo and holo structures are not able to provide a clear molecular understanding. In this regard, computational tools could be a very valuable asset. Most of the computational studies regarding hemophores have been performed on the hemophore HasA from Pseudomonas aeruginosa, a system that displays geometrical System Setup. The X-ray structures of the apo forms of features similar to those of HasAsm. In their first work, Rivera heme-containing proteins were obtained from the PDB.27 The et al. used targeted molecular dynamics simulations (TMD) PDB code 1YBJ was used for the apo form of HasAsm and and identified a series of interactions and motions in helices α2 4JES for the apo form of HasAyp. Calculations on the and α3 upon heme binding that could be involved in the experimental holo HasAsm were performed on the structure closing of loop L1.23 In a later study, the transition from the with the PDB code 1DKH. All structures were cleaned by apo to the holo conformation of HasAp was not observed in a removing crystallographic waters and small molecules present 100 ns classical MD. However, MD of the Arg33Ala mutant in the PDB file using UCSF Chimera.28 Hydrogen atoms were pointed that this residue could be important for controlling the added using Chimera, and webserver H++29 was also used to closing of loop L1. Strikingly, Arg33 is not present in double-check the protonation state of ionizable groups. 17070 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article Docking Calculations. Dockings were performed with performed followed by a production run of 800 ns. In all cases, GOLD5.2,30 using a simulation box of 10−15 Å and centered three replicas of each calculation were performed in order to at the binding site of each protein. Side-chain flexibility on the assure convergence and a full exploration of the conforma- binding site residues was considered when required using the tional space. default rotamer library. For higher accuracy, the minimum Analysis of GaMD Simulations. All GaMD trajectories number of operations was set to 100,000, and the number of were processed using cpptraj implemented in Ambertools18,34 the Genetic algorithm (GA) runs to 50. An optimized version and cluster analysis was performed with MDtraj.42 The of GoldScore as a scoring function capable of predicting trajectories were considered converged by assessing the metal−protein interactions was used.31 All the solutions were variation in root-mean-square deviation (rmsd) with respect analyzed using GaudiView.32 to the initial structure, all-to-all rmsd, RSMF, PCA, and a GaMD Molecular Dynamics Simulations. Uncon- clustering counting method.43 All calculations were performed strained enhanced sampling was performed with the Gaussian using carbon alpha. accelerated molecular dynamics (GaMD) method. This For interaction analysis, the getContacts.py44 script was used method consists in adding adaptively a harmonic boost and tuned in order to analyze hydrogen bonds, salt bridges, π− potential to smoothen the potential energy surface and to stacking, and hydrophobic interactions along the GaMD explore massive conformational landscape without predeter- simulations. All contact results from different replicas were mined reaction coordinates or any kind of restraints. combined to obtain the mean frequency of contacts through all Furthermore, GaMD simulations provide speedup simulations replicates. To study interaction differences between simulation and allow capturing events over a longer time. This method apo−holo pairs, the difference of interactions was calculated allows a larger exploration of the conformational space and and normalized. The results were simulated into UCSF allows us to identify different states of the biomolecules. A Chimera using pseudo-bond representation. To calculate free short conventional MD always proceeds with a GaMD energies, GaMD reweighting was performed with the calculation as a preparatory stage that collects potential PyReweighting toolkit and the Maclaurin method.45 statistics.33 This technique has been applied using AMBER18’s code using the coordinates extracted from a classical MD of 10−20 ns as a starting point.34 The parameter for the GaMD ■ RESULTS AND DISCUSSION In this section, first HasAyp and HasAsm hemophore are igamd was set to 3, which applies a force to dihedrals and to analyzed and discussed separately, and then, a comparison the total potential energy, and the threshold energy was set to between the two hemophores is presented. the lower bound, IE = 1. HasA from Y. pestis. As previously mentioned, the apo and GaMD simulations were prepared with the xleap35 using the the holo crystallographic forms of hemophore HasAyp have force field ff14SB for proteins, GAFF for non-standard similar structures, and loop L1 remains mostly static (Figure residues, ions94.lid for ions, and TIP3P for water. In the 2b). Still, we wanted to study better the conformational space case of metalloproteins, metal parameters were obtained using and assess if changes in L1 could exist with respect to the X-ray the MCPB.py approach.36 In MCPB.py, the charges were structures. For example, would it be possible that the loop calculated using RESP, and the force constants and equilibrium could eventually reach conformation like HasAsm in the parameters between the metal and the residues were calculated absence of heme? Three replicas of 800 ns GaMD simulations using the Seminario method.37 Optimization and frequency were performed. The combination of five analysis tools (PCA, calculations of heme at the DFT level were performed with clustering, rmsd, all-to-all rmsd, and RMSF) revealed that after Gaussian0938 in water solvent (SMD continuum model).39 800 ns, all three replica (2.4 μs) simulations had converged, The B3LYP hybrid functional including Grimme’s dispersion and the conformational space sampled minimal changes in the D340 was used with SDD + F (Fe) + 6-31G(d,p) as the main overall structure (Figure S1). group.41 Different Fe-oxidation states (Fe2+ and Fe3+) and The visual inspection of the trajectory and clusters of the multiplicity (low and high) were taken into account. For the GaMD simulations revealed how the tertiary structure of the sake of the presentation of the work, the results presented in protein, the β-sheet, and α-helices displays very little the main text correspond to calculations performed with the conformational changes (Figure 3). The loop and hairpin Fe-oxidation state, which was +3 high spin when pentacoordi- regions only show some degree of flexibility, as confirmed by nated and low spin when hexacoordinated; those are the RMSF analysis (Figure S2) in particular hairpin H3 apparently the most accepted oxidation and spin states from (average RMSF 2.86 Å) and loop L1 (average RMSF 1.87 Å). the literature. Still, Fe(II) simulations were also performed The latest oscillated during the simulation around its resting with the same spin state considerations (see ESI). position although with only small changes in its rearrangement All GaMD simulations were set up solvating the protein and could partially acquire the conformation of a small turn or using an explicit solvent approach, in which the protein was helix. embedded into a cubic box adding as counter ions Na+, Only in very brief sections of the trajectory, loop L1 reaches specifically 10 and 12 Na+ for HasAyp for apo and holo forms, either a more closed arrangement for which it could be respectively, and 13 and 15 Na+ for HasAsm for apo and holo foreseen that the heme could not be bound or a slightly forms, respectively. GaMD simulations were carried out under opened but never reaching the open conformation observed in periodic boundary conditions using AMBER18.34 For the the structure of HasA of S. marcescens. In this “closed” preparatory MD simulations, energy minimization was arrangement, hydrophobic interactions and hydrogen bonds performed to avoid steric clashes and relax the system, were observed during the whole simulation (>75%) in between followed by several equilibration steps in which the system was the α-helices and β-sheets. Interestingly, a network of several heated from 100 to 300 K. Finally, a production run of 100 ns π-contacts between aromatic residues of the pocket (Tyr55, was carried out. GaMD simulations started from coordinates of Phe43, Phe83, Tyr75, Phe50, and Phe77) maintains the heme- MD after 10−20 ns, and an equilibration of 50 ns was binding region in a closed conformation in the absence of 17071 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article the most flexible region in this hemophore due to the lack of strong intra-protein interactions. To gain further insights into the heme-binding process, protein−ligand docking of the heme cofactor was performed using the most representative protein conformations of the GaMD trajectory. The calculations were performed with the updated GoldScore scoring function for metalloligands.31 Taking into account that the GaMD of the apo-form shows substantial similarity with the holo structure available in the PDB, it was not surprising that all dockings showed excellent scoring values for heme binding (ca. 90 GoldScore units). The resulting complexes showed a structural arrangement similar to the experimental structure, including the presence of a coordination bound between the metal and the Tyr75 in all the cases. Finally, simulations with the heme bond were performed to ascertain how stable the heme−HasAyp complexes predicted in the dockings are. We envisioned that some of those Figure 3. Main clusters of GaMD simulation of the apo form of complexes may differ from the X-ray structure and expected to HasAyp. In red is represented the first cluster, and the remaining see how the structure of the protein is affected by the binding clusters are in gray. of heme. Three replicas of GaMD simulations starting from the heme-docking position were undertaken. The analysis shows that the theoretical holo forms of heme (Figure 4a). These interactions also may contribute to HasAyp do not present significant conformational changes. keeping the stability of the overall protein structure with The system tends to reach convergence only after 100 ns as almost no structural changes. When looking more in detail at demonstrated by the stability of rmsd, cluster counting, and loop L1, several interactions were identified as responsible for PCA analysis (Figures S3−S6). Neither the entire tertiary maintaining it in the crystallographic arrangement and structure nor secondary motives, including loop L1, presented preventing it from reaching an open form. Stable salt bridges significant variations. In general, the core structure of the and hydrogen bonds between residues Lys148 and Arg144 protein showed again low flexibility and very high stability, from helix-α A2 and residues Asp29 and 31 situated at the finding the only differences in loops L1 and L2 and the helix α- Nter of loop L1 (Figure 4a) were found during the entire 2 (Figures S4 and S6). These appeared more flexible in the apo simulation trajectory. Furthermore, 25−50% of the GaMD form than in the holo form, loop L1 being the one that showed trajectory, hydrogen bonds between polar residues of loop L1 the highest flexibility (average RMSF difference of 0.79 Å) (Lys38, Arg40, Ser60), and the backbone or side chains of the (Figures 5 and S7). It can be concluded, therefore, that the same loop were found. Overall, the hairpin H3 was found to be presence of heme restricts the movement of loop L1. Figure 4. Representation of interactions during GaMD simulations in HasAyp. Hydrogen bonds in blue and hydrophobic contacts in orange with main residues involved highlighted. (a) Interactions extracted from GaMD apo from HasAyp. (b) Interactions between heme and HasAyp extracted from GaMD holo. 17072 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article the neutral state with monoprotonation at Nδ will be described. The PCA and clustering analyses showed that the GaMDs converged by the end of the 800 ns (Figure S8). No major conformational changes were observed in the α- helix or β-sheets of the system. However, some variations were observed in the loop and hairpin regions. The RMSF analysis showed the highest values for loop L1, L2 and hairpin H1, H3 (Figures 6 and S9). Even though the flexibility of loop L1 is Figure 5. RMSF difference between apo and holo from GaMD of HasAyp. In blue, RMSF is higher in apo, and in red, RMSF is higher in holo. No significant RMSF differences are shown in yellow. Interestingly, the only region which was found to be slightly more flexible in the heme bound form was the hairpin H3 (average difference RMSF 0.13 Å). We ended our analysis comparing the network of interactions of the amino acids between the holo and apo forms. Although most hydrogen bond interactions are preserved in both systems, the holo form shows less frequent interactions inside loop L1 and between the loop and the protein than the apo form because these residues now interact Figure 6. RMSF of apo GaMD from HasAsm with relevant regions. with heme. Interestingly, despite the lack of coordination of the heme to any residue of L1, strong interactions appear quite high (average RMSF about 1.84 Å), no major movement between the prosthetic group and the loop, hence reducing its of loop L1 was observed in any of the replicas. Because loop L1 flexibility. Heme presents the main hydrogen bond and salt always remained close to the initial open conformation, bridge interactions between the propionates of the heme and simulations did not show the transition of loop L1 toward Arg144, Arg40, and His135 (Figure 4b). Since these more closed conformations. Therefore, nothing in the three interactions are maintained during the majority of the converged replica simulations of 800 ns suggests that the apo simulations, these may be responsible for maintaining loop form of HasAsm could naturally move toward a more closed L1 in a closed disposition even when heme is bound. In conformation without the presence of the heme cofactor. addition, the hydrophobic interactions found in the apo system To better understand why L1 is not able to escape its open involving Phe43, Phe83, and Tyr55 now interact directly with conformation and reach heme-bound like states, an analysis of the heme. However, these do not seem to affect the loop L1 the interactions during the GaMD simulations was carried out region. (Figure 7a). First, a network of hydrogen bonds can be Altogether, these calculations have revealed that the observed between loop L1 and three other regions: (1) its Nter hydrophobic residues Tyr55,75 and Phe43,83,50,77 of the part is fixed by interactions between Val30, Asn36, and Thr38 binding site are a platform for the binding of heme, while with hairpin H2 (Asn62 and Gln63); (2) central residues positively charged residues Arg144,40 and His135 stabilize Thr38 and Ser39 interact with β-sheets B5-6 (Ser99 or through salt bridges and made propionates face the solvent. Gln109), which keeps the most flexible part of the loop fixed; HasA from S. marcescens. In contrast to HasAyp, the apo and (3) its Cter has interaction of the backbone 40−44 with β- form of HasAsm has a completely different conformation of sheet B3 and the side chains of Ser42 with Ser58−59. All loop L1, presenting an open conformation in the heme-free interactions keep loop L1 fixed close to the initial position and structure (Figure 2a). We aimed at assessing the nature of the therefore in apo-like conformations. This is reinforced by a interactions occurring in both the holo and apo forms of network of hydrogen bonds in the heme-binding regions HasAsm and understanding the differences for heme uptake involving histidine 83, 128, and 133. Regarding the hydro- with respect to the HasAyp system. Following the same phobic interactions, in harpin H3, there are π−stacking protocol, GaMD simulations were carried out on the apo form interactions that stabilize close contacts with β-sheets B2−3 of HasAsm. Because His32, which is the heme axial ligand, and loop L1 (Figure 7b). Those involve mainly Tyr46 and could potentially be in different protonation states influencing Ala56 from β-sheet B2 that interact with hairpin H3 or β- the flexibility of loop L1, simulations were performed with sheets B5-6 with Pro105, Tyr106, and Leu98. However, the different protonation states of this residue. No significant most relevant hydrophobic interactions in the system are not differences were observed, and here, the system with His32 at related directly to loop L1 but are part of the heme-binding 17073 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article Figure 7. Representation of hydrogen bonds (HB) and hydrophobic (HP) interactions during GaMD of HasAsm in apo (a,b) and heme-bound before loop closing (c,d). Normalized difference between apo and holo forms of both hydrogen bonds (e) and hydrophobic interactions (f) is represented. See Figure S15 for Fe(II) analysis. 17074 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article region. Importantly, there are several π−stacking contacts has separated from the region of hairpin H3 and β-sheets B5− between the aromatic residues Phe45, Tyr75, Leu85, and 6. Only after 40 ns up to 700 ns, depending on the replica Tyr55, making a very robust network of interactions. simulation, the loop adopts a turn or small helix conformation As the GaMD simulations of the apo-HasAsm did not show and moves toward the heme-binding site and acquires a more any conformational changes of loop L1 consistent with the closed disposition (Figure 9a). heme-bound structure, we hypothesized that it may have been For the rest of the GaMDs, the loop oscillates between induced by the binding of the heme. We therefore pursued different conformations, with some of them consistent with with heme dockings followed by GaMD simulations. Similar to His32 facing the heme at distances consistent with what is observed in HasAyp, docking predicted good binding coordination to the metal, while in others, the nitrogen faces affinities (GoldScore values of ca. 68) and similar binding outside of the binding site. Despite the closing motion of the poses as in the X-ray structure of the heme-bound form, loop L1 appears once the heme has bound to the HasAsm apo showing a coordination bond between the heme and the form, a stable holo structure with His32 as the 6th ligand of the oxygen of the side chain of Tyr75. Surprisingly, the best iron is not maintained in the simulations. The distance of the docking solutions showed that the propionate groups were coordinating nitrogen of His32 to the iron tends to swing facing the inner part of the protein, and it was only after an between 3 and 12 Å depending on the replica with no stable equilibration with classical MD simulations of 10−20 ns that position of His32 in a coordination mode. This is due to the the heme rotated toward the external part of the binding force field conditions that only allow one coordination state at pocket and ended with a structure with excellent matching a time, and these simulations start with the sole Tyr75 bound with the experimental holo form. These results highlight that and the corresponding parameters of a pentacoordinated first binding of the heme to HasAsm could perfectly happen coordination sphere. Although some Fe(III) simulations show without the necessity of loop L1 to transit to the X-ray position shorter distances than in Fe(II) calculations, this occurs very and the His32 to coordinate the remaining axial site. These sporadically and indicates a limited impact on the net charge of simulations were followed by three replicas of 800 ns GaMD the metal, and the absence of an explicit Fe−His bond in the simulations. parametrization is the origin of the fluctuation. To further All the statistical indicators on the GaMD replica investigate this point, additional GaMD simulations were simulations (clustering, PCA, rmsd) revealed convergence performed. On one side, the holo X-ray structure (PDB code after 100−400 ns (Figures S10−S13). Several wells found in 1DKH) was submitted to GaMD simulations, with parameters the PCA analysis showed that the system visited distinct for the Fe−Tyr coordination and no coordination term for the conformations. In particular, the RMSF analysis showed that Fe−His bond. The histidine flies out of the heme binding the most flexible part of the protein corresponded to the loops pocket and fluctuates around 7.43 Å (going from 2.7 to 17.9); and hairpin regions, loop L1 being the most flexible one however, loop L1 remains mostly in a closed geometry (Figure (Figure S13). When comparing the RMSF with the apo form, S16). On the other hand, simulations starting from the once again a rigid core structure of the protein was observed snapshot of the GaMD of the HasAsm with the heme bound to with no changes in the tertiary structure, except for loop L1 the tyrosine and the smallest Fe−OTyr distance (2.9 Å) were (RMSF diff. of 4.1 Å) and the β-sheets B3 (RMSF diff. 0.98 Å) performed using a set of parameters for the hexacoordinated (Figures 8 and S14). Only loop L2 is more flexible in the apo metal with axial His-Fe-Tyr configuration. This simulation rapidly reaches a loop L1 conformation close to the experimental holo form (rmsd of 0.78 Å), and a stable His− Fe bond is observed during the simulation (Figure S16). This shows that reaching the final hexacoordinated structure requires the His coordination to be properly modeled. Still, the major determinants of the loop closure motion are the fluctuations of the protein once the heme has bound its cavity. PCA analysis of fragment trajectory where loop L1 closes with respect to the Fe−His32 distance reveals how PC1 involves the movement of loop L1 (Figure 9b), and reweighted PMF calculations in front of the PCA1 distance Fe−His32 shows that the barrier of the system to transit from one state to another is less than 8 kcal/mol (Figures 9c and S17 and S18). Analysis of interactions reveals how the binding of the heme clearly induces some changes in the network of hydrophobic and stacking interactions (Figure 7c,d). Basically, the binding Figure 8. RMSF difference between apo in blue and holo in orange disrupts all the previously mentioned π−stacking and hydro- extracted from apo and holo GaMD of HasAsm. phobic interactions between aromatic residues of the heme binding site (Tyr75, Tyr137, His83, Tyr55, and Phe45). form. This result sustains that the binding of heme, with Instead of interacting between them, these residues now extremity ends at the opposite side of the protein and in its interact with the heme. Furthermore, the interaction between pocket, induces a conformational change at loop L1, initiating Tyr46 and Tyr106 decreases, allowing much more flexibility of its closing mechanism. hairpin H3 that separates from β-sheets B2−3 and from loop Visual inspection of the trajectory of the three GaMD L1. This affects loop L1 and its hydrogen bond interactions; replicas shows a similar mechanism for loop L1 closing. during the first half of the simulations, all hydrogen bonds Detailed analysis shows that at the beginning of the GaMD between L1 and hairpin H3 drastically decrease, making it simulation, loop L1 tends to acquire an open conformation; it more flexible. Loop L1 can make more hydrogen bonds with β- 17075 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article Figure 9. GaMD simulation of HasAsm with heme−Fe(III) bound. (a) Frames of GaMD showing the loop L1 closing process colored according to the distance between Fe and His32 during GaMD. Graphic of distance Fe−His32 represented. (b) Cartesian coordinate PCA analysis colored according to the distance between Fe and His32 during GaMD. (c) Reweighted PMF calculations in front of PCA1 and distance Fe−His32. (b,c) are obtained using the fragment of the trajectory in which the loop is closing. sheets B5−6 and hairpin H2 (mainly with Thr60, Asn62, and resulting from the movement of loop L1, as described in Gln63). As mentioned before, the binding of the heme HasAsm experimental structures, is never observed despite decreases the π−stacking interactions between Phe45 from GaMD simulations allowing extensive conformational sam- loop L1 and the residues of the heme binding site; therefore, pling. Calculations show that this is due to specific salt bridges Phe45 can interact with residues Leu50 and Leu77 from and hydrogen bonds between negatively charged residues of hairpin H1. This change of interactions induces conforma- loop L1 (mostly Asp30 and Asp31) with positive residues from tional changes in the region of hairpin H1 and Cter of loop L1. helix from the core (the most important ones being Lys148 This could be the cause of the appearance of a turn or semi- and Arg144). The binding of the heme only slightly impacts helix in disposition on loop L1. At this point, the previously some hydrophobic interactions but not sufficiently to disrupt mentioned interactions with β-sheet B3 and H2 are substituted the geometry of loop L1. It can be concluded that the by a series of hydrogen bond interactions inside the same loop mechanism of heme-binding would correspond to a very light L1. This change also causes the separation of loop L1 from β- conformational selection. This part of the study shows that the sheets B5-6, and loop L1 starts to move toward the heme- HasAyp hemophore presents a well pre-organized geometry of binding site, while at the same time, α-helix A1 moves toward the receptor for heme binding as a classified transient heme- α-helix A2. From this point, where the loop L1 is in closed binding protein. disposition, there are a series of hydrogen bond interactions Simulations on HasAsm also show that both apo and holo between loop L1 and α-helix A2 that maintain it in a close structures are very stable. Interestingly, no structural rearrange- conformation. The main residues involved are Asn31 and ment of loop L1 is observed in any of them. Simulations on the Val30 from loop L1 with Ser141 from α-helix A2. If we apo form show that the binding site of the heme is well pre- compare the frequency of all these interactions with the apo organized and that binding should occur naturally. Moreover, form and normalize it, these tendencies can be clearly observed the motion of loop L1 and the coordination of His32 have no in Figure 7e,f. impact on the insertion of the heme in its cavity or the In-Depth Comparison between HasA from Y. pestis coordination of the oxygen of the side chain of Tyr75 to the and S. marcescens. In this study, the combination of GaMD heme. As a crucial result of this study, the movements of loop and docking highlights differing patterns of heme-binding L1 are only observed after the heme process takes place. mechanisms between the hemophores from Y. pestis and S. Analysis of contacts shows that the binding of heme disrupts marcescens. the network of hydrophobic and π−stacking interactions at the Simulations on the HasAyp clearly demonstrate that the apo heme-binding site and that this information is propagated to form remains in a geometry similar to that of the holo form. the other extreme of the protein. This change induces a The transition between a close and open conformation conformational change in hairpin H3 that acquires more 17076 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article flexibility and loses interactions with residues of loop L1. This acquires far wider flexibility. This ultimately leads to also leads to a disruption of the interactions between the CTer conformations where the coordinating histidine reaches the part of loop L1; both hydrophobic and hydrophobic contacts heme-binding pocket and forms an additional bond with the with hairpin H2 and β-sheets B5−6 are broken, which causes iron. loop L1 to be more flexible. Once L1 has been liberated from This study shows interesting information on heme recruit- this strong interaction, it becomes free to reorganize and move ment by different families of hemophores. In all cases, it would up to the heme-binding site. This ultimately leads to the appear that the good pre-organization of the heme pocket for correct organization of the loop so that His32 could find the the binding of the cofactor is consistent with transient heme and coordinate. mechanisms. However, for the S. marcescens, there is a This mechanism sustains that the general transfer of subsequent induced fit effect that allows the system to reach structural information along the structure of the protein is its final holo that ultimately leads to an additional coordination similar to those observed on hemophore HasAp using targeted bond. How the additional step observed in HasAsm could be dynamics. In particular, it suggested that the interaction an advantage or inconvenient for the organism from the between residues from helix A2 with heme induces a tilting evolutive point of view remains unclear. However, this study motion that perturbs interactions with helix A1 and loop L1, completes our knowledge on heme-binding complexity and initiating its closing movement. Compared to previous studies, how long-range interactions could be crucial for defining the we are here able to simulate the transition from apo to holo heme-bound geometry and the mechanism of acquisition. In without forcing or constraining the system. Long GaMD an age of increasing interest in the insertion of metal-based simulations of 800 ns simulations have been performed, which compounds into protein scaffold, this study highlights the have allowed all systems to converge and observe significant importance of weighing pre-organization and induced fit conformational changes. Furthermore, performing three effects, understanding their origin, and predicting their replicas for each case assured that these events are not casual. magnitude. The conformational changes of L1 observed in the experimental structures of the apo and holo forms are, in reality, a subsequent process where long-range interactions are ■ * ASSOCIATED CONTENT sı Supporting Information involved as a transmitted signal from the heme-binding pocket The Supporting Information is available free of charge at and can ultimately lead to an additional coordination with https://pubs.acs.org/doi/10.1021/acs.inorgchem.2c02193. respect to other hemophores. ■ CONCLUSIONS In this study, we investigated the heme-binding mechanism of Additional conformational and explorational analysis of all GaMD simulations, interaction analysis, and force field parameters for heme groups (PDF) two hemophores, one from Y. pestis and the other from S. marcescens. Despite similar fold and biochemical functions, the former presents little conformation changes between apo and holo forms, while the latter presents an important conforma- ■ AUTHOR INFORMATION Corresponding Author tional change characterized by a large movement of one loop, Jean-Didier Maréchal − Insilichem, Departament de Química, which ends in the formation of an additional coordination Universitat Autònoma de Barcelona, 08193 Bellaterra, bound between an axial histidine and the metal. To decode the Spain; orcid.org/0000-0002-8344-9043; origin of such differences, we combine molecular dynamics Email: jeandidier.marechal@uab.cat experiments with protein−ligand dockings. We take advantage Authors of the extensive conformational sampling that Gaussian accelerated molecular dynamics offer to ascertain if the Laura Tiessler-Sala − Insilichem, Departament de Química, receptors could easily reach sub-states compatible for heme- Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain binding and unbinding processes as well as docking techniques Giuseppe Sciortino − Insilichem, Departament de Química, that have been optimized to deal with metalloligands. Universitat Autònoma de Barcelona, 08193 Bellaterra, From this study, it can be concluded that in all cases, the apo Spain; Institute of Chemical Research of Catalonia (ICIQ), forms are very stable. Loop L1 is kept fixed due to a strong The Barcelona Institute of Science and Technology, 43007 network of hydrogen bonds with amino acids of surrounding Tarragona, Spain; orcid.org/0000-0001-9657-1788 loops and β-sheets. These unbound forms are also well pre- Lur Alonso-Cotchico − Insilichem, Departament de Química, organized for the binding of the heme. In both cases, the Universitat Autònoma de Barcelona, 08193 Bellaterra, binding of the heme occurs without any previous major Spain; Zymvol Biomodeling, 08018 Barcelona, Spain conformational changes including the loop transition. The apo Laura Masgrau − Insilichem, Departament de Química, is therefore very well pre-organized for heme binding. Universitat Autònoma de Barcelona, 08193 Bellaterra, The main differences between both species come from the Spain; Zymvol Biomodeling, 08018 Barcelona, Spain; cascade of molecular events happening after the heme binds orcid.org/0000-0003-4495-508X and the conformational changes associated to them. While very Agustí Lledós − Insilichem, Departament de Química, little perturbation of the overall map of contacts and Universitat Autònoma de Barcelona, 08193 Bellaterra, interactions happens in Y. pestis, a series of changes take Spain; orcid.org/0000-0001-7909-422X place in S. marcescens. The mechanism starts with modifica- Complete contact information is available at: tions of hydrophobic contacts in the heme-binding pocket https://pubs.acs.org/10.1021/acs.inorgchem.2c02193 followed by rearrangements of hairpin and loop contacts on the opposite side of the protein that end up modifying the Notes patch of the interactions that loop L1 is part of. As such, L1 The authors declare no competing financial interest. 17077 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article ■ ACKNOWLEDGMENTS L.T.-S., A.L., G.S., and J.-D.M. thank the Spanish MINECO (20) Arnoux, P.; Haser, R.; Izadi, N.; Lecroisey, A.; Delepierre, M.; Wandersman, C.; Czjzek, M. The Crystal Structure of HasA, a Hemophore Secreted by Serratia Marcescens. Nat. Struct. Biol. 1999, (project PID-2020-116861GB-I00). L.T.-S. thanks the Spanish 6, 516−520. Ministerio de Ciencia, Innovación y Universidades (grant (21) Kumar, R.; Lovell, S.; Matsumura, H.; Battaile, K. P.; Moënne- FPU18/05895) for the financial support. L.M. acknowledges Loccoz, P.; Rivera, M. The Hemophore HasA from Yersinia Pestis funding from the Ministerio de Ciencia, Innovación y (HasAyp) Coordinates Hemin with a Single Residue, Tyr75, and with Universidades (grant PGC2018-098592-B-100) and the Minimal Conformational Change. Biochemistry 2013, 52, 2705−2707. Universitat Autònoma de Barcelona Talent Program. (22) Arnoux, P.; Haser, R.; Izadi-Pruneyre, N.; Lecroisey, A.; Czjzek, M. Functional Aspects of the Heme Bound Hemophore HasA by ■ REFERENCES (1) Andreini, C.; Bertini, I.; Cavallaro, G.; Holliday, G. L.; Thornton, Structural Analysis of Various Crystal Forms. Proteins Struct. Funct. Genet. 2000, 41, 202−210. (23) Jepkorir, G.; Rodríguez, J. C.; Rui, H.; Im, W.; Lovell, S.; J. M. Metal Ions in Biological Catalysis: From Enzyme Databases to Battaile, K. P.; Alontaga, A. Y.; Yukl, E. T.; Moënne-Loccoz, P.; General Principles. J. Biol. Inorg Chem. 2008, 13, 1205−1218. Rivera, M. Structural, NMR Spectroscopic, and Computational (2) Ringenberg, M. R.; Ward, T. R. Merging the Best of Two Investigation of Hemin Loading in the Hemophore Hasap from Worlds: Artificial Metalloenzymes for Enantioselective Catalysis. Pseudomonas Aeruginosa. J. Am. Chem. Soc. 2010, 132, 9857−9872. Chem. Commun. 2011, 47, 8470−8476. (24) Kumar, R.; Qi, Y.; Matsumura, H.; Lovell, S.; Yao, H.; Battaile, (3) Kariyawasam, K.; Di Meo, T.; Hammerer, F.; Valerio-Lepiniec, M.; Sciortino, G.; Maréchal, J. D.; Minard, P.; Mahy, J. P.; Urvoas, A.; K. P.; Im, W.; Moënne-Loccoz, P.; Rivera, M. Replacing Arginine 33 Ricoux, R. An Artificial Hemoprotein with Inducible Peroxidase- and for Alanine in the Hemophore HasA from Pseudomonas Aeruginosa Monooxygenase-Like Activities. Chem.�Eur. J. 2020, 26, 14929− Causes Closure of the H32 Loop in the Apo-Protein. Biochemistry 14937. 2016, 55, 2622−2631. (4) Prier, C. K.; Zhang, R. K.; Buller, A. R.; Brinkmann-Chen, S.; (25) Exner, T. E.; Becker, S.; Becker, S.; Boniface-Guiraud, A.; Arnold, F. H. Enantioselective, Intermolecular Benzylic C−H Delepelaire, P.; Diederichs, K.; Welte, W. Binding of HasA by Its Amination Catalysed by an Engineered Iron-Haem Enzyme. Nat. Transmembrane Receptor HasR Follows a Conformational Funnel Chem. 2017, 9, 629−634. Mechanism. Eur. Biophys. J. 2020, 49, 39−57. (5) Wilson, C. J.; Apiyo, D.; Wittung-Stafshede, P. Role of Cofactors (26) Schlitter, J.; Engels, M.; Krüger, P. Targeted Molecular in Metalloprotein Folding. Q. Rev. Biophys. 2004, 37, 285−314. Dynamics: A New Approach for Searching Pathways of Conforma- (6) Poulos, T. L. Heme Enzyme Structure and Function. Chem. Rev. tional Transitions. J. Mol. Graph. 1994, 12, 84−89. 2014, 114, 3919−3962. (27) Burley, S. K.; Bhikadiya, C.; Bi, C.; Bittrich, S.; Chen, L.; (7) Zhang, L. Heme Biology: The Secret Life of Heme in Regulating Crichlow, G. V.; Christie, C. H.; Dalenberg, K.; Di Costanzo, L.; Diverse Biological Processes; World Scientific: Singapore, 2011. Duarte, J. M.; Dutta, S.; Feng, Z.; Ganesan, S.; Goodsell, D. S.; Ghosh, (8) Ma, B.; Kumar, S.; Tsai, C. J.; Nussinov, R. Folding Funnels and S.; Green, R. K.; Guranović, V.; Guzenko, D.; Hudson, B. P.; Lawson, Binding Mechanisms. Protein Eng. 1999, 12, 713−720. C. L.; Liang, Y.; Lowe, R.; Namkoong, H.; Peisach, E.; Persikova, I.; (9) Koshland, D. E. Application of a Theory of Enzyme Specificity to Randle, C.; Rose, A.; Rose, Y.; Sali, A.; Segura, J.; Sekharan, M.; Shao, Protein Synthesis. Proc. Natl. Acad. Sci. U.S.A. 1958, 44, 98−104. C.; Tao, Y. P.; Voigt, M.; Westbrook, J. D.; Young, J. Y.; Zardecki, C.; (10) Vogt, A. D.; Di Cera, E. Conformational Selection Is a Zhuravleva, M. RCSB Protein Data Bank: Powerful New Tools for Dominant Mechanism of Ligand Binding. Biochemistry 2013, 52, Exploring 3D Structures of Biological Macromolecules for Basic and 5723−5729. Applied Research and Education in Fundamental Biology, Bio- (11) Vogt, A. D.; Pozzi, N.; Chen, Z.; Di Cera, E. Essential Role of medicine, Biotechnology, Bioengineering and Energy Sciences. Nucleic Conformational Selection in Ligand Binding. Biophys. Chem. 2014, Acids Res. 2021, 49, D437−D451. 186, 13−21. (28) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; (12) Guengerich, F. P.; Wilkey, C. J.; Phan, T. T. N. Human Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera?A Cytochrome P450 Enzymes Bind Drugs and Other Substrates Mainly visualization system for exploratory research and analysis. J. Comput. through Conformational-Selection Modes. J. Biol. Chem. 2019, 294, Chem. 2004, 25, 1605−1612. 10928−10941. (29) Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: (13) Li, T.; Bonkovsky, H. L.; Guo, J. T. Structural Analysis of Heme Automating PK Prediction and the Preparation of Biomolecular Proteins: Implication for Design and Prediction. 2010 IEEE Structures for Atomistic Molecular Modeling and Simulations. Nucleic International Conference on Bioinformatics and Biomedicine Workshops, Acids Res. 2012, 40, W537−W541. 2010; pp 834−835. (30) Verdonk, M. L.; Cole, J. C.; Hartshorn, M. J.; Murray, C. W.; (14) Eliezer, D.; Wright, P. E. Is Apomyoglobin a Molten Globule? Taylor, R. D. Improved Protein-Ligand Docking Using GOLD. Structural Characterization by NMR. J. Mol. Biol. 1996, 263, 531− Proteins 2003, 52, 609−623. 538. (31) Sciortino, G.; Rodríguez-Guerra Pedregal, J.; Lledós, A.; (15) Falzone, C. J.; Mayer, M. R.; Whiteman, E. L.; Moore, C. D.; Garribba, E.; Maréchal, J.-D. Prediction of the Interaction of Metallic Lecomte, J. T. J. Design Challenges for Hemoproteins: The Solution Moieties with Proteins: An Update for Protein-Ligand Docking Structure of Apocytochrome B5. Biochemistry 1996, 35, 6519−6526. Techniques. J. Comput. Chem. 2018, 39, 42−51. (16) Falzone, C. J.; Wang, Y.; Vu, B. C.; Scott, N. L.; Bhattacharya, (32) Rodríguez-Guerra Pedregal, J.; Sciortino, G.; Guasp, J.; S.; Lecomte, J. T. J. Structural and Dynamic Perturbations Induced by Municoy, M.; Maréchal, J.-D. GaudiMM: A Modular Multi-Objective Heme Binding in Cytochrome b5. Biochemistry 2001, 40, 4879−4891. Platform for Molecular Modeling. J. Comput. Chem. 2017, 38, 2118− (17) Feng, Y.; Sligar, S. G.; Wand, A. J. Solution Structure of 2126. Apocytochrome B562. Nat. Struct. Biol. 1994, 1, 30−35. (33) Miao, Y.; Feher, V. A.; McCammon, J. A. Gaussian Accelerated (18) Arnesano, F.; Banci, L.; Bertini, I.; Faraone-Mennella, J.; Molecular Dynamics: Unconstrained Enhanced Sampling and Free Rosato, A.; Barker, P. D.; Fersht, A. R. The Solution Structure of Energy Calculation. J. Chem. Theory Comput. 2015, 11, 3584−3595. Oxidized Escherichia coli Cytochrome b562. Biochemistry 1999, 38, (34) Case, D. A.; Ben-Shalom, I. Y.; Brozell, S. R.; Cerutti, D. S.; 8657−8670. Cheatham, T. E., III; Cruzeiro, V. W. D.; Darden, T. A.; Duke, R. E.; (19) Schneider, S.; Sharp, K. H.; Barker, P. D.; Paoli, M. An Induced Ghoreishi, D.; Gilson, M. K.; Gohlke, H., Goetz, A. W.; Greene, D.; Fit Conformational Change Underlies the Binding Mechanism of the Harris, R.; Homeyer, N.; Izadi, S.; Kovalenko, A.; Kurtzman, T.; Lee, Heme Transport Proteobacteria-Protein HemS. J. Biol. Chem. 2006, T. S.; LeGra, S.; York, D. M.; Kollman, P. A. Amber 2018; University 281, 32606−32610. of California: San Francisco, 2018. 17078 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079 Inorganic Chemistry pubs.acs.org/IC Article (35) Cornell, W. D.; Cieplak, P.; Gould, I. R.; Bayly, C. I.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (36) Li, P.; Merz, K. M. MCPB.py: A Python Based Metal Center Parameter Builder. J. Chem. Inf. Model. 2016, 56, 599−604. (37) Seminario, J. M. Calculation of Intramolecular Force Fields from Second-Derivative Tensors. Int. J. Quantum Chem. 1996, 60, 1271−1277. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ort, J. V. Gaussian 09, Revision A.02, 2016. (39) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (40) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (41) Ehlers, A. W.; Böhme, M.; Dapprich, S.; Gobbi, A.; Höllwarth, A.; Jonas, V.; Köhler, K. F.; Stegmann, R.; Veldkamp, A.; Frenking, G. A set of f-polarization functions for pseudo-potential basis sets of the transition metals ScCu, YAg and LaAu. Chem. Phys. Lett. 1993, 208, 111−114. (42) McGibbon, R. T.; Beauchamp, K. A.; Harrigan, M. P.; Klein, C.; Swails, J. M.; Hernández, C. X.; Schwantes, C. R.; Wang, L. P.; Lane, T. J.; Pande, V. S. MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. Biophys. J. 2015, 109, 1528−1532. (43) Alonso-Cotchico, L.; Rodríguez-Guerra Pedregal, J.; Lledós, A.; Maréchal, J.-D. The Effect of Cofactor Binding on the Conforma- tional Plasticity of the Biological Receptors in Artificial Metal- loenzymes: The Case Study of LmrR. Front. Chem. 2019, 7, 211. (44) Venkatakrishnan, A. J.; Fonseca, R.; Ma, A. K.; Hollingsworth, S. A.; Chemparathy, A.; Hilger, D.; Kooistra, A. J.; Ahmari, R.; Madan, M. B.; Kobilka, B. K.; Dror, R. O. Uncovering Patterns of Atomic Interactions in Static and Dynamic Structures of Proteins, 2019. bioRxiv:10.1101/840694. (45) Miao, Y.; Sinko, W.; Pierce, L.; Bucher, D.; Walker, R. C.; McCammon, J. A. Improved Reweighting of Accelerated Molecular Dynamics Simulations for Free Energy Calculation. J. Chem. Theory Comput. 2014, 10, 2677−2689. 17079 https://doi.org/10.1021/acs.inorgchem.2c02193 Inorg. Chem. 2022, 61, 17068−17079