Comparative simulation studies of native and single-site mutant human beta-defensin-1 peptides

Human defensins play important roles in a broad range of biological functions, such as microbial defense and immunity. Yet, little is known about their molecular properties, i.e. secondary structure stability, structural variability, important side chain interactions, surface charge distribution, and resistance to thermal fluctuations, and how these properties are related to their functions. To assess these factors, we studied the native human β-defensin-1 monomer and dimer as well as several single-site mutants using molecular dynamics simulations. The results showed that disulfide bonds are important determinants in maintaining the defensins’ structural integrity, as no structural transitions were observed at 300 K and only minor structural unfolding was detected upon heating to 500 K. The α-helix was less thermally stable than the core β-sheet structure held together by hydrogen bonds and hydrophobic interactions. The monomer α-helix stability was directly correlated, whereas the end-to-end distance was inversely correlated to the experimentally measured β-defensin-1 chemotactic activity, in the order: mutant 2 (Gln24Glu) > mutant 3 (Lys31Ala) = wild type > mutant 1 (Asn4Ala). The structural stability of the β-defensin-1 dimer species exhibited an inverse correlation to their chemotactic activity. In dimers formed by mutants 2 and 3, we observed sliding of one monomer upon the surface of the other in the absence of unbinding. This dynamic sliding feature may enhance the molecular oligomerization of β-defensin-1 peptides contributing to their antibacterial activity. It could also help these peptides orient correctly in the CC chemokine receptor 6 binding site, thereby initiating their chemotactic activity. In agreement with this notion, the remarkable sliding behavior was observed only for the mutants with the highest chemotactic activity.

The basic β-defensin structure is composed of an N-terminal α-helix and three antiparallel β-strands, the latter held together by three intramolecular disulfide bonds. The most important conserved motif is the GXC in the second β-strand, where G is Gly, C is the fourth Cys in the peptide, and X is any amino acid (Hoover, Chertov, & Lubkowski, 2001). In the case of HBD-1, this motif is represented by Gly25-Thr26-Cys27. The proline residue in the β1-β2 turn of β-defensins, most of the Gly residues (Hoover et al., 2000), and the β-bulge at the conserved Gly found in the β2-strand, are also conserved (Hoover et al., 2001). Based upon X-ray studies, the reported HBD-1 structure is comprised of 36 amino acids having six cysteines forming three conserved disulfide bonds, namely Cys5/Cys34, Cys12/ Cys27, and Cys17/Cys35 (Hoover et al., 2001). HBD-1 is mostly reported to be a monomer, but a structure of four HBD-1 monomers in one unit cell, not biologically significant, has also been resolved (Hoover et al., 2001). The X-ray structures of single-site mutants of HBD-1 have been reported as well (Pazgier et al., 2007). NMR diffusion measurements of the radius of hydration (Schibli et al., 2002) and Nuclear Overhauser effect experiments showed that HBD-1 exists in the monomeric form in aqueous solution (Bauer et al., 2008).
Despite the established physiological importance of the human β-defensins, only a few computer-based modeling studies for human β-defensins have been reported to date (Fung, Floudas, Taylor, Zhang, & Morikis, 2008;Sharadadevi & Nagaraj, 2010;Suresh & Verma, 2006). In one paper (Sharadadevi & Nagaraj, 2010), a molecular dynamics (MD) simulation study using GROMACS package (Berendsen, van der Spoel, & van Drunen, 1995;Lindahl, Hess, & van der Spoel, 2001;van der Spoel et al., 2005) has been carried out to study HBD-1 conformational behavior with and without disulfide bridges. The authors reported the unfolding behavior of HBD-1 and the distances between sulfur atoms during the course of the MD simulations (Sharadadevi & Nagaraj, 2010). Although HBDs are good candidates for computational studies due to their small size, the literature is lacking MD simulation studies for HBDs, and there are no simulation studies addressing HBDs chemotactic properties. Only one study (Pazgier et al., 2007) reported single-site mutants and compared their biological activities to the wild type (WT). Here, we carried out comprehensive simulation studies of the structural and dynamic properties of the native HBD-1 peptide in its monomeric and dimeric forms as well as three single-site mutants of HBD-1 using all-atom MD simulations. We have chosen these specific mutants because they exhibit important functional differences. They vary in their level of chemotactic and antibacterial activities compared to the WT. The monomer and dimer structures have been studied because they are significant biologically.
Here, we report the main results of our studies, in which we resolved and contrasted the structural and dynamic properties of the native HBD-1 peptide with those of the single-site mutants. For each system, we analyzed the secondary structure dynamics using secondary structure propensity and Ramachandran plots. We examined dynamics of binary contacts between side chains utilizing contact maps and maps of stable contacts, molecular flexibility using root-mean-square deviations (RMSD) and root-mean-square fluctuations (RMSF), surface charge distributions employing solvent accessible surface areas (SASA), and the end-to-end distances. In this comparison, we related the observed structural and dynamic changes to the altered side chains of the mutants, and discussed their implications for the experimentally measured functional differences in chemotactic activity reported for these mutants (Pazgier et al., 2007). The most significant findings are the following. We found correlations between the measured chemotactic activity and some of the molecular properties, such as the α-helix stability, end-to-end distance, and stability of the dimeric form (Table 1). The dimers formed by mutants Mut2 (Gln24Glu) and Mut3 (Lys31Ala) exhibited sliding and rolling behavior of the two monomer chains over each other. The basic residues were found to be exposed to solvent, which points to their importance in antibacterial and chemotactic activity of HBD-1.

Structural models
The molecular models have been constructed using the X-ray structures from the Protein Data Bank (PDB) (Berman et al., 2000) (see Table S1). The initial structures for the HBD-1 mutants in their dimeric form were obtained by substituting the appropriate residue in the HBD-1 WT dimer structure using Visualization Molecular Dynamics (VMD) mutator (Mutator Plugin, Version 1.3: http://www.ks.uiuc.edu/Research/vmd/plugins/mutator/). These substitutions are Asn4Ala (Mut1), Gln24Glu (Mut2), and Lys31Ala (Mut3). The rationale for choosing these three mutants was that their chemotactic activity was significantly different when compared to the WT; Mut1 (Mut2) showed lower (higher) chemotactic activity, whereas Mut3 was almost as active as the WT (Pazgier et al., 2007) (see Table 1). The mutation sites were also diverse; point mutations were in the α-helix (Mut1), in the β2-strand (Mut2), and in the β2-β3 turn (Mut3) (Pazgier et al., 2007).

Computer simulations
Equilibrium MD simulations: All-atom MD simulations in explicit water were performed using the NAMD package (Phillips et al., 2005) with CHARMM22 force field (MacKerell et al., 1998) and TIP3 model for water molecules (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983). For each peptide, a water box has been created so that in each direction x, y, and z the distance between the molecule and the box boundaries was P12 Å. For example, the box size for the WT monomer was 48.5 Å Â 50.4 Å Â 45.8 Å. All the systems have been neutralized using chloride counterions. The density of each system was maintained at 1 g/cm 3 to mimic ambient water conditions. After initial minimization, all the systems were heated to T = 300 K at constant volume and equilibrated at T = 300 K for 4 ns. Long 30 ns equilibrium simulation runs have been performed for each system, using the Number of Particles, Pressure, Temperature (NPT) ensemble and a 1 fs integration step. The van der Waals interactions were gradually turned off at a distance between 10 and 12 Å. The list of nonbonded interactions was truncated at 13.5 Å. The short-range nonbonded interactions and the long-range electrostatic forces were computed for every step. The particle mesh Ewald (PME) summation was employed to describe the long-range electrostatics (Darden, York, & Pedersen, 1993;Toukmaji & Board, 1996). Heating simulations: We carried out 10 ns heating simulations, in which the obtained structures equilibrated at T = 300 K have been heated from T in = 300 K to T fin = 500 K. GROMACS simulations: Additional 25 ns simulation runs for the native WT monomer (Table S1) have been performed using the GROMACS package (Berendsen et al., 1995;Lindahl et al., 2001;van der Spoel et al., 2005) (GROMOS96 45a3 (Schuler, Daura, & van Gunsteren, 2001) and Optimized Potentials for Liquid Simulations (OPLS) (Jorgensen, Maxwell, & Tirado-Rives, 1996) force fields to compare our findings with the results reported in Sharadadevi and Nagaraj (2010) and to exclude possible correlations between the results obtained and the force fields used. In these runs, we employed the same simulation setup used in Sharadadevi and Nagaraj (2010), albeit with a larger water box (49.4 Å Â 51.3 Å Â 46.7 Å). In simulations using GRO-MOS96 and OPLS force fields, we also utilized PME method (Darden et al., 1993;Toukmaji & Board, 1996). A summarized description of all MD simulation runs is given in Table S1.

Data analyses
We analyzed the secondary structure dynamics using time evolution of secondary structure propensity, i.e. propensities (probabilities) of amino acid residues to form different secondary structure motifs (e.g. α-helices, β-strands, etc.), and the Ramachandran plots (distribution of dihedral angles Φ and Ψ). We used the structure overlap function (χ) (Camacho & Thirumalai, 1993), 0 6 χ(t) 6 1, to estimate the extent of similarity of the transient conformation Table 1. A summary comparison of experimental chemotactic activity with the simulation results for the native HBD-1 peptide (WT) and the three single-site mutants (Mut1, Mut2 and Mut3).
Experimental chemotactic activity (from Klüver et al., 2005) III (100 ± 11%) IV (40 ± 8%) I (130 ± 8%) II (109 ± 14%) Stability of α-helix in monomeric form ( For each attribute, the four molecules are ranked in descending order using Roman numerals with I being the highest and IV being the lowest. A slight ranking difference between our findings and the reported data is shown in italics. observed at time t to the reference structure (initial state) at time t = 0. When χ = 0, the two structures are completely dissimilar, and when χ = 1 they are exactly the same. The structure overlap function is given by where, N is the number of atoms in the system, h is the Heaviside step-function (h(x) = 1 when x > 0 and zero otherwise), and δ = 0.2 is the tolerance to the thermal fluctuations (at T = 300 K). Here, r ij (t) is the distance between the ith and jth atoms at time t, and r 0 ij is the same distance in the reference state (t = 0). The contact maps of interacting amino acids represented the trajectory averaged picture of binary contacts between side chains (both inter and/or intramonomer interactions). A contact between centers of mass (COM) of side chains of a pair of residues, for which the COM-COM distance is <6.5 Å cut-off, is counted as a contact. Dynamic maps of stable contacts is a powerful tool to follow the pattern of formation and rupture of stable binary contacts between coupled residues (both native and non-native) in a dimer. This tool has been used to study protein-protein interactions in the HBD-1 dimers. All the stable contacts, which persisted for at least 1000 ps in a simulation run, have been displayed as a function of time. The contact maps and plots of the total number of intermonomer contacts have been used to analyze the transient structures and kinetics of protein-protein interactions. The RMSD, which quantify the average displacements of C α -particles relative to their positions in the reference state, were used to probe the global peptide fluctuations. The RMSF allow one to access local conformational fluctuations. We used a "running average" structure as a reference state in each 100 ps time interval to eliminate rotational and translational contributions (Kabsch, 1976(Kabsch, , 1978. The SASA per residue provide insights into solvent accessibility of the side chains on the surface of the molecule (Ferrara, Apostolakis, & Caflisch, 2002). This measure has been used to describe potential surface interactions and receptor binding properties. The plots of the total SASA vs. time have been used to measure these properties of the peptide molecule as a whole. We monitored the end-to-end distance as a function of time, as a suitable measure of the global conformational transitions in the HBD-1 peptide and its derivatives, such as unfolding, thermal denaturation, and dissociation (in HBD-1 dimers).

Structural features
The compact structure of the 36 amino acid native HBD-1 WT (Figure 1  In Panels (a)-(c), the WT monomer is shown in cyan, and chains A and B are shown in blue and red, respectively; the cysteine residues and disulfide bonds (Cys5/Cys34, Cys12/Cys27, and Cys17/Cys35) are represented as yellow lines and bonds, respectively. In Panel (d), the hydrophobic, basic, acidic, and polar residues are shown, respectively, in white, blue, red, and green color.
(b)) was stabilized by the β-sheet and three conserved disulfide bonds between Cys residues 5 and 34, 12 and 27, and 17 and 35 (Hoover et al., 2001). In the native WT peptide, residues His2-Ser7 formed an α-helix, while residues Gln11-Leu13, Ile23-Cys27, and Ala32-Cys35 formed three β-strands connected by two turns localized to residues Pro18-Thr21 and Tyr28-Lys31 (Hoover et al., 2001). The native HBD-1 dimer structure, characterized by an asymmetric arrangement of chains A and B (Figure 1(b)), was stabilized by the intermonomer interactions among amino acids in the β2-strand of chain A and in the β2-β3 turn of chain B. All the conserved residues and structure elements (Hoover et al., 2000(Hoover et al., , 2001 discussed earlier were retained in all defensins we simulated in this study. Comparing the WT monomer equilibrated for 4 ns (Figure 1(c)) to its initial structure (Figure 1(a)), we found that the α-helix and β-sheet showed some structural differences in terms of their orientation with respect to the bulk of the molecule (Figure 1). Two surface representations of different orientations of the equilibrated WT monomer (Figure 1(d)) showed extensive patches of residues with similar charge and polarity. Specifically, we observed basic and hydrophobic patches that are known to be crucial for binding and inserting into bacterial membranes (Lehrer et al., 1989;Morgera et al., 2008;Wimley et al., 1994). These findings correlate well with the known antibacterial activity of defensins (García, Jaumann, et al., 2001;Klüver et al., 2005;Krishnakumari et al., 2003;Mandal et al., 2002), and possibly their role in chemotaxis (Pazgier et al., 2007). The results of heating simulations for the native peptide WT and mutants Mut1, Mut2, and Mut3 showed that, for all the systems, the number of hydrogen bonds decreased and the radius of gyration increased with increasing temperature from T in = 300 K to T fin = 500 K (data not shown). This should be expected as proteins tend to expand with increasing temperature while losing some or all of their hydrogen bonds that stabilize their native fold. However, this effect was minimal due to the constraints posed by the three disulfide bonds.

Secondary structure dynamics
We examined the stability of the secondary structure elements by analyzing the time evolution of secondary structure propensity of amino acid residues and the Ramachandran plots. The β-sheet structure formed from the three antiparallel β-strands showed essentially similar Figure 2. Time evolution of the secondary structure propensity for the α-helical portion (residues 2-7) of the monomers of native HBD-1 peptide (WT) and mutated HBD-1 peptides (Mut1, Mut2, and Mut3) from equilibrium simulations at constant temperature T = 300 K (left Panel), and from thermal unfolding simulations at increasing temperature (from T in = 300 K to T fin = 500 K; right Panel). The α-, π-, and 3 10 -helices are shown in purple, red, and blue color, respectively; unstructured turn and coil regions are presented, respectively, in cyan and white color. The unstructured regions, observed at the end of the simulations of thermal unfolding for WT and Mut3 are encircled. behavior in the WT, Mut1, Mut2, and Mut3. This is hardly surprising since the three conserved disulfide bond topology exists in all the peptides, giving rise to a highly stable β-sheet structure. The time evolution of secondary structure propensity for the α-helical region in WT, Mut1, Mut2, and Mut3 showed some minor differences ( Figure 2). No pronounced distortion of the secondary structure was observed in equilibrium MD simulation runs, although His2 and/or Tyr3 in WT and α-helix in Mut1 were found to be partially unstructured ( Figure 2). All simulation runs for all the peptides showed the presence of some π-helix character except for Mut2, which had the most stable α-helix. Mut1 had the least stable α-helix (Figure 2), and the order of decreasing stability for the α-helix was Mut2 > Mut3 ≈ WT > Mut1. This correlated well with the experimentally measured chemotactic activity of HBD-1: Mut2 (130 ± 8%) > Mut3 (109 ± 14%) = WT (100 ± 11%) > Mut1 (40 ± 8%) (Pazgier et al., 2007) (Table 1).
We compared the structural integrity of the native HBD-1 peptide and its mutants in the monomeric and dimeric forms ( Figure 3). Mut1 (Mut2) dimer was the most (least) stable in terms of secondary structure, thus, giving an inverse correlation with the reported HBD-1 chemotactic activity (Pazgier et al., 2007) (Table 1). This implies that HBD-1 monomers are more active CCR6 ligands than their dimeric counterparts, similar to the MIP-3α monomer that is reported to be more adept (than the dimer) at binding to the CCR6 receptor (Chan, Hunter, Tack, & Vogel, 2007). The α-helix in the Mut2 monomer was more stable than the same structural elements in both chains comprising the Mut2 dimer, showing less π-helical character ( Figure 3). The structurally asymmetric dimer chains A and B exhibited different dynamic behavior. For example, the dimer β1-strand of chain A was the least stable secondary structure element, while β3-strand of chain A was the most stable.
We carried out heating simulations for WT and the three mutants to probe the thermal stability of these peptides, and to test the strength of the disulfide bonds. As shown in Figure 2, at T = 460 K, the α-helix in the WT molecule started to unravel from Asp1 through Cys5 (i.e. up to the first disulfide bond position), and remained unstructured up to T = 500 K. Mut1 was relatively stable during heating; at T = 500 K residues 2-4 formed a 3 10 -helix ( Figure 2). Hence, introducing an Asn4Ala substitution (Mut1) in WT thermally stabilized the α-helix. Upon heating Mut1, transient β1-strand deformation was observed in the 430-450 K temperature range (data not shown). During heating of Mut2, no significant changes to the α-helix were observed ( Figure 2); it was only transiently deformed. In heating simulations for Mut3, at T = 475 K, the α-helix formed a turn and at T = 460 K the β1-strand completely lost structure ( Figure 2).
The results of analysis of secondary structure propensities and Ramachandran plots for all HBD-1 monomers and dimers showed that no significant structural changes occurred. At higher temperature, residues 2-7 in WT and Mut3 showed slightly lower α-helix propensities and some loss of structure (see also Figure 2). The Ramachandran plots for Mut2 monomer, obtained from equilibrium and heating simulations, exhibited similar distributions of amino acids in the β-sheet region but different distribution for the α-helical region (Figure 4(a)). As expected for Mut2 monomer, residues 3-7 showed lower α-helical character at elevated temperature (Figure 4(c)) compared to room temperature (Figure 4(b)).

Tertiary structure analysis
For WT and the three mutants, we compared the final structures obtained at the end of 30 ns equilibrium simulation runs with the corresponding initial structures using the structure overlap function χ (see "Methods" section). We found that the former and the latter were not very dissimilar. We ranked HBD-1 peptides in decreasing order of their χ value, averaged over three independent trajectories, as follows: Mut3 (0.96) > Mut2 (0.94) > WT (0.92) > Mut1 (0.9). Although these differences were not statistically significant, their order agreed with the ranked chemotactic activities measured experimentally (Pazgier et al., 2007), Mut2 > Mut3 = WT > Mut1 (Table 1), indicating that a direct correlation may exist between the extent of conformational fluctuations in HBD-1 monomers and their chemotactic activity. We also carried out a similar analysis for the dimers of the WT and the mutants and observed a greater variation in χ as compared to their monomeric forms. The highest χ was observed for the WT dimer (0.96) and Mut1 dimer (0.95), and the lowest χ was observed for the Mut2 dimer (0.83). These results were opposite to those of the monomers and agreed with the differences in secondary structure stability discussed above, where Mut1 (Mut2) dimer was the most (least) stable (Table 1).
We next analyzed transient and final structures obtained in the course of heating simulations. The WT showed a large decrease in χ ( Figure 5), which corresponds to partial unraveling of the α-helix starting at T = 460 K (Figure 2). The curve of χ for Mut1 showed a transient drop at T = 430-450 K due to deformation of the β1-strand mentioned above. For Mut2 and Mut3, χ varied considerably and after reaching T = 460 K decreased sharply to χ ≈ 0.8 and 0.7, respectively (Figure 5). For Mut3, the decrease in χ at T = 460 K was due to the β1-strand becoming more unstructured, which was concomitant with unraveling of the α-helix discussed previously (see also Figure 2). We used the final χ values observed in 10 ns heating simulations ( Figure 5) to rank structural stability of the HBD-1 monomers against the thermal perturbation. We obtained the following order: Figure 3. Time evolution of the secondary structure propensity for the mutated HBD-1 peptide Mut2 monomer (Panel (a)) and dimer (Panel (b)) obtained from equilibrium simulations at room temperature T = 300 K. The β-sheets are shown in yellow. The color denotation for α-, π-, and 3 10 -helices, and for unstructured turn and coil regions is the same as in order. These quantities were positively correlated with the experimental chemotactic activities reported for these molecules (Pazgier et al., 2007) (Table 1).

Dynamics of binary contacts
An examination of the dynamics of formation and rupture of intra and intermonomer contacts allowed us to gather residue level information to help us explore a temporal evolution of protein-protein interactions in the HBD-1 peptides. Indeed, contact maps, extracted from equilibrium and heating simulations, clearly illustrated the specific interacting residue sites and their stability at equilibrium (T = 300 K) and at elevated temperature (from T in = 300 K to T fin = 500 K). We analyzed the dynamics of intramonomer contacts for monomers of the WT and the mutants, and the dynamics of intermonomer contacts for dimers.
For the Mut1 (Asn4Ala) monomer, Ala4 formed stronger and longer lived contacts than Asn4 of the WT monomer. The polar contact of Asn4 with residue Gln24 was present in the WT monomer but was absent in Mut2 (Gln24Glu). The contacts formed in Mut3 (Lys31Ala) by Ala31 with both Tyr14 and Ser15 were weaker than those the Ramachandran plot (i.e. distribution of dihedral angles Φ and Ψ) for the mutated HBD-1 peptide Mut2 obtained from equilibrium simulations at room temperature (T = 300 K; black dots), and from heating simulation at increasing temperature (from T in = 300 K to T fin = 500 K; green squares). Each data point represents one residue whose dihedral angles were averaged over the simulation run. Panel (b): The propensity of residues 2-7 in the mutated HBD-1 peptide Mut2 to form α-helix obtained from equilibrium simulations. Panel (c): same as in Panel (b) but obtained from heating simulation at increasing temperature from T in = 300 K to T fin = 500 K. Figure 5. The structure overlap function χ vs. time for the native HBD-1 peptide (WT) and mutated peptides (Mut1, Mut2, and Mut3) obtained from heating simulations at increasing temperature from T in = 300 K to T fin = 500 K. Color and plot denotations are presented in the graph. formed by Lys31 in the WT monomer. The weaker interactions were due to the shorter methyl group side chain of Ala compared to the three methylene bridge in the Lys side chain. All the contacts formed by the basic residues were monitored due to their importance in the antibacterial and chemotactic activities (García, Jaumann, et al., 2001;Klüver et al., 2005;Krishnakumari et al., 2003;Mandal et al., 2002;Pazgier et al., 2007). In all the monomer equilibrium simulations, we observed Lys33 and Lys36 establishing binary contacts, while Lys22, Arg29, and Lys31 formed no contacts. No significant differences were detected in the distributions of contacts obtained from equilibrium simulations and heating simulations. However, in heating simulations a few minor changes were detected in the α-helix in WT and Mut3 (Figure 2; see also Figure 6), and in the β2-strand of Mut2 (Figure 3(a)). In Mut3, the contacts between residues in the α-helix (residues 2-7) and β2-strand (residues 23-27) became progressively weaker as shown in the contact maps and structural snapshots for T = 320, 400, and 480 K (Figure 6). At 320 K, only Tyr3 and Ser7 in the α-helix were lacking contacts; the other four residues all formed strong contacts. At T = 400 K, the α-helix moved away from the β2-strand ( Figure 6; see also Figure 2); Asn4, Cys5, and Val6 were still forming strong contacts. At 480 K, the α-helix unfolded ( Figure 2) and moved away from the β2-strand resulting in significantly weaker contacts ( Figure 6). As expected, the exception to this trend was observed for Cys5 (disulfide bond).
Equilibrium simulations for the Mut1 dimer revealed stable contacts in the α-helical region, especially at Asn4 position (data not shown). These contacts were not detected in the WT dimer, which could be due to the Asn4Ala mutation creating new contacts at position 4. This might be one of the reasons why the Mut1 dimer was more stable than the WT dimer. The curve of χ for the Mut2 dimer showed a decrease in χ values at 30 ns ($0.60); the decreases in χ values at $6, 18, 22, and 24 ns were accompanied by decreases in the number of stable contacts (Figure 7(a) and (b)). At the end of the simulation run, the Mut2 dimer had different coupled residues compared to the initial structure (Figure 7(b)), strongly suggesting that the two chains had been sliding over Figure 6. Side chain contact maps for mutated HBD-1 peptide Mut3 obtained from heating simulation at T = 320, 400, and 480 K. The strength of intramolecular binary contacts between interacting amino acids varied from the weakest (shown in white color) to the strongest (black); the yellow squares with the letter 'C' point out the Cys5/Cys34 disulfide bonds in the α-helix. Contacts formed between residues Asp1-Val6 in the α-helix and residues Gln24-Tyr28 in the β2-strand are circled. In the rightmost Panel, the residues forming the secondary structure elements are boxed; α-helix: residues 2-7, and β-sheet: residues 11-13, 23-27, and 32-35. The amino acid residues mediating the α-helix interaction with the β2-strand are shown explicitly. each other. For the Mut2 dimer, we observed formation of the stable contacts Glu24A/Thr26B, Glu24A/ Arg29B, Glu24A/Tyr3B, and Glu24A/Asn4B, which were not detected for the WT dimer. This was due to the different physicochemical properties of the Gln and Glu residues. The Mut3 mutation (Lys31Ala) influenced stable contacts in the Mut3 dimer, where residue Ala31 in chain B formed weaker stable contacts than residue Lys31 in chain B in the other HBD-1 dimers. In all runs for the Mut3 dimer, we observed formation of contacts Thr26A/Ala31B, Ile23A/Ala31B, and Thr26A/Ala31B (Figure 7(d)). For Mut3, some transition regions in the time-dependent map of stable contacts are obvious. However, we see that some contacts might suddenly appear or disappear, or might exist only for a short period of time (Figure 7(d)). The two successive snapshots (Figure 7(c)) correspond to the structures observed at 15 and 18 ns on either side of the abrupt 17 ns drop in χ to $0.65. In these two representations, the two monomer chains are clearly in different relative orientations. This conformational transition is reflected in the changes in the map of stable contacts (Figure 7(d)), which can be interpreted in terms of the monomers sliding to their new positions without any attempt to dissociate.
native state at the end of a 30 ns run (Figure 8). The WT dimer was stable maintaining $15-20 contacts in the course of a 30 ns run. The Mut2 dimer was the least stable in terms of the total number of intermonomer contacts. We also examined the Q vs. t pattern for the WT and all mutant dimers, and arranged them in the following order of their decreasing stability: Mut1 > Mut3 > WT > Mut2. We found that the higher dimer stability corresponded to the weaker measured chemotactic activity (Pazgier et al., 2007) (with a slight difference in WT and Mut3 ranking).

Root-mean-square fluctuations and deviations
The RMSF for each ith amino acid residue, where x i (t j ) andx i are the positions for C αatom i at time t j and average value of the reference position of C α atom i, respectively, and t tot is the length of the simulation run (in ns), provides a measure of molecular fluctuations at the local scale. We calculated RMSF using the results of equilibrium simulations at T = 300 K and the PDB structures for the monomers and dimers of the WT and the three mutants. For all monomeric systems, the RMSF profiles were similar and the RMSF values varied in the 1.8-2.3 Å range (data not shown). The highest RMSF values were observed for the residues forming the turns Ser8-Gly10, Tyr14-Lys22, and Tyr28-Lys31, the N-terminus (Tyr3) and the C-terminus (Lys36). By comparing the results of equilibrium and heating simulations, we found that at elevated temperatures the WT monomer had lower RMSF values ( Figure S1(a)). This may be due to the hydrophobic effect, where solvent exclusion provides increasing stability with rising temperatures. This effect was observed for all the HBD-1 peptides except for Mut3 ( Figure S1(b)). The RMSF analysis for dimers revealed that the WT dimer was the most stable compared to Mut1, Mut2, or Mut3 (smallest RMSF; data not shown). Only the first 10 residues of Mut1 dimer showed lower RMSF values than those of WT, including the mutation site Asn4Ala. This was due to the shorter side chain of Ala vs. the Asn replacement, thereby showing a smaller fluctuation at that position. In the Mut2 dimer, the Gln24Glu mutation did not result in a larger RMSF value at that position compared to the WT dimer (data not shown), perhaps because both Glu and Gln had similar side chain size. The two chains in a dimer interacted thereby affecting RMSF for each chain, unlike two stand-alone monomers. Accordingly, the RMSF profile for the dimer in question was not identical to the RMSF profile for the monomers, and the two dimer chains were also different especially at the dimerization sitethe region of functional asymmetrywhere fluctuations are expected to be minimal. This together with the time evolution of secondary structure propensity Figure 8. The total number of intermonomer contacts between interacting amino acid residues in chains A and B (see Figure 1) obtained from equilibrium simulations at T = 300 K for the dimers of native HBD-1 peptide (WT; average of two trajectories) and mutated peptide (Mut2). The three representative simulation runs, displayed by solid black, dotted red, and solid green curves, show stochastic variability in the number of intermonomer association interactions. The first trajectory for the Mut2 dimer (solid black curve), also discussed in Figure 6, had noticeable structural changes that occurred at $6, 18, 22, and 24 ns (marked by vertical arrows). In the second run for the Mut2 dimer (dotted red curve), a complete dissociation transition was observed at t ≈ 22 ns, which was followed by redimerization at t ≈ 25 ns. (Figure 3(b)), corroborated the theory of dimer asymmetry. In the equilibrium simulations of the WT dimer, the dimerization site was identified to be the β2-strand of chain A interacting with the residues in β2and β3strands, and β2-β3 turn in chain B (residues 23-27; see Figure 9). This site had low fluctuations (Figure 9(b)) and, consequently, higher interaction levels on the side chain contact map (Figure 9(a)). Aside from the dimerization site, all dimers showed the minimum RMSF value at all Cys residues due to the disulfide bond con- Figure 9. Panel (a): the time-averaged map of interacting residues (side chain contact map) in chain A (shown in blue color) and chain B (red) for the dimer of WT HBD-1 peptide obtained from equilibrium simulation at T = 300 K. The side chain contact map shows structural evidence for the asymmetric dimerization site between chains A and B, encircled in the structure and on the map. The contacts' representations ranged from black (strongest coupling) to white (weakest interaction) depending on the contact strength throughout the simulation run; the yellow squares point out the position of the disulfide bonds (Cys 5/34, 12/27, and 17/35). Panel (b): the RMSF for amino acid residues in chains A and B for the same run displayed in Panel (a), showing similar yet nonidentical pattern for the two dimer chains. Higher RMSF values were observed at the C-and N-termini and in the turn regions, while the minimum RMSF values were observed for the α-helical portion, disulfide bonds, and the dimerization site localized to residues 23-26 and 26-33 of chain A and B, respectively (encircled in red). straints, and at the β-strands due to the formation of the β-sheet architecture.
We also analyzed the time dependence of the RMSD given by as a molecular scale measure of dynamic transitions in the HBD-1 peptides, where x i (t j ) is the position of the ith C α atom at time t j in the average structure; x i (t 0 ) is the reference position for the same atom in the initial structure (t = 0), and N is the total number of C α atoms. In Figure 10, we display representative trajectories of RMSD for the dimers of WT, Mut2, and Mut3 (one of the three simulation runs for Mut2 and Mut3 are same as in Figure 7). The RMSD profile for the WT dimer remained flat for the entire 30 ns simulation run at T = 300 K, indicating that the native WT dimer was stable ( Figure 10(a)). We see that the onset of the conformational transition in Mut2 and Mut3 coincided with sudden changes in their RMSD values (Figure 7(b) and (c)). For example, there was a noticeable increase in the RMSD value for Mut3 dimer after 16 ns (Figure 10(c)), which negatively correlated with the sharp drop in the χ value at $17 ns (Figure 7(c)), due to the redistribution of stable bonding contacts in the dimer (Figure 7(d)). Here, the interaction pattern, in which the β2-strand of chain A formed contacts with the β2-strand and β2-β3 turn of chain B, changed to the new pattern, in which the β2-β3 turns in both chains interacted together strongly through Figure 10. The RMSD for the dimers of native (WT) HBD-1 peptide (Panel (a)), and mutants Mut2 (Panel (b)) and Mut3 (Panel (c)) obtained from equilibrium simulations at T = 300 K. The first, second, and third simulation runs are represented by the solid black, dotted red, and solid green curves, respectively. The solid black (dotted red) trajectory for Mut2 (Mut3) dimer are discussed in Figure 6. The WT dimer was found to be the most stable (RMSD ≈ 2 Å), while the curves for Mut2 and Mut3 showed the highest RMSD (≈ 7-10 Å). In Panels (b) and (c), the vertical arrows mark the times at which noticeable structural changes had occurred as discussed in Figure 6.
the Gly30A/Arg29B and Arg29A/Tyr28B contacts, and β2-strand coupled to β2-β3 turns through the Thr26A/ Arg29B contact (see also Figure 7(c) for relevant structures). Hence, chain A appeared to be sliding over chain B without any dissociation attempt, as discussed previously. Similarly, the RMSD value for the Mut2 dimer showed stepwise increases (Figure 10(b)), which correlated well with the previously described sudden drops in χ value at $6, 18, 22 and 24 ns (Figure 7(a)), and weakening of the dimerization site (see Figure 7(b)).

Solvent accessible surface area
We analyzed the SASA per residue averaged over three equilibrium simulation runs for the monomer WT and mutants. This measure provides information about the accessibility of the side chains of hydrophilic amino acids exposed to water, and, hence, it can be used to single out and monitor the residues participating in associations with other biomolecules. Surprisingly, we observed roughly similar SASA values for all residues and for all monomers during equilibrium simulations at T = 300 K except for the mutation sites (data not shown). Specifically, substitutions Asn4Ala (Mut1) and Lys31Ala (Mut3) showed lower while Gln24Glu substitution (Mut2) showed slightly higher SASA values compared to the WT. We attributed this to the length of the side chain; the shorter the side chain, the less it would be exposed to the solvent. The results for SASA calculations for heating simulations (Figure 11) revealed the same qualitative pattern. All the mutation sites were less solvent exposed than their WT counterparts. The relative similarity between the results of equilibrium and heating simulations imply that, although at high temperature the HBD-1 peptides loose some of their secondary structure, the overall tertiary structure, which determines the surface arrangement of polar residues, is preserved. Equilibrium simulations at T = 300 K for dimers revealed that the SASA patterns for each monomer chain in dimeric and monomeric forms were very similar (data not shown). This stands in contrast to the results obtained from RMSF analysis (Figure 9(b)). Some residues, e.g. Ala16, Ile19, Lys33 and Lys36 in the case of WT dimer, showed noticeable differences between the two monomeric chains in the dimer, which also corroborated our previous findings for the secondary structure propensity (Figure 3(b)) and RMSF (Figure 9(b)) regarding asymmetrical dimers. We also analyzed the total (molecular) SASA as a function of time. A representative example, obtained from heating simulations for the WT monomer ( Figure S2(a)) showed higher exposure to water at the end of the heating simulations (≈4200 Å 2 ) compared to the equilibrium simulations (≈4000 Å 2 ). This coincided with the α-helix unraveling in the WT monomer at high temperature (Figure 2), making the Nterminal more water exposed. A similar SASA pattern was observed for Mut3 monomer (data not shown). For the Mut1 monomer, equilibrium and heating simulations showed the same SASA values because, as we described previously, the molecule remained the most stable. During the heating simulation for the Mut2 monomer (Figure S2(b)), the SASA value increased at T = 470 K due to the β2-strand deformations, and at T = 475 K due to the transiently deformed α-helix (Figure 2).

End-to-end distance
The global transitions in proteins, such as unfolding, can be accessed by analyzing their end-to-end distance. We examined this quantity for the WT and the mutants monomers using the results of equilibrium and heating Figure 11. The SASA per residue obtained from heating simulations of the native HBD-1 peptide (WT, shown in solid black color) and mutants Mut1 (dotted red), Mut2 (solid green), and Mut3 (dashed blue) at increasing temperature (from T in = 300 K to T fin = 500 K). Each mutant had smaller SASA value at its mutation site (encircled) as compared to the WT peptide with the smallest SASA difference for Mut2 (substitution Gln24Glu).
simulations. At higher temperature, WT and Mut3 showed larger end-to-end distance values than their average end-to-end distance obtained from three equilibrium simulation runs at T = 300 K, mostly due to α-helix unfolding at T = 460 K for WT and T = 475 K for Mut3 displayed in Figure 2 (see results for WT in Figure S2 (c)). By contrast, the end-to-end distance for Mut1 remained constant at $14 Å all through the course of heating simulations (data not shown), which was due to the relative stability already noted for that monomer (Figure 2). Upon heating, the end-to-end distance for Mut2 initially increased gradually to $17 Å ( Figure S2 (d)), following the changes in the secondary structure, i.e. partial deformations of the β2-strand and α-helix at around T = 476 K, and then dropped to $12 Å due to restoration of the α-helix ( Figure S2(d); see, also, Figure 2). We ranked the monomer WT and mutants in the order of their increasing end-to-end distances obtained at the end of heating simulations (at T = 500 K). Interestingly, the end-to-end distances, Mut2 (12 Å) < Mut3 (13 Å) ≈ WT (13.5 Å) < Mut1 (14.5 Å), were inversely correlated with the magnitude of the experimental chemotactic activities (Pazgier et al., 2007) (Table 1).

Discussion
In this study, we examined the dynamics and stability of HBD-1 WT and three single-site mutants. This enabled us to provide new insights into the structural differences between these species in order to understand their biological functions as antibacterial and chemotactic agents. Understanding the molecular basis for the antibacterial structure-activity relationships of HBD-1 WT and mutants might help in designing and developing novel antibiotics, as well as help in preventing the development of bacterial resistance to current antibiotics. This initiative has already been taken by Polymedix to develop nonpeptide antibiotics that resemble β-defensins (http:// polymedix.com/pdf/polymedix-antibiotic-whitepaperSept-08.pdf). Understanding the molecular basis for HBD-1 WT and mutants' chemotactic activity might also aid in discovering potential treatments for pathologies involving the immune system.
When overlaid with the WT monomer, the X-ray structures of all mutants studied had the same overall structures and secondary structures as the native WT; the only structural differences were in the turns (Pazgier et al., 2007). However, dynamic structural changes could result from the specific mutations that lead to the measured differences in biological activity. The investigators who created and reported the mutants' structures (Pazgier et al., 2007) purposefully left unaltered the conserved residues such as the cysteines (disulfide positions), GXC motif (Hoover et al., 2001), and glycine and proline residues (Hoover et al., 2000), due to their critical roles in determining the defen-sins' 3D structure and function. The β-bulge (Hoover et al., 2001) was conserved among WT and all mutants, where its first residue (number 24) adopted α-helical parameters, e.g. dihedral angles Φ and Ψ in the Ramachandran plots as reported for classical β-bulges (Chan, Hutchinson, Harris, & Thornton, 1993). Our Ramachandran plot showed that for residue 24 of all monomers, Φ ranged from À94 to À85 and Ψ from À57 to À44, i.e. in the range for classical α-helices (data not shown).
The dimers studied were all asymmetric, as reported for the PDB structure for the WT peptide (Hoover et al., 2001) (Figure 1(b)), where interaction occurs through the β2-strand of chain A and the β2-β3 turn of chain B. This asymmetry is also reflected in the time evolution of the dimer's secondary structure propensity (Figure 3(b)), side chain contact maps (Figure 9(a)), RMSF (Figure 9(b)), and SASA (data not shown). A result of this asymmetry is that all the dimers' individual monomer chains behave differently when compared to the stand-alone monomers. This is because the two chains in a dimer interact with a different portion of their sequence and each monomer affects the other's dynamics. An example of this is clearly seen when comparing the time evolution of secondary structure propensity for Mut2 monomer (Figure 3 (a)) to its dimer (Figure 3(b)). Here, in order of α-helix stability, Mut2 monomer > chain A of Mut2 dimer > chain B of Mut2 dimer; Mut2 monomer shows less π-helix character within the α-helix region (Figure 3).
The initial monomer and dimer structures together with those observed in equilibrium simulations at T = 300 K had their positively charged side chains surface exposed and grouped into separated patches of cationic and hydrophobic residues. This was clearly observed, e.g. in the SASA per residue plots and could be visually observed in VMD, looking similar to the WT equilibrated structure illustrated in Figure 1(d). The positively charged groups on the surface of the molecule have been established to be an important feature for maintaining chemotactic activity (Pazgier et al., 2007). It is thought that the spatially separated cationic and hydrophobic regions allow defensins to insert into the bacterial membranes. The basic residue patches are thought to interact with the anionic phospholipids allowing an oligomeric defensin structure to rearrange, insert hydrophobically into the membrane, and create a lethal hole in the bacterial membrane, causing antibacterial activity (Bauer et al., 2008;Hoover et al., 2000). Comparing the equilibrated WT monomer (Figure 1(c)) with its surface representation (Figure 1(d)), the five basic residues (Arg29, and Lys22, 31, 33, and 36) are observed to be more concentrated towards the C-terminus. These residues are present at the β3-strand, and at the end of β1-β2 and β2-β3 turns. On the other hand, the hydrophobic cluster is more concentrated at the β1-β2 turn. One hydrophobic residue is present per secondary structure element: Val6, Leu13, Ile23, and Ala32 in the α-helix, β1-, β2-, and β3strands, respectively. These general features of 3D structure described for the WT monomer were also true for all the mutants used in our study due to their similar 3D structures (Pazgier et al., 2007).
The three conserved disulfide bonds play an important role in maintaining the tertiary structure motif of defensins. Any effect of the specific mutations on the disulfide bond dependent dynamic properties therefore is a legitimate concern. During all simulations, the three disulfide bond lengths showed minimum change at ambient and elevated temperatures. This is not surprising and highlights the crucial role of disulfide bonds in maintaining a correct constrained tertiary structure backbone. Disulfide bond based tertiary structure is thought to exist in order to provide defensins with resistance to proteolytic degradation in the environments where they are secreted and active (Selsted & Ouellette, 2005). Despite the presence of the disulfide bonds, heating the molecules resulted in their swelling. This was reflected in the decrease in the number of intramonomer hydrogen bonds and an increase in the radius of gyration, as expected for heating simulations (Daggett & Levitt, 1992. Whenever the α-helix lost structure, as in the case of WT and Mut3 monomers, this occurred towards the end of the trajectory and always in regions away from the disulfide bonds. For example, the WT α-helix was partially unfolded from Asp1 through Cys5 at T = 458 K ( Figure 2).
The single-site mutations were found to affect the molecular properties. During the equilibrium simulations of Mut1 (Asn4Ala) monomer, residue 4 had on average a lower SASA than the WT, which is due to the shorter side chain of Ala compared to Asn. In Mut1, residue 4 also formed a strong intramonomeric contact with Gln24, that was absent in the WT. This could be attributed to the difference in side chain polarity; Asn is polar while Ala is not. The Asn4Ala mutation also caused Mut1 dimer to form more stable contacts, as evidenced from equilibrium simulation at T = 300 K, compared to the WT, giving it a greater overall stability. During all equilibrium simulations for Mut1 dimer, the first 10 residues exhibited lower RMSF than those of the WT. This is due to the shorter side chain of Ala vs. Asn causing less fluctuation. Heating the Mut1 (Asn4Ala) monomer from T in = 300 K to T fin = 500 K resulted in a more stable α-helix than that found in the WT (Figure 2), suggesting that at higher temperature this mutation stabilized the αhelix. This stability is likely due to either the hydrophobic nature of the Ala substitution, or because Ala is known to overstabilize the helix (Tanizaki, Clifford, Connelly, & Feig, 2008).
In the case of equilibrium simulations for Mut2 (Gln24Glu) monomer, because Gln and Glu amino acids have almost the same chain length, there was no significant difference compared to the WT regarding their average RMSF or SASA at residue 24 throughout the simulation. Upon heating, the Gln24Glu mutation did not only affect Mut2 β2-strand structure stability (wherein residue 24 lies), but also affected the α-helix. This is because these two secondary structure elements interact with each other. All equilibrium simulations for Mut2 (Gln24Glu) dimer showed the Glu24A/ Thr26B stable contact (Figure 7(b)) that was absent in the WT dimer. Other new stable contacts, absent in all other HBD-1 dimers tested due to Gln24A mutation, included Glu24A/Tyr3B, Glu24A/Asn4B, and Glu24A/ Arg29B. The Glu24A/Arg29B contact involved the negatively charged side chain of Glu interacting with the cationic Arg. The rest of the new contacts were established due to hydrogen bond formation involving the terminal carboxylate group of Glu or due to the resulting difference in contacts between the WT and Mut2 species.
Equilibrium simulations for Mut3 (Lys31Ala) monomer demonstrated lower SASA at the mutation site compared to the WT, due to the shorter side chain of Ala compared to Lys. Heating Mut3 resulted in distortions in both β1and β3-strands, because they interact directly with each other in a manner similar to the Mut2 discussed earlier. The Mut3 contacts between Ala31 and both Tyr14 and Ser15 were weaker than those contacts formed by Lys31 in the case of the WT. This was due to the side chain of Ala being too short to maintain the contact during the simulation run. It could also be due to the polar nature of Tyr, Ser, and Lys, making these residues incapable of interacting with the nonpolar Ala. Ala31 contacts in Mut3 dimer were also relatively weak compared to Lys31 in the other tested mutant dimers. Some of the Mut3 trajectories showed a Thr26A/Ala31B contact (Figure 7(d)), while one had a transiently stable Ile23A/Ala31B contact and a strong Thr26A/Ala31B contact. This again could be attributed to the short Ala side chain that cannot reach more residues in the other chain.
Analysis of the time evolution of secondary structure propensity agreed well with the results from the structure overlap function, propensities, and Ramachandran plots, where the monomers showed almost equally stable αhelices throughout the course of their equilibrium simulations at T = 300 K (Figure 2). Heating the molecules from T in = 300 K to T fin = 500 K caused some noticeable changes in the secondary structure ( Figure 2). The αhelix was unstructured at T = 460 K in the WT, at T = 500 K in the Mut1, transiently in the 476-488 K range in the Mut2, and at T = 475 K in the Mut3 (Figure 2). On the other hand, β-strands showed only minor changes compared to the extent of alterations in the αhelices. Specifically, the β1-strand was transiently unstructured between T = 430-450 K in Mut1 and at T = 460 K in Mut3, while the β2-strand was unstructured at T = 470 K in the case of Mut2.
During heating simulations, the loss of secondary structure elements was reflected in various analysis methods. Temporal changes in the end-to-end distance were directly correlated with the secondary structure dynamic changes. When the molecules started to unravel, their two ends moved away from each other, thereby increasing the end-to-end distance, accompanied by an increase in SASA (see, e.g. Figure S2). The gradual loss of Mut3 monomer side chain contacts during heating from T in = 300 K to T fin = 500 K (Figure 6) is also a result of the loss of secondary structure elements. Ramachandran and propensity plots (see Figure 4) reflected secondary structure changes throughout the simulation runs. All these secondary structure changes also coincided with decreases in the structure overlap function (Figure 5), and showed that the order of monomer stability was Mut1 > WT > Mut2 > Mut3. This order was only slightly different from the α-helix dynamics results (Figure 2 right panel) having a decreasing order Mut2 > Mut1 > Mut3 > WT. This minor difference likely results from the fact that the χ results involved differences between the entire structure, involving all secondary structure elements' dynamics, rather than just the α-helix dynamics. The heating simulation for Mut1 was the most stable trajectory, possessing an almost intact α-helix (only deformed at T = 500 K) and a transient β1-strand deformation between 430 and 450 K. This was followed by WT that exhibited an unstructured α-helix at T = 460 K, followed by Mut2 having β2-strand and transient α-helix destructuring at T = 470 K and T = 475-490 K, respectively. Then, Mut3 was least stable due to its unstructured β1-strand and α-helix at T = 460 K and T = 475 K, respectively, that both persisted until the end of the trajectory.
In a previously reported MD simulation study (Sharadadevi & Nagaraj, 2010), the authors observed an early α-helix unfolding with total loss of hydrogen bonds at T = 300 K. We did not observe α-helix unfolding until $22 ns of the WT monomer equilibrium simulation at the same temperature ( Figure S3(a)), and no total loss of hydrogen bonds was detected. Our MD simulations also suggest that different force fields used in simulations could affect the output results and thus should be carefully chosen. The outcome of using different force fields in our simulations suggested that the difference between the αand π-helix structures is so minor that it could be neglected. This makes sense when one considers the minor difference between the dihedral angles observed for the two helix types.
Because we aimed at investigating HBD-1 antibacterial and chemotactic activities, all the analyzed terms were related to these biological functions whenever possible. All the contacts formed by the basic residues were monitored as we expected them to have minimal intra and/or intermolecular interactions. This feature is consis-tent with the low number of acidic residues found in the sequence capable of favorable charge interactions. Free patches of basic residues are likely a requirement for the defensins performing their antibacterial and chemotactic activities (García, Jaumann, et al., 2001;Klüver et al., 2005;Krishnakumari et al., 2003;Mandal et al., 2002;Pazgier et al., 2007). For all monomer equilibrium simulations at T = 300 K, only Lys33 and Lys36 exhibited cumulative side chain contacts throughout the simulation runs, while Lys22, Arg29, and Lys31 formed almost none. This would potentially make them free to bind to sites on the CCR6 receptor, allowing HBD-1 to exert its chemotactic effect, and also make these positive charges available to carry out the antibacterial effect, whose effectiveness depends on the molecule's total surface positive charge (García, Jaumann, et al., 2001;Klüver et al., 2005;Krishnakumari et al., 2003;Mandal et al., 2002).
Increasing the molecule's surface positive charge density via oligomerization is believed to increase the defensins' antibacterial effect at the bacterial membrane (Suresh & Verma, 2006). Although the WT has mostly been reported experimentally as a monomer (Hoover et al., 2001;Pazgier et al., 2007;Schibli et al., 2002), it was also crystallized as a dimer (Hoover et al., 2001). The single-site mutants used in our study were only reported crystallographically as monomers (Pazgier et al., 2007). However, we were able to demonstrate the potential multiple monomer interaction sites for the Mut2 and Mut3 species. This hypothesis was based on the observed sliding behavior without any dimer dissociation. This sliding phenomenon could help promote the oligomerization of the tested mutants. One of the equilibrium simulation runs for Mut2 dimer showed a stepwise decrease in χ at 6, 18, 22, and 24 ns (Figure 7(a)) indicating loss of similarity with the initial PDB structure. This was accompanied by a decrease in the total number of intermonomer contacts (Figure 8 solid black), and an increase in RMSD (Figure 10(b) solid black). Its dynamic maps of stable contacts (Figure 7(b)) revealed weaker contacts together with a variation in the interaction site switching from interaction between chain A αhelix with chain B β2-strand to a new interaction between chain A β2-strand and chain B β2-β3 turn. At 6 ns some contacts between α-helix and β2-strand were lost, then at 22 ns the interaction site weakened where only Arg29B formed three stable contacts with Lys22-Ile23-Gln24 (Figure 7(b)). The weakest Mut2 dimer trajectory ( Figure 8) lost contacts at 6 ns, totally unbound at $22 ns, then redimerized via the formation of a few contacts at $25 ns. This change in interaction sites for the HBD-1 Mut2 and Mut3 dimer species (intermolecular sliding behavior) was manifested in the χ (Figures 7  (a) and 6(c)), dynamic maps of stable contacts (Figures 7 (b) and 6(d)), and in the intermonomer contacts plot (Figure 8).
This behavior suggests that HBD-1 molecules have multiple potential intermonomer interacting sites in vivo, and this could be a biological strategy that evolution has devised to aid the molecules in performing their functions as antibacterial and chemotactic agents. Perhaps the sliding behavior enables HBD-1 to oligomerize at the bacterial membrane, increasing their overall surface positive charge density, thus enhancing their antibacterial effect, as hypothesized in an earlier modeling study (Suresh & Verma, 2006). It could also facilitate HBD-1 binding to the CCR6 receptor site via sliding dynamics that 'finds' the correct dimer conformation in order to bind, thereby initiating chemotactic activity. If the active species initiating chemotaxis is a monomer, and dimer vs. monomer as the active species is not a settled issue in the literature (Chan et al., 2007), then its ability to exhibit sliding behavior with CCR6 within the binding region may be the mechanism whereby receptor binding occurs. The observed sliding behavior of Mut2 and Mut3, but not WT and Mut1, could explain their higher chemotactic activity compared to WT and Mut1, where the reported measured chemotactic activities were as follows: Mut2 > Mut3 = WT > Mut1 (Pazgier et al., 2007) (Table 1). Longer trajectories need to be generated to see whether WT and Mut1 possess this sliding behavior; being more stable than Mut2 and Mut3 dimers does not mean that they lack this trait because this could be attributed to insufficient sampling time. If WT and Mut1 dimers were proved to lack this behavior, this could mean that not only the mutation possibly caused this sliding phenomenon, but also that the site of mutation is crucial for whether sliding occurs in the dimer. Mut1 has the mutation site at the α-helix, while Mut2 has the mutation at the β2-strand and does not seem to be as stable as Mut3, whose mutation occurs at the end of the β2-β3 turn.
Comparing the HBD-1 monomers reported chemotactic activity (Pazgier et al., 2007) to the α-helix stability observed in equilibrium simulations at T = 300 K, we concluded that the chemotactic activity is directly proportional to their α-helix stability ( Table 1). The α-helix has been reported to be crucial for binding to CCR6 and therefore important for chemotaxis (Hoover et al., 2002;Pérez-Cañadillas et al., 2001). This is undoubtedly the reason for the low chemotactic activity of Mut1, since it is mutated within the α-helix leading to a destabilized structure. We also observed that the dimer stability was inversely proportional to the chemotactic activity (Table 1). In terms of the total number of intermomomer contacts, Mut1 dimer was the most stable molecule followed by Mut3, then the WT, with Mut2 the least stable. This decrease in dimer stability was observed to correlate with a greater chemotactic activity with a slight difference between the WT and Mut3 ranking. Together, these results present an interesting model that implies HBD-1 monomers are more active as CCR6 ligands than the dimers in binding to CCR6, as is the case for MIP-3α (Chan et al., 2007).
The order of the calculated average χ values for monomer equilibrium simulations, Mut3 (0.96) > Mut2 (0.94) > WT (0.92) > Mut1 (0.9), was directly correlated to the experimental chemotactic activities with a slight difference in Mut2 and Mut3 rankings (Table 1). On the other hand, end-to-end distances for the same equilibrium simulations showed an ordered relationship opposite to that of the magnitude of the experimental chemotactic activities: Mut2 (12 Å) < Mut3 (13 Å) ≈ WT (13.5 Å) < Mut1 (14.5 Å) ( Table 1). These data suggest that the more compact the structure, the better the surface exposure/orientation of the side chains needed to bind to the CCR6 receptor. This correlation also supports the theory that the monomer, rather than the dimer, could be the active in vivo binding ligand for CCR6. This is as far as we can go in correlating our results to measured chemotatic responses for the following reason. Chemotoxis is a complex event phenotype involving many cellular recognition/signaling events, not all of which are known. Therefore, we cannot suggest interpretations beyond the simple direct comparison of rank effects of the WT and mutants' properties with measured chemotaxis outcomes.

Conclusion
In this study, we examined the dynamics and stability of HBD-1 WT and three single-site mutants (Mut1, Mut2 and Mut3) using all-atom MD simulations of their monomers and dimers and theoretical analysis. Our results shed some light on the underlying basis of the antibacterial and chemotactic properties of HBD-1 that is rarely studied in the literature. During all the equilibrium simulations at room temperature, the reported conserved structural features thought important for HBD-1 to perform their functions were maintained with minor or no changes. Among these are the GXC motif, glycine and proline residues, and the β-bulge. The results also confirmed the conservation of disulfide bonds and their role in maintaining the defensins' structural integrity where only minor length changes were observed even upon heating the molecules.
The single-site mutation types and positions affected the molecules' side chain contacts and stability, and accordingly their biological functions. HBD-1 monomers' α-helix dynamic stability was directly related to their reported experimental chemotactic activity (Pazgier et al., 2007), whereas the dimers' α-helix dynamic stability was inversely related (Table 1; see also Figure 2). There was an observed sliding behavior in the case of Mut2 and Mut3 that is thought to give them the higher chemotactic activity when compared to Mut1 and WT (Table 1). A possible explanation is that the sliding alters the dimer contacts and tertiary structure, and could aid in the HBD-1 oligomerization that potentially enhances their chemotactic activity. The same mechanism is thought to promote the HBD-1 antibacterial activity by enhancing the surface positive charge density. The importance of the surface positive charge to the molecule's biological activity was also mirrored by their effects on the structure and dynamic properties. All the basic amino acids were exposed to the surface as observed in the 3D structure representations (Figure 1(d)), SASA, and minimal side chain contacts for basic residues. The potential dependence of chemotactic activity on the monomer and dimer calculated average χ, and end-to-end distances (Table 1) suggests that closer resemblance to the native structure and a more compact 3D structure results in better surface exposure and orientation to bind to the CCR6 receptor. Our results also suggest that the monomer, rather than the dimer, is the active in vivo chemotactic agent. Other structural observations proved the α-helix to be less thermally stable than the β-sheet. The dimer structure asymmetry was seen in 3D representations (Figure 9(a)), time evolution of the secondary structure propensity (Figure 3(b)), side chain contact maps (Figure 9(a)), RMSF (Figure 9 (b)), and SASA (data not shown).
Understanding the molecular basis for the SAR of HBD-1 could help in the process of drug discovery of novel drugs for targeting the immune and pathogen defense systems, while preventing any potential development of drug resistance.