An In Silico Analysis of the Binding Modes and Binding Affinities of Small Molecule Modulators of PDZ-Peptide Interactions

Inhibitors of PDZ-peptide interactions have important implications in a variety of biological processes including treatment of cancer and Parkinson’s disease. Even though experimental studies have reported characterization of peptidomimetic inhibitors of PDZ-peptide interactions, the binding modes for most of them have not been characterized by structural studies. In this study we have attempted to understand the structural basis of the small molecule-PDZ interactions by in silico analysis of the binding modes and binding affinities of a set of 38 small molecules with known K_i or K_d values for PDZ2 and PDZ3 domains of PSD-95 protein. These two PDZ domains show differential selectivity for these compounds despite having a high degree of sequence similarity and almost identical peptide binding pockets. Optimum binding modes for these ligands for PDZ2 and PDZ3 domains were identified by using a novel combination of semi-flexible docking and explicit solvent molecular dynamics (MD) simulations. Analysis of the binding modes revealed most of the peptidomimectic ligands which had high K_i or K_d moved away from the peptide binding pocket, while ligands with high binding affinities remained in the peptide binding pocket. The differential specificities of the PDZ2 and PDZ3 domains primarily arise from differences in the conformation of the loop connecting βB and βC strands, because this loop interacts with the N-terminal chemical moieties of the ligands. We have also computed the MM/PBSA binding free energy values for these 38 compounds with both the PDZ domains from multiple 5 ns MD trajectories on each complex i.e. a total of 228 MD trajectories of 5 ns length each. Interestingly, computational binding free energies show good agreement with experimental binding free energies with a correlation coefficient of approximately 0.6. Thus our study demonstrates that combined use of docking and MD simulations can help in identification of potent inhibitors of PDZ-peptide complexes.     

High‐energy water sites determine peptide binding affinity and specificity of PDZ domains

PDZ domains have well known binding preferences for distinct C‐terminal peptide motifs. For most PDZ domains, these motifs are of the form [S/T]‐W‐[I/L/V]. Although the preference for S/T has been explained by a specific hydrogen bond interaction with a histidine in the PDZ domain and the (I/L/V) is buried in a hydrophobic pocket, the mechanism for Trp specificity at the second to last position has thus far remained unknown. Here, we apply a method to compute the free energies of explicit water molecules and predict that potency gained by Trp binding is due to a favorable release of high‐energy water molecules into bulk. The affinities of a series of peptides for both wild‐type and mutant forms of the PDZ domain of Erbin correlate very well with the computed free energy of binding of displaced waters, suggesting a direct relationship between water displacement and peptide affinity. Finally, we show a correlation between the magnitude of the displaced water free energy and the degree of Trp‐sensitivity among subtypes of the HTRA PDZ family, indicating a water‐mediated mechanism for specificity of peptide binding.     

Signaling pathways of PDZ2 domain: a molecular dynamics interaction correlation analysis

PDZ domains are found in many signaling proteins. One of their functions is to provide scaffolds for forming membrane-associated protein complexes by binding to the carboxyl termini of their partners. PDZ domains are thought also to play a signal transduction role by propagating the information that binding has occurred to remote sites. In this study, a molecular dynamics (MD) simulation-based approach, referred to as an interaction correlation analysis, is applied to the PDZ2 domain to identify the possible signal transduction pathways. A residue correlation matrix is constructed from the interaction energy correlations between all residue pairs obtained from the MD simulations. Two continuous interaction pathways, starting at the ligand binding pocket, are identified by a hierarchical clustering analysis of the residue correlation matrix. One pathway is mainly localized at the N-terminal side of helix alpha1 and the adjacent C-terminus of loop beta1-beta2. The other pathway is perpendicular to the central beta-sheet and extends toward the side of PDZ2 domain opposite to the ligand binding pocket. The results complement previous studies based on multiple sequence analysis, NMR, and MD simulations. Importantly, they reveal the energetic origin of the long-range coupling. The PDZ2 results, as well as the earlier rhodopsin analysis, show that the interaction correlation analysis is a robust approach for determining pathways of intramolecular signal transduction.     

Beyond the Binding Site: The Role of the β2 – β3 Loop and Extra-Domain Structures in PDZ Domains

A general paradigm to understand protein function is to look at properties of isolated well conserved domains, such as SH3 or PDZ domains. While common features of domain families are well understood, the role of subtle differences among members of these families is less clear. Here, molecular dynamics simulations indicate that the binding mechanism in PSD95-PDZ3 is critically regulated via interactions outside the canonical binding site, involving both the poorly conserved loop and an extra-domain helix. Using the CRIPT peptide as a prototypical ligand, our simulations suggest that a network of salt-bridges between the ligand and this loop is necessary for binding. These contacts interconvert between each other on a time scale of a few tens of nanoseconds, making them elusive to X-ray crystallography. The loop is stabilized by an extra-domain helix. The latter influences the global dynamics of the domain, considerably increasing binding affinity. We found that two key contacts between the helix and the domain, one involving the loop, provide an atomistic interpretation of the increased affinity. Our analysis indicates that both extra-domain segments and loosely conserved regions play critical roles in PDZ binding affinity and specificity.     

Structural basis for NHERF1 PDZ domain binding

The Na(+)/H(+) exchange regulatory factor-1 (NHERF1) is a scaffolding protein that possesses two tandem PDZ domains and a carboxy-terminal ezrin-binding domain (EBD). The parathyroid hormone receptor (PTHR), type II sodium-dependent phosphate cotransporter (Npt2a), and β2-adrenergic receptor (β2-AR), through their respective carboxy-terminal PDZ-recognition motifs, individually interact with NHERF1 forming a complex with one of the PDZ domains. In the basal state, NHERF1 adopts a self-inhibited conformation, in which its carboxy-terminal PDZ ligand interacts with PDZ2. We applied molecular dynamics (MD) simulations to uncover the structural and biochemical basis for the binding selectivity of NHERF1 PDZ domains. PDZ1 uniquely forms several contacts not present in PDZ2 that further stabilize PDZ1 interactions with target ligands. The binding free energy (ΔG) of PDZ1 and PDZ2 with the carboxy-terminal, five-amino acid residues that form the PDZ-recognition motif of PTHR, Npt2a, and β2-AR was calculated and compared with the calculated ΔG for the self-association of NHERF1. The results suggest that the interaction of the PTHR, β2-adrenergic, and Npt2a involves competition between NHERF1 PDZ domains and the target proteins. The binding of PDZ2 with PTHR may also compete with the self-inhibited conformation of NHERF1, thereby contributing to the stabilization of an active NHERF1 conformation.     

Beyond the binding site: the role of the β₂-β₃ loop and extra-domain structures in PDZ domains

A general paradigm to understand protein function is to look at properties of isolated well conserved domains, such as SH3 or PDZ domains. While common features of domain families are well understood, the role of subtle differences among members of these families is less clear. Here, molecular dynamics simulations indicate that the binding mechanism in PSD95-PDZ3 is critically regulated via interactions outside the canonical binding site, involving both the poorly conserved β?-β? loop and an extra-domain helix. Using the CRIPT peptide as a prototypical ligand, our simulations suggest that a network of salt-bridges between the ligand and this loop is necessary for binding. These contacts interconvert between each other on a time scale of a few tens of nanoseconds, making them elusive to X-ray crystallography. The loop is stabilized by an extra-domain helix. The latter influences the global dynamics of the domain, considerably increasing binding affinity. We found that two key contacts between the helix and the domain, one involving the β?-β? loop, provide an atomistic interpretation of the increased affinity. Our analysis indicates that both extra-domain segments and loosely conserved regions play critical roles in PDZ binding affinity and specificity.     

Characterization of PDZ domain‐peptide interaction interface based on energetic patterns

PDZ domain is one of the abundant modular domains that recognize short peptide sequences to mediate protein–protein interactions. To decipher the binding specificity of PDZ domain, we analyzed the interactions between 11 mouse PDZ domains and 217 peptides using a method called MIECSVM, which energetically characterizes the domain‐peptide interaction using molecular interaction energy components (MIECs) and predicts binding specificity using support vector machine (SVM). Cross‐validation and leave‐one‐domain‐out test showed that the MIEC‐SVM using all 44 PDZ‐peptide residue pairs at the interaction interface outperformed the sequence‐based methods in the literature. A further feature (residue pair) selection procedure illustrated that 16 residue pairs were uninformative to the binding specificity, even though they contributed significantly (～50%) to the binding energy. If only using the 28 informative residue pairs, the performance of the MIEC‐SVM on predicting the PDZ binding specificity was significantly improved. This analysis suggests that the informative and uninformative residue interactions between the PDZ domain and the peptide may represent those contributing to binding specificity and affinity, respectively. We performed additional structural and energetic analyses to shed light on understanding how the PDZ‐peptide recognition is established. The success of the MIEC‐SVM method on PDZ domains in this study and SH3 domains in our previous studies illustrates its generality on characterizing protein‐ peptide interactions and understanding protein recognition from a structural and energetic viewpoint.     

Structural and functional analysis of the PDZ domains of human HtrA1 and HtrA3

High-temperature requirement A (HtrA) and its homologs contain a serine protease domain followed by one or two PDZ domains. Bacterial HtrA proteins and the mitochondrial protein HtrA2/Omi maintain cell function by acting as both molecular chaperones and proteases to manage misfolded proteins. The biological roles of the mammalian family members HtrA1 and HtrA3 are less clear. We report a detailed structural and functional analysis of the PDZ domains of human HtrA1 and HtrA3 using peptide libraries and affinity assays to define specificity, structural studies to view the molecular details of ligand recognition, and alanine scanning mutagenesis to investigate the energetic contributions of individual residues to ligand binding. In common with HtrA2/Omi, we show that the PDZ domains of HtrA1 and HtrA3 recognize hydrophobic polypeptides, and while C-terminal sequences are preferred, internal sequences are also recognized. However, the details of the interactions differ, as different domains rely on interactions with different residues within the ligand to achieve high affinity binding. The results suggest that mammalian HtrA PDZ domains interact with a broad range of hydrophobic binding partners. This promiscuous specificity resembles that of bacterial HtrA family members and suggests a similar function for recognizing misfolded polypeptides with exposed hydrophobic sequences. Our results support a common activation mechanism for the HtrA family, whereby hydrophobic peptides bind to the PDZ domain and induce conformational changes that activate the protease. Such a mechanism is well suited to proteases evolved for the recognition and degradation of misfolded proteins.     

Proteome‐wide analysis of phospho‐regulated PDZ domain interactions

A key function of reversible protein phosphorylation is to regulate protein–protein interactions, many of which involve short linear motifs (3–12 amino acids). Motif‐based interactions are difficult to capture because of their often low‐to‐moderate affinities. Here, we describe phosphomimetic proteomic peptide‐phage display, a powerful method for simultaneously finding motif‐based interaction and pinpointing phosphorylation switches. We computationally designed an oligonucleotide library encoding human C‐terminal peptides containing known or predicted Ser/Thr phosphosites and phosphomimetic variants thereof. We incorporated these oligonucleotides into a phage library and screened the PDZ (PSD‐95/Dlg/ZO‐1) domains of Scribble and DLG1 for interactions potentially enabled or disabled by ligand phosphorylation. We identified known and novel binders and characterized selected interactions through microscale thermophoresis, isothermal titration calorimetry, and NMR. We uncover site‐specific phospho‐regulation of PDZ domain interactions, provide a structural framework for how PDZ domains accomplish phosphopeptide binding, and discuss ligand phosphorylation as a switching mechanism of PDZ domain interactions. The approach is readily scalable and can be used to explore the potential phospho‐regulation of motif‐based interactions on a large scale.     

Binding free energy landscape of domain-peptide interactions

Peptide recognition domains (PRDs) are ubiquitous protein domains which mediate large numbers of protein interactions in the cell. How these PRDs are able to recognize peptide sequences in a rapid and specific manner is incompletely understood. We explore the peptide binding process of PDZ domains, a large PRD family, from an equilibrium perspective using an all-atom Monte Carlo (MC) approach. Our focus is two different PDZ domains representing two major PDZ classes, I and II. For both domains, a binding free energy surface with a strong bias toward the native bound state is found. Moreover, both domains exhibit a binding process in which the peptides are mostly either bound at the PDZ binding pocket or else interact little with the domain surface. Consistent with this, various binding observables show a temperature dependence well described by a simple two-state model. We also find important differences in the details between the two domains. While both domains exhibit well-defined binding free energy barriers, the class I barrier is significantly weaker than the one for class II. To probe this issue further, we apply our method to a PDZ domain with dual specificity for class I and II peptides, and find an analogous difference in their binding free energy barriers. Lastly, we perform a large number of fixed-temperature MC kinetics trajectories under binding conditions. These trajectories reveal significantly slower binding dynamics for the class II domain relative to class I. Our combined results are consistent with a binding mechanism in which the peptide C terminal residue binds in an initial, rate-limiting step.     

Comprehensive analysis of all evolutionary paths between two divergent PDZ domain specificities

To understand the molecular evolution of functional diversity in protein families, we comprehensively investigated the consequences of all possible mutation combinations separating two peptide‐binding domains with highly divergent specificities. We analyzed the Erbin PDZ domain (Erbin‐PDZ), which exhibits canonical type I specificity, and a synthetic Erbin‐PDZ variant (E‐14) that differs at six positions and exhibits an atypical specificity that closely resembles that of the natural Pdlim4 PDZ domain (Pdlim4‐PDZ). We constructed a panel of 64 PDZ domains covering all possible transitions between Erbin‐PDZ and E‐14 (i.e., the panel contained variants with all possible combinations of either the Erbin‐PDZ or E‐14 sequence at the six differing positions). We assessed the specificity profiles of the 64 PDZ domains using a C‐terminal phage‐displayed peptide library containing all possible genetically encoded heptapeptides. The specificity profiles clustered into six distinct groups, showing that intermediate domains can be nodes for the evolution of divergent functions. Remarkably, three substitutions were sufficient to convert the specificity of Erbin‐PDZ to that of Pdlim4‐PDZ, whereas Pdlim4‐PDZ contains 71 differences relative to Erbin‐PDZ. X‐ray crystallography revealed the structural basis for specificity transition: a single substitution in the center of the binding site, supported by contributions from auxiliary substitutions, altered the main chain conformation of the peptide ligand to resemble that of ligands bound to Pdlim4‐PDZ. Our results show that a very small set of mutations can dramatically alter protein specificity, and these findings support the hypothesis whereby complex protein functions evolve by gene duplication followed by cumulative mutations.     

Peptide selectivity between the PDZ domains of human pregnancy‐related serine proteases (HtrA1, HtrA2, HtrA3, and HtrA4) can be reshaped by different halogen probes

The human HtrA family of serine proteases (HtrA1, HtrA2, HtrA3, and HtrA4) are the key enzymes associated with pregnancy and closely related to the development and progression of many pathological events. Previously, it was found that halogen substitution at the indole moiety of peptide Trp‐1 residue can form a geometrically satisfactory halogen bond with the Drosophila discs large, zona occludens‐1 (PDZ) domain of HtrA proteases. Here, we attempt to systematically investigate the effect of substitution with 4 halogen types and 2 indole positions on the binding affinity and specificity of peptide ligands to the 4 HtrA PDZ domains. The complex structures, interaction energies, halogen‐bonding strength, and binding affinity of domain‐peptide systems were modeled, analyzed, and measured via computational modeling and fluorescence‐based assay. It is revealed that there is a compromise between the local rearrangement of halogen bond involving different halogen atoms and the global optimization of domain‐peptide interaction; the substitution position is fundamentally important for peptide‐binding affinity, while the halogen type can effectively shift peptide selectivity between the 4 domains. The HtrA1‐PDZ and HtrA4‐PDZ as well as HtrA2‐PDZ and HtrA3‐PDZ respond similarly to different halogen substitutions of peptide; –Br substitution at R2‐position and –I substitution at R4‐position are most effective in improving peptide selectivity for HtrA1‐PDZ/HtrA4‐PDZ and HtrA2‐PDZ/HtrA3‐PDZ, respectively; –F substitution would not address substantial effect on peptide selectivity for all the 4 domains. Consequently, the binding affinities of a native peptide ligand DSRIWWV^–COOH as well as its 4 R2‐halogenated counterparts were determined as 1.9, 1.4, 0.5, 0.27, and 0.92 μM, which are basically consistent with computational analysis. This study would help to rationally design selective peptide inhibitors of HtrA family members by using different halogen substitutions.     

Stabilization of the integrase‐DNA complex by Mg2+ ions and prediction of key residues for binding HIV‐1 integrase inhibitors

The HIV‐1 integrase is an attractive target for the therapeutics development against AIDS, as no host homologue of this protein has been identified. The integrase strand transfer inhibitors (INSTIs), including raltegravir, specifically target the second catalytic step of the integration process by binding to the DDE motif of the catalytic site and coordinating Mg²⁺ ions. Recent X‐ray crystallographic structures of the integrase/DNA complex from prototype foamy virus allowed to investigate the role of the different partners (integrase, DNA, Mg²⁺ ions, raltegravir) in the complex stability using molecular dynamics (MD) simulations. The presence of Mg²⁺ ions is found to be essential for the stability, whereas the simultaneous presence of raltegravir and Mg²⁺ ions has a destabilizing influence. A homology model of HIV‐1 integrase was built on the basis of the X‐ray crystallographic information, and protein marker residues for the ligand binding were detected by clustering the docking poses of known HIV‐1 integrase inhibitors on the model. Interestingly, we had already identified some of these residues to be involved in HIV‐1 resistance mutations and in the stabilization of the catalytic site during the MD simulations. Classification of protein conformations along MD simulations, as well as of ligand docking poses, was performed by using an original learning method, based on self‐organizing maps. This allows us to perform a more in‐depth investigation of the free‐energy basins populated by the complex in MD simulations on the one hand, and a straightforward classification of ligands according to their binding residues on the other hand. Proteins 2014; 82:466–478. © 2013 Wiley Periodicals, Inc.     

Origins of PDZ domain ligand specificity. Structure determination and mutagenesis of the Erbin PDZ domain

The LAP (leucine-rich repeat and PDZ-containing) family of proteins play a role in maintaining epithelial and neuronal cell size, and mutation of these proteins can have oncogenic consequences. The LAP protein Erbin has been implicated previously in a number of cellular activities by virtue of its PDZ domain-dependent association with the C termini of both ERB-B2 and the p120-catenins. The present work describes the NMR structure of Erbin PDZ in complex with a high affinity peptide ligand and includes a comprehensive energetic analysis of both the ligand and PDZ domain side chains responsible for binding. C-terminal phage display has been used to identify preferred ligands, whereas binding affinity measurements provide precise details of the energetic importance of each ligand side chain to binding. Alanine and homolog scanning mutagenesis (in a combinatorial phage display format) identifies Erbin side chains that make energetically important contacts with the ligand. The structure of a phage-optimized peptide (Ac-TGW(-4)ETW(-1)V; IC(50) = approximately 0.15 microm) in complex with Erbin PDZ provides a structural context to understand the binding energetics. In particular, the very favorable interactions with Trp(-1) are not Erbin side chain-mediated (and therefore may be generally applicable to many PDZ domains), whereas the beta2-beta3 loop provides a binding site for the Trp(-4) side chain (specific to Erbin because it has an unusually long loop). These results contribute to a growing appreciation for the importance of at least five ligand C-terminal side chains in determining PDZ domain binding energy and highlight the mechanisms of ligand discrimination among the several hundred PDZ domains present in the human genome.     

Sequence-specific long range networks in PSD-95/discs large/ZO-1 (PDZ) domains tune their binding selectivity

Protein-protein interactions mediated by modular protein domains are critical for cell scaffolding, differentiation, signaling, and ultimately, evolution. Given the vast number of ligands competing for binding to a limited number of domain families, it is often puzzling how specificity can be achieved. Selectivity may be modulated by intradomain allostery, whereby a remote residue is energetically connected to the functional binding site via side chain or backbone interactions. Whereas several energetic pathways, which could mediate intradomain allostery, have been predicted in modular protein domains, there is a paucity of experimental data to validate their existence and roles. Here, we have identified such functional energetic networks in one of the most common protein-protein interaction modules, the PDZ domain. We used double mutant cycles involving site-directed mutagenesis of both the PDZ domain and the peptide ligand, in conjunction with kinetics to capture the fine energetic details of the networks involved in peptide recognition. We performed the analysis on two homologous PDZ-ligand complexes and found that the energetically coupled residues differ for these two complexes. This result demonstrates that amino acid sequence rather than topology dictates the allosteric pathways. Furthermore, our data support a mechanism whereby the whole domain and not only the binding pocket is optimized for a specific ligand. Such cross-talk between binding sites and remote residues may be used to fine tune target selectivity.     

Distinct binding specificity of the multiple PDZ domains of INADL, a human protein with homology to INAD from Drosophila melanogaster

PDZ domains are protein-protein interaction modules that typically bind to short peptide sequences at the carboxyl terminus of target proteins. Proteins containing multiple PDZ domains often bind to different trans-membrane and intracellular proteins, playing a central role as organizers of multimeric complexes. To characterize the rules underlying the binding specificity of different PDZ domains, we have assembled a novel repertoire of random peptides that are displayed at high density at the carboxyl terminus of the capsid D protein of bacteriophage lambda. We have exploited this combinatorial library to determine the peptide binding preference of the seven PDZ domains of human INADL, a multi-PDZ protein that is homologous to the INAD protein of Drosophila melanogaster. This approach has permitted the determination of the consensus ligand for each PDZ domain and the assignment to class I, class II, and to a new specificity class, class IV, characterized by the presence of an acidic residue at the carboxyl-terminal position. Homology modeling and site-directed mutagenesis experiments confirmed the involvement of specific residues at contact positions in determining the domain binding preference. However, these experiments failed to reveal simple rules that would permit the association of the chemical characteristics of any given residue in the peptide binding pocket to the preference for specific amino acid sequences in the ligand peptide. Rather, they suggested that to infer the binding preference of any PDZ domain, it is necessary to simultaneously take into account all contact positions by using computational procedures. For this purpose we extended the SPOT algorithm, originally developed for SH3 domains, to evaluate the probability that any peptide would bind to any given PDZ domain.     

Detailed regulatory mechanism of the interaction between ZO-1 PDZ2 and connexin43 revealed by MD simulations

The gap junction protein connexin43 (Cx43) binds to the second PDZ domain of Zonula occludens-1 (ZO-1) through its C-terminal tail, mediating the regulation of gap junction plaque size and dynamics. Biochemical study demonstrated that the very C-terminal 12 residues of Cx43 are necessary and sufficient for ZO-1 PDZ2 binding and phosphorylation at residues Ser (-9) and Ser (-10) of the peptide can disrupt the association. However, only a crystal structure of ZO-1 PDZ2 in complex with a shorter 9 aa peptide of connexin43 was solved experimentally. Here, the interactions between ZO-1 PDZ2 and the short, long and phosphorylated Cx43 peptides were studied using molecular dynamics (MD) simulations and free energy calculation. The short peptide bound to PDZ2 exhibits large structural variations, while the extension of three upstream residues stabilizes the peptide conformation and enhanced the interaction. Phosphorylation at Ser(-9) significantly weakens the binding and results in conformational flexibility of the peptide. Glu210 of ZO-1 PDZ2 was found to be a key regulatory point in Cx43 binding and phosphorylation induced dissociation.     

Lipid a from lipopolysaccharide recognition: Structure,dynamics and cooperativity by molecular dynamics simulations

Molecular dynamics simulations of Lipid A and its natural precursor Lipid IVA from E.coli have been carried out free in solution, bound to the myeliod differentiation protein 2 (MD2) and in the complex of MD2 with the toll like receptor 4 (TLR4). In addition, simulations of the ligand free MD2 and MD2‐TLR4 complex were performed. A structural and energetic characterization of the bound and unbound states of Lipid A/IVA was generated. As the crystal structures depict, the main driving force for MD2‐Lipid A/IVA are the hydrophobic interactions between the aliphatic tails and the MD2 cavity. The charged phosphate groups do strongly interact with positively charged residues, located at the surface of MD2. However, they are not essential for keeping the lipids in the cavity, indicating a more prominent role in binding recognition and ionic interactions with TLR4 at the MD2/TLR4 interface. Interestingly, in the absence of any ligand MD2 rapidly closes, blocking the binding cavity. The presence of TLR4, though changing the dynamics, was not able to impede the aforementioned closing event. We hypothesize that fluctuations of the H1 region are essential for this phenomenon, and it is plausible that an equilibrium between the open and closed states exists, although the lengths of our simulations are not sufficient to encompass the reversible process. The MD2/Lipid A‐TLR4 complex simulations show that the presence of the ligand energetically stabilizes the complex relative to the ligand‐free structures, indicating cooperativity in the binding process. © Proteins 2013. © 2012 Wiley Periodicals, Inc.     

Impact of calcium on N1 influenza neuraminidase dynamics and binding free energy

The highly pathogenic influenza strains H5N1 and H1N1 are currently treated with inhibitors of the viral surface protein neuraminidase (N1). Crystal structures of N1 indicate a conserved, high affinity calcium binding site located near the active site. The specific role of this calcium in the enzyme mechanism is unknown, though it has been shown to be important for enzymatic activity and thermostability. We report molecular dynamics (MD) simulations of calcium‐bound and calcium‐free N1 complexes with the inhibitor oseltamivir (marketed as the drug Tamiflu), independently using both the AMBER FF99SB and GROMOS96 force fields, to give structural insight into calcium stabilization of key framework residues. Y347, which demonstrates similar sampling patterns in the simulations of both force fields, is implicated as an important N1 residue that can “clamp” the ligand into a favorable binding pose. Free energy perturbation and thermodynamic integration calculations, using two different force fields, support the importance of Y347 and indicate a +3 to +5 kcal/mol change in the binding free energy of oseltamivir in the absence of calcium. With the important role of structure‐based drug design for neuraminidase inhibitors and the growing literature on emerging strains and subtypes, inclusion of this calcium for active site stability is particularly crucial for computational efforts such as homology modeling, virtual screening, and free energy methods. Proteins 2010. © 2010 Wiley‐Liss, Inc.