Binding of Copper Ions with Octapeptide Region in Prion Protein

Article Cite This: J. Phys. Chem. B 2019, 123, 5216−5228

pubs.acs.org/JPCB

Binding of Copper Ions with Octapeptide Region in Prion Protein: Simulations with Charge Transfer Model Ke Chen, Wenfei Li, Jun Wang,* and Wei Wang* National Laboratory of Solid State Microstructure, Collaborative Innovation Center of Advanced Microstructures, and School of Physics, Nanjing University, Nanjing 210093, P.R. China

Downloaded via NOTTINGHAM TRENT UNIV on July 22, 2019 at 02:11:31 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Copper ions are important cofactors of many metalloproteins. The binding dynamics of proteins to the copper ion is important for biological functions but is less understood at the microscopic level. What are the key factors determining the recognition and the stabilization of the copper ion during the binding? Our work investigates the binding dynamics of the copper ion with a simple system (the N-terminus of PrP) using simulation methods. To precisely characterize the protein−ion interaction, we build up an effective copper−peptide force field based on quantum chemistry calculations. In our model, the effects of charge transfer, protonation/deprotonation, and induced polarization are considered. With this force field, we successfully characterize the local structures and the complex interactions of the octapeptide around the copper ion. Furthermore, using an enhanced sampling method, the binding/unbinding processes of the copper ion with the octapeptide are simulated. Free-energy landscapes are generated in consequence, and multiple binding pathways are characterized. It is observed that various native ligands contribute differently to the binding processes. Some residues are related to the capture of the ion (behaving like “arm”s), and some others contribute to the stabilization of the coordination structure (acting like “core”s). These different interactions induce various pathways. Besides, a nonnative binding ligand is determined, and it has essential contributions and modulations to the binding pathways. With all these results, the picture of copper−octapeptide binding is outlined. These features are believed to happen in many ion−peptide interactions, such as the cooperative stabilization between the coordinations with neighboring backbone nitrogens and an auxiliary intermediate coordination with the neighboring oxygen from the N-terminal direction. We believe that our studies are valuable to understand the complicated ion−peptide binding processes.

■

INTRODUCTION Copper ions (both the reduced and oxidized forms) are involved in complicated metabolism in biological systems and are important cofactors of many metalloenzymes.1 It is observed that the intracellular concentration of copper(II) ions is much lower than the extracellular one,2 which ensures that copper ions are mostly bound to proteins.3,4 The binding modes, processes, and dynamics would be important to understand the functions of copper ions together with the proteins. For example, there are complicated binding modes between copper ions and octapeptides of prion proteins (PrPs).5−15 The switching between various modes is essentially related to the binding dynamics of copper ions. It is observed experimentally that the interaction of the N-terminus with copper ions can trigger PrP to misfold into a highly proteaseresistant form with enhanced β-sheet contents.16−20 It is suggested that the His96−Cu2+−His111 cross-link takes an important role in the formation of the misfolded form.21,22 Besides, it is found that the binding of copper ions can trigger the endocytosis of PrP23,24 and can block effectively the potential toxicity of PrP N-terminus.25,26 Quantitative © 2019 American Chemical Society

characterization of copper binding with the octapeptide regions may help to understand the behaviors of PrP. In metalloproteins, there are complicated physical interactions, including hydrophobic interactions, van der Waals and coordination interactions, and so on. There are a wide range of strengths for these interactions, which introduce diverse time scales in binding dynamics, related to the degrees of freedom of electrons, ions, and peptide chains. Besides, because of the connectivity and competition in the peptides, the coordination to the related ligands might not be independent random processes. This is one of the physical reasons to make the binding dynamics difficult to study. How to characterize the binding patterns, the binding geometry, and the binding dynamics is the basic question to study the interaction between copper ions and proteins. To understand the binding, a model system with simple interactions as well as an efficient and precise implementation for the concerned interactions are Received: March 15, 2019 Revised: May 31, 2019 Published: June 4, 2019 5216

DOI: 10.1021/acs.jpcb.9b02457 J. Phys. Chem. B 2019, 123, 5216−5228

Article

The Journal of Physical Chemistry B expected. In present work, we employ the Cu-HGGGW complex as the model system, since only coordination with copper ion and peptide are involved. On the basis of the crystal structure of Cu-HGGGW complex, we build up a force field for copper−peptide interactions which mimic all major effects (such as charge transfer, polarization, deprotonation, and so on) related to coordination bonds. This force field enables us to use classical molecular dynamics to simulate the binding dynamics efficiently. With replica exchange molecular dynamics (so-called REMD),27,28 the landscape of the model system is extensively sampled. The trajectories of a total of 8.4 μs are obtained. The free-energy landscape of the binding process is derived in consequence. It is observed that the formation of the coordination bonds guides the binding processes. Together with several constant-temperature simulations, the formation order of the coordination bonds are analyzed. Different roles of various ligands are observed. To further understand the binding processes, the transition network between conformation clusters (characterized by their coordination bonds) is generated. Multiple pathways and mechanisms are observed. The effects of ligands during various stages of binding processes are discussed. It is found that the native coordinations with two middle residues are responsible for the stability of the system, while the ligands in two terminal residues help the capture of the ion at the early stage of the binding processes. Additionally, a nonnative coordination is observed to be essential for the binding processes, stabilizing the system at the early stage and modulating the pathways of the binding pathways. As a conclusion, through the analysis on this simple system, several general features for metalloproteins are suggested. This may help us to understand the contributions and behaviors of various ligands during coordination.

Figure 1. Copper−octapeptide 1:1 binding mode with structure of the Cu-HGGGW motif determined by X-ray crystallography14 (top) and the amino acid sequence of human prion protein 21−124 with the 4-fold octarepeat region highlighted and residues contributing coordinating atoms colored (bottom). The color of G3O is omitted in the sequence as it shares the residue with G3N.

main interactions in the metalloproteins. These kinds of interactions are generally strong. A long time period is necessary to observe the binding/unbinding transition. This kind of interaction is generally related to the electron degrees of freedom, which requires quantum mechanics (QM) to characterize. The long time scale and the consequent large computational demand prohibits large-scale simulations of the structural variations of ion−peptide systems with QM. The calculation based on QM becomes the bottleneck of the simulations. Considering the large noises from the solvent environment as well as the approximations related to other interactions, high precision description for coordination interaction based on QM may be replaced with a classical force-field. A proper classical force field may help to provide an efficient characterization with acceptable precision. In our work, we derived a practical force field between copper ion and peptides. How can we build a proper force field to describe the interactions between the copper ion and the related ligands? The elimination of the degrees of freedom of electrons would introduce some many-body effects into the effective force field. Thus, the coordination interactions between copper(II) ions and proteins are more complicated compared with the regular electrostatic interactions or van der Waals interactions. For example, the polarization effect is pronounced in copper(II)containing systems since the highly charged ion induces highly polarized electron cloud.31 Rick and Stuart indicated that the polarization effect is not additive when there are more than one ligand in the system,32 which indicates that the polarization (and charge delocalization especially) of copper ions and proteins cannot be neglected to quantitatively characterize such a system.33 Besides, there are many other microscopic processes involved, such as (de)protonation, charge redistribution of imidazole ring, and so on.34 These all require more detailed models of charge distributions in the system. As first step, in the force field related to ion−peptide systems, the electrostatic interaction is undoubtedly fundamental. Besides, considering the contribution of the electron degrees of freedom, some additional terms are necessary to

■

METHODS Copper−Peptide System. In our work, the interaction between copper(II) ion and a peptide with the sequence PHGGGWGQ is considered. This peptide is one octarepeat segment in PrP. There is a specific binding structure between the copper ion and the peptide. The crystallographic structure of the complex Cu-HGGGW is experimentally determined.14 In this structure, the copper ion is bound equatorially by Nδ of histidine side chain (H1N), two deprotonated N atoms from the next two glycines (G2N and G3N), and the carbonyl O of the second glycine (G3O), while water molecules take one or both of the axial binding positions, see Figure 1. The structural data are given in Supporting Information, which is provided by the Millhauser lab. The crystal structure is used as a basic constraint for parameter optimization and model evaluation. Indeed, there are several other binding modes related to the coordinations to multiple octapeptides. These binding modes are often postulated on the basis of local coordination information on copper ions.13,15,29,30 Because of the complexity of interactions and structures, only the structure of the system with an octapeptide and a copper ion is determined presently. In this sense, we pick the Cu-HGGGW complex as the model system. It is worth noting that all kinds of ligands (such as nitrogen ligands, oxygen donors, and histidine imidazoles) are all involved in the model system. Thus, the results for the model system may help to understand other binding modes between octapeptides and copper ions. Interactions for Copper−Peptide System. The coordination interactions between ions and peptides are one of the 5217


Article

The Journal of Physical Chemistry B

Correspondingly, a quadratic form is employed to characterize the transferred charges between copper ion and oxygen in water molecule, as ΔqCu−i = Cir2Cu−i + AirCu−i + Bi. Here Ci, Ai, and Bi are concerned parameters, and can be determined by fitting the results of quantum chemical calculations. It is worth noting that these ligands are calculated separately on the basis of model systems, and the coupling between these ligands during coordination is ignored in this formalism. Following the modification for the electrostatic interactions, the parameters for the van der Waals interactions of copper ion are also modified (different from the assignments in the original AMBER ff03 force field). In detail, the parameters for van der Waals interactions, ϵCu−i and σCu−i, are taken from the literature39 to reproduce the experimental hydration free energy of copper(II) ions, with the values ϵCu−i = 0.0734 kcal/ mol and σCu−i = 1.55 Å. Another effect related to the binding of copper ions is the redistribution of the charges in the related imidazole ring because the electrons in the imidazole ring are highly delocalized. A significant polarization of the imidazole ring may be observed in the case with the binding between copper ion and histidine. Clearly, this effect cannot be depicted by the AMBER ff03 force field. To mimic this effect, we employ a distance-dependent modification for the charges of the ring atoms. When the copper ion goes close to one of the N atoms in the side-chain imidazole ring, the charge of the concerned N atom increases linearly with the distance, and the charges of the other atoms in the ring are accordingly scaled so that the total charge of the imidazole ring is unchanged (nearly zero). This modification describes the redistribution of charge in the ring, and produces a dipole in the ring, which mimics the copper-induced ring polarization. Quantitatively, the amount of charge increase on the coordinating N atom is determined by fitting the experimental coordinating geometry. Similar to the charge-transfer effect, the large-distance limit should indicate the vanishing of the polarization, and the interaction would converge to the original force field. In our case, we assume that the convergence to the force field AMBER ff03 is achieved when the distance is twice as large as the summation of VdW radii of copper ion and N atom. Note that there are two kinds of binding sites in the imidazole ring, that is, Nδ or Nϵ atom, especially in the case with low copper concentration.40 In our work, the parameters for Nδ are derived. As a remark, in our case, only the polarization of the imidazole ring in histidine is considered because there are no delocalized electrons for other kinds of residues. Besides the electron degree of freedom, the inclusion of copper ion may alter the protonation of certain ligands in peptides. For example, in the copper−octapeptide 1:1 binding mode, two backbone N atoms of glycine residues are deprotonated. In the state without coordination, these N atoms are generally protonated. The binding with copper ion would produce a transfer of proton from peptide chain to solvent. Complicating this fact, the protonation depends on other interactions between various units of peptide. It is reported that the histidine−aromatic interaction, such as the His−Trp interaction, can raise the pKa of histidine imidazole and contributes to protein stability.41,42 On the basis of potentiometric measurements, the protonation constant of the imidazole side chain in Ac-PHGGGWGQ-NH 2 is determined as pKa of 6.2−6.4.43 This result indicates that the side chain of the histidine in the peptide is partially protonated at bioactive pH condition. All these exemplify that the charge largely varies

characterize the stability and local geometry. This is important for the transitional metallic ions (such as copper) because of the large variations of their electron clouds. Motivated from previous studies,35,36 some special terms are considered for the interactions between copper(II) ion and the peptides to describe the effects of charge transfer, protonation/deprotonation, and polarization of units in proteins. They are described in the following. It is worth noting that we expect the original force fields for peptides and solvents are minimally modified after the inclusion of copper ions. This kind of setup would not introduce any unexpected biases for the interaction and dynamics of peptide and solvents. In this work, only the terms related to copper ions are modified. The interactions between all other atoms (including those inside the peptides and between peptide and solvent molecules) are characterized with the typical force field (AMBER ff0337). This kind of setup matches our expectations. Considering the effect of charge transfer, variable charges are employed in electrostatic interactions related to copper(II) ions. That is, the (nonbonded) interactions between copper ion and the ligand i is depicted as É ÅÄÅ 12 6Ñ ÅÅij σ qCuqi yz ij σCu − i yz ÑÑÑÑ Cu − i Å zz − jj zz ÑÑ + 4ϵCu − iÅÅÅjjj VCu − i = jj r zz Ñ ÅÅjk rCu − i zz{ 4π ϵ0rCu − i k Cu − i { ÑÑÑÖ (1) ÅÇ in which the charges qCu and qi are defined as the functions of the distance between copper ion and ligand i, rCu−i. When the distance is large, the charges are equal to the quantities of the isolated ion and ligand. While the relative distance between these units approaches that in the native conformation, that is, these units are rather close, there is an apparent charge transfer between these units. The corresponding charges would be apparently different from the isolated case. To discard the high-order momentums as well as the dynamic effects, the interaction can be described as the electrostatic interactions with distance-dependent charges. This kind of force field describes the charge transfer between copper ion and ligands in a quasi-static manner and has been successfully applied in the simulation of zinc-ion coupled folding of Cys2His2 zinc-finger protein.35 Note that this format of interaction is similar to electrostatic interaction in the original force field (Amber ff03) and thus eases the implementations. The detailed parametrization is similar to the previous implementation for zinc ions based on quantum chemistry calculations. As a first step, for a combination of ligands, quantum chemical geometry optimizations at the B3LYP/631+G* level with MDF10 pseudo potential applied to copper38 are performed on the systems (copper and its ligand atoms), and the details for the systems are shown in Supporting Information. After the optimization, the copper− ligand distance and the natural bond orbital (NBO) charge distribution of the concerned system are determined for further processing of transferred charges. To incorporate these data into a force field for molecular dynamics, the distance dependencies of transferred charges are fitted with simplified functions in practice. For the ligands Nδ of histidine, backbone O, backbone N, and side-chain O of glutamine, the distance dependence of charges is assumed to be linear, as ΔqCu−i = AirCu−i + Bi. For the ligand from water molecules, the binding mode is a little different. The water molecules are likely to occupy either equatorial or axial position of the copper ion, which are not equivalent due to Jahn−Teller effect. 5218


Article


Simulation Setup. In this work, the copper−octapeptide complex is investigated on the basis of the above modified force field. The complex is solvated by 2320 TIP3P waters in a 46 Å × 40 Å × 39.5 Å box. To mimic the solvent buffer, three NaCl molecules are added, which corresponds to the concentration of 51 mM. During protonation, the charges of Na ions decrease to balance the charge. To calculate the thermodynamic properties, REMD simulations are carried out. In total, 36 replicas with their temperatures ranging from 275 K to 490 K27,28,45−47 are employed. For each replica, a 1 ns simulation is carried out for heating and equilibrating. After that, a trajectory with a length of 114.3 ns is simulated for each replica. It is found that after the evolution of the first 10 ns, diverse configurations of coordination can be observed between replicas. It is reasonable to skip the first 10 ns trajectory to remove the effects of initial setup. The later part with the length of 104 ns is used to calculate thermodynamic properties. Distance restraints are applied between copper and native ligands (a copper−ligand distance smaller than 6 Å for H1N and G3O and 5 Å for G2N and G3N) to reduce the translational entropy during the simulations and to accelerate the simulations. These upper limits for the distances between the copper ion and ligands are large enough to permit access to the fully unbounded state of the peptide. Additionally, to investigate the dynamic processes, 260 unbinding trajectories with the lengths of 20 ns at 350 K are calculated, with a 30 ps heating and a 20 ps restrained equilibration being carried out for each unbinding trajectory. The starting structure of these simulations is prepared on the basis of the local geometry of Cu-HGGGW binding site determined by X-ray experiment.14 The structure of other residues are obtained through energy minimization. Trajectories are run with AMBER11 package with the concerned charges modified on-the-fly. It is worth pointing out that the combination of AMBER package and the charges modified onthe-fly may break the detailed balance of the dynamics because of the omission of the many-body forces related to the contributions from the variation of charges. Yet, for our cases, the deviation is acceptable because the variation of charges is small and slow (see further discussions in Supporting Information) and will not deteriorate our conclusions. The Modular Reweighting package,48 an implementation of the weighted histogram analysis method (WHAM),49 is employed in the calculation of the free energy landscape. The molecular structures are visualized in PyMOL.50

related to protonation/deprotonation state. Indeed, this kind of processes cannot be described in the original AMBER force field. Rather than introducing the protons explicitly, we describe the protonation/deprotonation effects in an implicit way. Considering the hydrogen atom is small and largely mobile, the exclusive volume of the hydrogen atom can often be omitted. The only effect of hydrogen atom is its charge. In this sense, the differences between protonated and deprotonated states can be described with the charge quantity and distribution of the concerned residue. On the basis of this idea, we omit the hydrogen atom (proton) in the protonated state, and the charge of hydrogen atom is merged with that of the concerned N atom and the dipole is discarded. For the case with partial protonation, the partial charge of the proton is added to the residue. In our work, a partial protonation charge of 0.5e is assigned to the histidine to describe the solvated condition without the binding of the copper ion. This protonation corresponds to pH about 6.4, which is within the bioactive pH range. Copper−octapeptide complex is well observed at this pH by both potentiometric measurements and Raman spectroscopy.43,44 Differently, there is no existing model for deprotonated glycine. To determine the corresponding charge distribution, quantum chemical optimizations at the B3LYP/6-31+G* level are employed for both protonated and deprotonated glycines. The NBO charge distributions for both states are derived. The difference between the two distributions is applied to the ff03 glycine model to derive the charge distribution of deprotonated glycine. The obtained charges are fine-tuned to reproduce the experimental copper− O distance in the energy-minimized structure before being used in this work. With these two states on hand, the charge variation of atoms during the protonation/deprotonation process is approximated as a linear interpolation related to the distance between the copper ion and the N atom. It is assumed that the system adopts the solvated state when the distance dCu, N is larger than twice the sum of the VdW radii of copper ion and N atom and is in the fully deprotonated state when the distance dCu, N is not larger than the equilibrium distance. It is worth noting that there are many other polarizable force fields in practice (see the review by Li and Merz33 and the references therein). These force fields are also designed to characterize variations of charge distribution coupled with the evolution of conformations. Compared with the regular polarizable force fields, our model has a different design philosophy. In our model, there are no explicit definitions of additional degrees of freedom. Several essential states are determined (such as the fully coordinated state and the unbounded state in present case). The charge variations are extrapolated on the basis of these presumed states. Clearly, our method requires sufficient numbers of input states, which is often available for small and simple systems. Once these states are assigned, the charge distribution can be determined easily and precisely with our model by extrapolation. Various kinds of effects (such as protonation/deprotonation) are all included, and the extrapolations are generally defined locally on the landscape of charge. Besides, some computation-demanding self-consistent calculations are often needed for regular polarizable force fields,33 which limits their application in long-time simulations. Our model is simple and efficient in computation and is suitable for the coordination interactions with clear and unique energetic minimum.

■

RESULTS AND DISCUSSION Validity of the Modified Force Field Related to Copper Ion. On the basis of the idea described in Methods section, the force field related to the copper ion is determined. The related parameters are given in Supporting Information. This force field can describe the charge variation during coordination processes with acceptable precision. Can this force field describe the coordination geometry around the copper ion? This is not a trivial problem, since the parameters of the force field are generated on the basis of simplified systems and the validity to describe the heterogeneity around copper ions still requires investigations. Faced with this problem, several tests are carried out, including solvation structure of copper ion, the coordination with multiple histidines, and the coordination geometry around the copper ion in the copper−octapeptide system. In all of these systems, compared with the quantum chemistry calculation or

5219


Article


calculate this distance with DFT optimization using the B3LYP/6-31+G* and CPCM implicit solvation model, which is believed to be more precise compared with the classical force field. In the DFT calculation, the equatorial Cu−Ow distance is in the range of 2.017−2.019 Å (as indicated by the blue line in Figure 2a). Clearly, the result based on ModFF is more comparable with the DFT result. It implies that ModFF is better than FF03 in the quantitative aspect. Additionally, the distance was reported previously as 2.02 Å in B3LYP-based QM/MM MD simulations on the [Cu(H2O)6]2+51,52 and 1.95 Å to 2.04 Å in EXAFS measurements.53,54 All these results further support the rationality of the ModFF. Another typical coordination structure in the copper− peptide system is the interaction between copper ion and histidines. Generally, there are multiple binding sites for histidines, including Nδ and Nϵ in the imidazole ring. Here, two kinds of models are considered in the comparison between ModFF and FF03. They are the coordinations between the copper ion and (1) one histidine protonated only at its Nδ site (Cu-HIE) and (2) four histidines protonated at their Nδ sites (Cu−4HIE). The Cu-HIE binding can be observed in the 1:1 copper−octapeptide coordination mode,14 and the Cu-4HIE model is one of the putative configurations of the 1:4 copper− octapeptide coordination mode but has been reported to be less favored at the high copper concentration.40 For these models, the initial structures are generated on the basis of the configurations of corresponding DFT calculations.40 Starting from these structures, a 2000-step minimization, a 50 ps heating simulation, and a 10 ns equilibrium simulation are carried out sequentially. The structural characteristics of these models are analyzed on the basis of these simulations, and the related results are listed in Table 1 For these two models, both coordination numbers are six: one Nδ and five waters for Cu-HIE and four Nδ and two waters for Cu-4HIE. Note that the site Nδ is deprotonated after coordination and the site Nϵ is still protonated. This is consistent with general knowledge about the coordination with copper ions. There are quantitative differences in their geometry and stability for the force fields ModFF and FF03. First, the bond length dCu, N determined by ModFF is smaller than that by FF03. This is true in both models. Especially, for the Cu-4HIE system, the bond lengths from the ModFF are quantitatively consistent with those from the DFT calculations previously in ref 40. This further confirms the rationality of the force field ModFF. Second, through the molecular dynamics at

experimental results, better consistency in coordination structures is achieved with our modified force field (named as ModFF hereafter) than with AMBER ff03 force field (as FF03 for short). The details are described in the following section. Solvation is the basic process for peptide systems. To correctly characterize the solvation provides a reliable background to describe the ion−peptide interactions. We solvate a copper ion in a water environment. After a 2000-step minimization and a 50 ps heating simulation, the system is relaxed with a 20 ns molecular dynamics at 300 K. As a comparison, this simulation is carried out for both ModFF and FF03, respectively. On the basis of the trajectories of these relaxation simulations, the radial distribution function (RDF) for the distance between copper ion and the oxygen atom of water molecule (Ow) is generated (as shown in Figure 2a). It

Figure 2. Radical distribution function (RDF) of copper solvation system derived from ModFF and FF03 (a) and the minimization conformation of model Cu-HGGGW by ModFF and FF03, together with the crystallographic structure (b). Only coordinating water molecules are shown in (b) for clarity.

is observed that, for both force fields, the hexacoordinate octahedral structures are achieved, as suggested from the areas of the first peak of the RDF. The coordination number and the structure are consistent with the general knowledge from experiments. Quantitatively, for these two force fields, the equilibrium distances between copper and water oxygen are different. For ModFF, the Cu−Ow distance (dCu, Ow) is 1.99 Å at the first RDF peak, and for FF03, the dCu, Ow at the first RDF peak is 1.91 Å. A larger bond length dCu, Ow for the ModFF originates from the weakening of electrostatic interaction because of the charge-transfer effect. As a reference, we

Table 1. Derived Copper−Ligand Distances from Model Systems minimization (Å) model system Cu-HIE ModFF FF03 Cu-4HIE ModFF QMa FF03 Cu-HGGGW ModFF X-rayb FF03

N1

N2

MD average (Å)

N3

N4

2.04 2.53 2.02

N1

N2

N3

N4

2.07 2.00

1.99

2.05

2.10

-

-

2.07

2.19 G3N 1.97 1.92 3.15

2.32 G3O 2.07 2.07 1.82

H1N 2.02

G2N 2.10

G3N 1.99

G3O 2.13

-

-

-

2.12

2.00−2.02 2.08 H1N 2.04 1.99 1.96

2.22 G2N 2.02 2.01 3.53

a

Data taken from ref 40. bData taken from ref 14. 5220


Article

The Journal of Physical Chemistry B room temperature, the model systems with ModFF and FF03 show different thermodynamic stability. For FF03, both CuHIE and Cu-4HIE cannot sustain stable coordination during the simulations. All the coordinations are destroyed. This is not consistent with the preference of coordination between copper ion and histidine.55−57 While for ModFF, the Cu-HIE structure is rather stable, and for the Cu-4HIE system, there are two histidines remaining bound with the copper ion during simulations. This kind of coordination with two Nδ histidines has been proposed experimentally13,15 and evaluated to be energetically favored by DFT calculations40 for the copper− octapeptide 1:2 coordination mode. Compared with FF03, the enhancement of the stability comes from the proper considerations for polarization and deprotonation. The cooperation between these interactions make the minimized bond length dCu, N in the Cu-4HIE system even smaller than that in Cu-HIE system. Solely electrostatic interactions in FF03 are not sufficient to describe the coordination bonds. As a conclusion, the differences of both geometry and stability suggest that the force field ModFF can give out reasonable stability for the coordination with copper ions. Going beyond the model systems, we also compare these two kinds of force fields for the copper−octapeptide system. The core part of the complex, that is, the coordination between the copper ion and the segment pentapeptide HGGGW, is determined with the similar procedures (both minimization and equilibration). It is observed that both force fields can achieve the coordination with the proper ligands, but a better description for the coordination distance can be obtained with ModFF force field (with relative errors smaller than 5% at maximum). This is consistent with the previous results for ModFF. Differently, the coordination is quickly destructed for regular force field during the equilibrium, while the ModFF ensures the stability of the complex during the whole simulations. As a remark, in the minimization and simulations with FF03 force field, the Cu−O interaction seems to be overestimated (with a short minimized distance and strong binding preference compared with other ligands). This is similar to the above observation about the Cu-water interaction in FF03. This demonstrates the improper balance for the interactions between copper ion and various ligands in FF03 force field. The results clearly exemplify the necessity to employ the ModFF in the studies of copper−octapeptide system. Energy Landscape of Copper−Octapeptide System. Since the coordination interaction is strong, to observe the binding/unbinding equilibrium requires long regular simulations, which is prohibited for present computational resources. Some enhanced sampling techniques are expected. In present work, REMD is employed to sample the whole energy landscape of the copper−octapeptide system. The detailed method is described in the Methods section. In Figure 3a, the PMF at a certain temperature (350 K) is given. The reaction coordinate is defined as the summation of the lengths of the coordination bonds related to the copper ion. It is observed that, after 40 ns (as the red line shown in Figure 3a), the system can approach the region with large dtot (even with dtot > 18 Å, the condition with no more than one coordinated native ligand), but the profile fluctuates largely for dtot > 18 Å, which suggests that equilibrium is not achieved. Comparing the results from various lengths of trajectories, it is observed that the part of PMF with dtot < 13 converges fast after tens of nanoseconds, and the convergence of the whole

Figure 3. Reweighted PMF profiles from different simulation time lengths (a) or temperatures (b). The sum of distances between native coordinating atoms and the copper ion (dtot) is used as CV. Profiles in (a) are derived at temperature 350 K. The first 10 ns of sampling is dropped in PMF calculation.

PMF is accomplished after 110 ns. This is because the part with larger dtot has a larger conformational space, and more conformational heterogeneity may be involved. On the basis of these kinds of variations, it can be postulated that the transition state is around the region between dtot = 13 and 14 Å, which separates two different areas of conformational space. On the basis of the trajectories of REMD, the thermodynamic features of the coordination processes can be characterized. Note that, in our following thermodynamic analysis, the first 10 ns trajectory is not included to remove the effects of the initial conformations. The remainder of the simulations is about 104.3 ns at each temperature. The thermodynamic features are determined on the basis of the reweighting method. Among various features, the free-energy landscapes under different temperatures are fundamental to understand this copper−octapeptide system. It is observed that all the coordinations are rather stable at 340 K and the landscape looks like a funnel with the native structure as the minimum. There are few dissociated (unbounded) conformations observed. Following the increase of temperature, the native conformation (or basin) still acts as the global minimum when the temperature is lower than 400 K. This observation indicates that a fully coordinated state still takes a dominant role for the thermodynamics of the system, though some intermediates emerge. However, at the temperature higher than 400 K, the landscape apparently biases toward the regions with large dtot, with the global minimum gradually increasing from 8 to19 Å, implying a rough two-state transition from the native conformation to the fully unbounded state. Meanwhile, the variation of the coordination number of the copper− octapeptide system has a similar trend. Clearly, the dissociation is not an all-or-nothing process, which indicates that the dissociation of various coordination bonds are not cooperative. Multiple states are necessary to describe the processes. To have a detailed understanding on the conformation variations of the copper−octapeptide system, the lengths of various coordination bonds related to copper ion (namely the lengths related to bonds Cu−H1N, Cu-G2N, Cu-G3N, and Cu-G3O, respectively) are used as the coordinates to characterize the free-energy landscape. That is, a fourdimensional landscape can be generated. Since there are few other nonbonded interaction besides the coordination in our system, we focus on the local minima on the landscape spanned by the lengths of coordination bonds. With the density-based projective clustering method EPCH,58,59 the conformations are clustered, and various local minima on the landscape are identified. Practically, there are a total of six two5221


Article

The Journal of Physical Chemistry B dimensional projected spaces, each of which spans the lengths of two bonds. We first collect the conformations around the local minima in a certain projected space. Then, we check if these conformations are located at the local minima in all other projected spaces. Only the conformations which are around the local minima in any projected space are collected. The conformations which are related to the same set of local minima in projected space share the same coordination properties and are classified into the same cluster. By scanning all the snapshot conformations in one or more trajectories, a series of conformation clusters are obtained. For the copper− octapeptide system, to visualize various coordination modes during its dynamics, three trajectories at the temperatures 340.8 K, 352.2 K, and 358.1 K are employed, under the consideration that the system can fluctuate largely in the whole conformational space and various coordination modes can be observed. As a result, the largest 15 clusters are picked out, and the typical structures of these clusters are shown in Figure 4c.

Table 2. Coordination of the 15 Largest Clusters cluster

H1N

G2N

G3N

G3O

4A 3A 3B 2C 2A 2B 2C 1A, 1B 1C 1D 1E, 1F 0A 0B

• • •

• • • • • •

• • • • • • •

•

H1Oa

water bridgeb •

• • •

• •

• • •

•

• • •c •

a

Nonnative ligand. bThe interaction between copper and the ligand in which the ligand is bridged to a copper-coordinating water by a hydrogen bond. Bridged to G3O if not specially noted. cBridged to G3O and G3N simultaneously.

these clusters are largely populated on the landscape especially at the temperature 352.2 K. It is found that not all the combinations of coordination bonds are observed in these clusters. For the case NC = 3, only two kinds of combinations of coordination bonds are observed: one with the bonds Cu-G2N, Cu-G3N, and Cu− H1N (for clusters 3A and 3B); and the other with the bonds Cu-G2N, Cu-G3N, and Cu-G3O (for cluster 3C). It is also found that the populations related to these combinations vary largely. The populations for cluster 3B is about 10 times those for cluster 3C. This reflects that the strength of the coordination bond Cu−H1N is apparently stronger than that of the bond Cu-G3O. This may be attributed to the distortion of the coordination geometry. For the first case without the bond Cu-G3O, the local arrangement of ligand atoms are not around their native positions. The bond Cu−H1N would like to be perpendicular to the Cu-G2N-G3N plane. This kind of structure can maximize the flexibility of the imidazole group and the solvation of the copper ion. Differently, for the case with Cu−H1N dissociated, the coordination geometry is similar to the native conformation. The comparison of these conformations implies that there is intrinsic frustration in the native coordination structure. This demonstrates that the frustration widely exists in nature.60−62 On the other hand, the common bonds (Cu-G2N and Cu-G3N) in these clusters contribute importantly to the stability of the native coordination structure to overcome the frustration. This reflects an essential contribution of G2N-G3N coordination to the overall stability, as a “core” of the coordination structure. For the clusters with NC = 1 or 2, the nonuniform occurrences of the coordination bonds is also observed. The bond Cu-G3N is the most populated in clusters with NC = 2 and Cu−H1N in clusters with NC = 1. Interestingly, a nonnative coordination bond Cu−H1O is frequently observed, as in clusters 1A, 1B, 1C, and 2C. This bond contributes essentially to the stability of the structures related to the concerned clusters, though this bond is not counted in the number NC in the naming procedure of clusters. These observations suggest that various coordination bonds may contribute differently to the stability and dynamics of the copper−octapeptide system.

Figure 4. Energy landscape and representative configurations of clusters. (a): the free energy landscape at temperature 350 K. (b): population of each cluster at temperatures 340.8, 352.2, and 358.2 K. (c): representative configurations of each cluster. The copper ion is colored from blue (hydrated) to orange (dehydrated) according to the hydration status.

Considering the generation of these clusters, these clusters can be categorized on the basis of their coordination number NC. For example, the conformations in the cluster 4A all have four native coordination bonds. Clearly, there is only one cluster corresponding to NC = 4 because there is only one possible combination for four native coordination bonds, which is consistent with the uniqueness of the native coordination structure. Indeed, this cluster is the largest one because the native conformation is the most stable one comparing to other clusters. In more detail, all the concerned coordination bonds related to these 15 clusters are recorded in Table 2. In all these clusters, the copper is penta- or hexa-coordinated by amino acid ligands and water molecules. The coordinations with water molecules are not explicitly declared in the table. The corresponding populations at various temperatures of these clusters are also given in Figure 4b. This result indicates that all 5222


Article

The Journal of Physical Chemistry B To visualize the relations between these structure clusters, these clusters are assigned on a two-dimensional landscape spanning with the distances dCu−H1N + dCu−G3O and dCu−G2N + dCu−G3N at the temperature of 350 K. The locations of all the clusters are shown in Figure 4a. Almost all the minima are identified, except for two minima located at the left part of the panel, because the populations of the clusters related to these two minima is rather small. Locally, some clusters (1D, 1E, and 1F) are degenerate on this projection of the landscape. These minima (clusters) are separated with a rough gap along dCu−G2N + dCu−G3N ≈ 7 Å, which divides the landscape into two regions. When the distance dCu−G2N + dCu−G3N < 7 Å (>7 Å), all the clusters have their NC ≥ 2 (≤1). In each region, the barriers among the minima (clusters) are relatively small compared with the barrier to separate two regions. With this type of feature, a rough two-state behavior of the system would be expected as observed in Figure 3b. It is clear that this kind of transition is accompanied by the variation of coordination number, which indicates that the coordination number is a reasonable order parameter. This is consistent with the physics that the coordination interaction is the strongest one during the structure formation. Different from regular two-state folding, there are multiple minima in each basin. This reflects the competition between other interactions and coordination. As a result, there are many relatively weak transitions following the variation of temperature, both before and after the main two-state transition, which deteriorates the two-state behavior. This illustrates the microscopic reason for the complex thermodynamic behaviors. Since the barrier is not uniform along the boundary dCu−G2N + dCu−G3N = 7 Å, the transition between two regions can happen only in the condition with specific dCu−H1N + dCu−G3O. The transition would occur with certain cooperativity between various coordination bonds. This suggests that there are some specific pathways which imply that various coordination bonds contribute differently to the binding processes. Pathway of Coordination. The above landscape suggests that the formation of various coordination bonds is not random. A detailed characterization on the binding order of the coordination bonds may help to understand the binding of the ion to peptide. First, we make statistics on the occurrence of various coordination bonds following the increase of the coordination number. As shown in Figure 5, the occurrence of each native coordination bond (marked with the ligand atom) is measured as the probability it occurred in the relevant snapshot among all snapshots with a same Nc. Here, a threshold length of copper-ligand bond is assigned as 2.8 Å. It is observed that the probabilities related to G2N and G3N increase monotonically during the progress of the coordination. These two bonds vary concurrently and take the role of the stabilization core of the native coordination structure. Without these interactions, the native coordination pattern cannot be established. We carried out a series of unbinding simulations (a total of 260 trajectories at 350 K). The unbinding trajectories often obstructed dominantly at the clusters 2A/2B. This indicates the role of coordinations with G2N/G3N for stability. Differently, the probabilities related to H1N and G3O vary nonmonotonically. This implies that there are frustrations related to these two bonds which may experience unbinding−rebinding processes or competitions of multiple pathways. Especially, the bond related to H1N plays an important role in the early stage of the coordination.

Figure 5. Probability of interactions (coordination and water bridge) between copper and each native ligand atom and coordination probability of nonnative ligand H1O. An illustration of water bridge is shown in the cell figure on top left corner of H1N plot.

This suggests that this bond may help for the recognition between ion and the unfolded peptide. This is partially because of the large side-chain of the concerned histidine. The flexibility related to the side chain may facilitate the fly casting63 processes and thus the recognition. On the other hand, the bond related to G3O, a water-mediated interaction between ion and the atom G3O (namely water bridge) is observed (as shown by the high green bar related to G3O with NC = 0). This kind of hydrogen-bond bridge gives a large distance and a weak affinity between ion and G3O atom. Note that the interactions between ion and H1N/G3O are related to the side-chain/water. To achieve the native conformation of the backbone, the deformations related to the side chain or the breaking of water bridge is necessary. This may introduce some frustrations and cause the variations of the binding percentage of these bonds. This information helps us to establish a picture about the coordination: the bonds related to G2N and G3N act as a “core” responsible for the stability, while the bonds related to H1N and G3O act as “arm” related to recognition. Some of the coordination bonds are positively or negatively correlated, which can be observed from the unbinding simulations. For the cases with two coordinations with G2N and G3N, the dissociations of these two ligands are almost concurrent. The mean lag time between dissociation of G2N and G3N is 13.74 ps. This clearly demonstrates the cooperative dynamics of the coordinations with G2N and G3N. When the coordination with G3O exists, the mean lag time between dissociation of G2N and G3N increases apparently, that is, 214.5 ps. Besides, the dissociation of G2N is always ahead that of G3N. This reveals a positive correlation between G3O and G3N and a negative correlation between G3O and G2N. When the coordination of H1O is fixed, that same order of dissociations of G2N and G3N is observed, but the lag time between the dissociations of G2N and G3N reaches 2647 ps. This implies that the effect of H1O is similar to (but stronger than) that of G3O. Especially, the 5223


Article

The Journal of Physical Chemistry B large lag time indicates a strong bias toward different ligands, which suggests that H1O may take an important role to modulate the binding processes in copper−octapeptide interaction. These events exemplify the competitive interactions between coordinations with various ligands. The occurrence and correlation of the coordination bonds exemplify the complexity of the coordination processes. There are different behaviors for “core” and “arm” parts, and nonmonotonic variations are observed for the ligands H1N and G3O. The coordinations of ligands are not independent, and these cannot be explained with a simple progressive model. To capture the underlying physics behind these observations, the detailed pathways are analyzed on the basis of the trajectories. Along the trajectories, the transitions between conformational clusters are recorded on the basis of the REMD trajectories. These transitions build up a transition network, which describes the possible evolution of conformations during binding/unbinding. This idea is similar to the Markov Model64 to rebuild the landscape from local dynamics (transitions). The difference from the Markov Model is that the cluster of conformations are presumed on the basis of the information of coordination bonds. The transition network is shown in Figure 6a. Several small clusters (not appearing in Figure 4) are included because they take irreplaceable roles in the transition network. The coordination information of all clusters can be found in Supporting Information. In the transition network, the counts of transitions between clusters are marked along with the transition direction. It is observed that the transitions between any two clusters during the simulations are roughly equal. This result suggests that the detailed balance is roughly satisfied and implies the convergence of our simulations. On the transition network, the conformation clusters are further merged together on the basis of the horizontal transitions (namely the transitions between the clusters with same NC). It is found that there are roughly several categories of pathways (as shown in Figure 6a). All these dominant pathways are projected onto the PMF landscape (as shown in Figure 6b). At the beginning of the binding, there are three categories of pathways (as marked blocks with different colors), and these pathways merged gradually following the progress of the binding processes. Among these pathways, the red and green pathways are initiated with the coordinations with “arm” residues, and the blue pathway starts from the binding of one “core” residue. Based on the transition counts, the “arm” residues are likely to be bound during the initial kinetic partition. This can be understood on the basis of the fly casting mechanism. The “arm” residues are generally flexible. Especially for the histidine, the coordination site H1N is on its side chain. The rotation or swing of the side chain does not affect the conformation of the backbone of the peptide chain and introduces more flexibility for the binding site. Besides, the interaction wth the site G3O is extended through an intermediate water molecule. This water-bridged interaction is also observed in the clusters with NC = 0. During the binding with the copper ion, these “arm” residues can be elongated to increase their capture radii and produce a higher probability to seize the ion. Differently, the “core” residues are inside the peptide chain. The related binding sites are rather rigid. The recognition of these residues by copper ions would be difficult because of the small capture radii and the screening effect by the “arm” residues. If we further reduce the distance constraint between the copper ion and the peptide, a larger difference for these pathways may be observed. Following the progress of the

Figure 6. Transition statistics between major clusters (a). Only transitions counted over 10 in at least one direction are shown, except those between 2D and 3D for the connectivity of network. For clusters at the same NC, bidirectional arrows are used where both bidirectional counts exceed 100. Zero-coordinated clusters 0A, 0B, clusters 1A, 1B, 1C, 1E, 1F, 1G, 1I corresponding to H1N coordination, and clusters 1J, 1H corresponding to G3N coordination are grouped, i.e. the in and out transitions are counted for the whole group. Clusters that have a water bridge between copper and G3O are marked by double box frame, and those that have H1O coordinated are marked by * on the cluster name. Putative major binding pathways are projected to the PMF landscape (b) with several representative transient conformations (c). Only representative clusters are labeled in (b) for clarity.

binding process, there are some interpathway transitions around some clusters (such as the clusters 2A and 2C). Especially, the pathways starting from H1N (green) and from G3O (red) merge into a single pathway through the cluster 3D. Finally, all pathways can reach the native conformations. This outlines a competitive picture among various pathways. In detail, these pathways are different not only in their starting positions. For the pathways staring from “core” residues (blue), the transition counts vary weakly for various NCs. A barrier comes from the initial binding of the ion. This indicates that the binding along this pathway is in a progressive manner. The large transition between the clusters NC = 1 (1J and 1H) and NC = 2 (2B) suggests that there is a positive cooperativity between two “core” residues (G2N and G3N), that is, the formation of one “core” residue would help the binding of the other residue. On the other hand, for the pathways starting from “arm” residues, there are bottlenecks between the clusters 2C/2D and 3D. The transitions are rare between these clusters compared with other transitions along the pathways. This indicates that the transitions between these clusters are the rate-limiting steps, which are related to the binding of the ligands G3N or H1N in the cluster 3D. The 5224


Article


nonnative coordination between copper ion and the ligand H1O, which is widely observed in the clusters with NC ≤ 2. This interaction can also be realized through a water bridge, which increases the range of interaction. As an example, some small clusters at the early stage of the binding are extracted (namely T1 to T5 in Figure 6c). The interaction with the ligand H1O (or through a water bridge) exists widely. At the late stage, a nonnative-to-native switch (i.e., the coordination switching from H1O to G2N) is observed. This kind of phenomenon is previously declared in many metalloprotein systems, such as azurin,66 cytochrome c,66 zinc finger,35 calmodulin,65 and amyloid-β.67 Our simulation also observes such kind of microscopic behavior. Moreover, in our simulations, it is suggested that the binding with the nonnative ligand is obligatory for the binding processes. Clearly, the nonnative interaction is usually accompanied with local frustration.66 How does the frustration help the binding? This is an interesting question to understand the ion-binding processes. Microscopically, the nonnative ligand increases the local density of ligands to coordinate with the copper ion, which makes the ions be attracted around the coordination center and eases the further recognition. Additionally, the binding of nonnative ligands H1O can help to lower the freeenergy barrier related to ion dehydration.65 This view was also observed in the previous study on the DFT energy evaluations on successive deprotonation of Cu-HGGGW complex,68 which pointed out the favored H1O participation at one deprotonation stage. This quantitatively explains the large occurrence of the coordination with H1O at the early stage of binding and supports our picture of binding. Meanwhile, together with several other native coordination bonds formed, the chargetransfer effect reduces the energetic penalty to break the concerned nonnative coordination. The switching happens easily in consequence. The many-body features of the coordination interaction smooth the frustration introduced by the nonnative coordination. The small energetic barrier can be inferred from the large flux during the concerned transitions (such as from 1H/1J to 2B, and from 3D to 3A/4A). Furthermore, the existence of this nonnative interaction modulates the role of the “core” ligands. Because of the competition of the ligand H1O, the binding of the copper ion with the “core” ligands is generally initiated from the ligand G3N. This explains the asymmetry of two “core” residues. On the other hand, the interaction between copper ion and the ligand H1O partly restricts the directions of the backbone. It makes the concurrent bindings with the ligands H1N and G3N become difficult. This hidden degree of freedom strengthens the anticooperativity between the ligands H1N and G3N. Indeed, the nonnative coordination illuminates the whole picture of ion-binding processes.

difficulty to form these two coordination bonds concurrently implies that the coordinations related to G3N and to H1N are exclusive to each other. With this kind of kinetic bottleneck, the binding behaves more like a two-state process. Clearly, different kinetic behaviors may happen in the simple copper− peptide system. There are some effective interactions between various ligands. Physically, considering the microscopic interactions in the peptide system, this kind of effective interaction is originated from the entropic penalty between residues. It is easy to find out that the neighboring residues are often positively correlated, such as the cooperation between G2N and G3N. The ligands in these residues are generally spatially close to each other. Differently, for the distant residues which are not neighboring, the conformational entropy reduces the preference to bind with the same copper ion. Through the simulations on the peptide without the copper ion, the probabilities to have the distance smaller than the average distance in the coordinated state for various pairs of ligands are calculated. It is observed that the probability for neighboring ligands is generally at least 3 times that for distant ligands. Even in such a small system, the conformational properties are still important to affect the dynamics of the binding with ions. The existence of multiple pathways can explain the complexity observed in the evolution of coordination bonds. First, the biased kinetic partition toward “arm” and “core” residues introduces the differences between these kinds of ligands. The preference toward the ligands H1N and G3O due to the fly casting mechanism makes them more popular at the beginning of the binding process. Thus, a larger occurrence for H1N (or G3O) is observed for the case NC = 1. In this stage, the capture radius dominates the recognition of the coordination. Similarly, some ligands interact with copper ions through a bridge of water molecules. This increases the range of interaction and thus the capture radius. In consequence, the interactions through water bridges are favored. It is the reason there are more water-bridged interactions for small NC. Previously, it has been determined that water-mediated coordination contributes significantly to lowering the free-energy barrier of ion dehydration in the calmodulin system.65 This result is well in line with our observations. Second, the nonmonotonic trends for the coordination bonds related to “arm” residues is originated from the superposition of different pathways. Clearly, for each pathway, the coordination bonds are formed one by one. There is no breakage of any coordination bonds along the binding pathways. The nonmonotonic behaviors are antiintuition at first sight. Because of the existence of multiple pathways with different properties, the occurrences of various coordination bonds are the superposed results of all possible pathways. For the pathways starting from “arm” residues, the bottlenecks during the pathways make some clusters (mainly around NC = 2) rather rare. While for the pathway starting from “core” residues, the cooperativity apparently enriches the content of certain clusters (such as the clusters 2B and 2A). The superposition of these pathways produce a large occurrence for “core” residues and suppresses the manifestation of the “arm” residues when NC = 2. This results in the decrease of the occurrence of “arm” residues and gives out a nonmonotonic variation for the concerned coordination bonds. Among these pathways, there is another interaction which contributes essentially to the binding processes. This is the

■

CONCLUSIONS To simulate the ion-binding dynamics, we develop a model which enables classical molecular dynamics to deal with the interaction between protein systems and the copper ion by including charge transfer, protonation/deprotonation, and polarization effects into force field. This model is verified by a series of tests and is applied to the case of Cu-PHGGGWGQ complex. Through simulations, the free-energy landscape and binding pathways are determined. Various residues show their different roles in stability and ion arresting. The contributions of a nonnative ligand H1O are discussed. These diverse mechanisms in such a small peptide system can be observed in 5225


The Journal of Physical Chemistry B various protein systems. Their combinations and the underlying physical processes in copper−octapeptide system are instructive for us to understand PrP misfolding and other complex protein dynamics. For example, for the multiplehistidine binding mode of octapeptides, the critical concentration is apparently lower than that for the case we investigated here, though the multiple-histidine mode is less stable compared with the full coordination to monomeric octarepeats.40 This can be ascribed to the large capture radius related to the flexibility of the histidine. The large occurrence of histidine in metalloproteins may be related to this feature and might be a result of natural evolution. The large capture radius of histidines may also contribute to the intra- and intermolecular cross-links (such as His96−Cu2+−His111 interaction21,22). Through these interactions between copper ions and octapeptide segments, the flexibility of the N-terminal part may be reduced. This kind of effect can reduce the conformational entropy to make the formation of the pincerlike structure with RNA69 and other factors easier and thus speed up the consequent misfolding of PrP. Besides, the copper ions and RNA generally share some binding sites in PrP (for instance the basic amino acid histidine). The coordinations between copper ions and histidines can block some unnecessary binding sites. For example, the His96−Cu2+− His111 cross-link may possibly impede the binding between RNA and polybasic region 111−121, which may enhance the effect of RNA on PrP misfolding. Additionally, the intermolecular cross-linking can effectively increase the local concentration of PrP molecules, which would introduce more interactions between PrP molecules and facilitate the growth of aggregations.19 We believe that a comprehensive understanding on the cooperation and competition between copper ions, RNA, and various other factors would be necessary to uncover the mystery behind the misfolding and aggregations of PrP. It is worth pointing out that our model is simple and precise enough and can be applied further to other metalloproteins. Our force field, analyzing methods, and resultant pictures would be a starting point for more complicated coordination processes.

■

■

ACKNOWLEDGMENTS

■

REFERENCES

This work was supported by the National Natural Science Foundation of China (No. 11774157, and 11334004). The authors acknowledge HPCC of Nanjing University and of NLM for the computational support. We also thank Prof. Glenn L. Millhauser for the data of Cu-HGGGW complex.

(1) Festa, R. A.; Thiele, D. J. Copper: An Essential Metal in Biology. Curr. Biol. 2011, 21, R877−R883. (2) Rae, T. D.; Schmidt, P. J.; Pufahl, R. A.; Culotta, V. C.; O’Halloran, T. V. Undetectable Intracellular Free Copper: The Requirement of a Copper Chaperone for Superoxide Dismutase. Science 1999, 284, 805−808. (3) Banci, L.; Bertini, I.; Ciofi-Baffoni, S.; Kozyreva, T.; Zovo, K.; Palumaa, P. Affinity Gradients Drive Copper to Cellular Destinations. Nature 2010, 465, 645−648. (4) Kaplan, J. H.; Maryon, E. B. How Mammalian Cells Acquire Copper: An Essential but Potentially Toxic Metal. Biophys. J. 2016, 110, 7−13. (5) Millhauser, G. L. Copper and the Prion Protein: Methods, Structures, Function, and Disease. Annu. Rev. Phys. Chem. 2007, 58, 299−320. (6) Quintanar, L.; Rivillas-Acevedo, L.; Grande-Aztatzi, R.; GomezCastro, C. Z.; Arcos-Lopez, T.; Vela, A. Copper Coordination to the Prion Protein: Insights from Theoretical Studies. Coord. Chem. Rev. 2013, 257, 429−444. (7) Aronoff-Spencer, E.; Burns, C. S.; Avdievich, N. I.; Gerfen, G. J.; Peisach, J.; Antholine, W. E.; Ball, H. L.; Cohen, F. E.; Prusiner, S. B.; Millhauser, G. L. Identification of the Cu2+ Binding Sites in the NTerminal Domain of the Prion Protein by EPR and CD Spectroscopy. Biochemistry 2000, 39, 13760−13771. (8) Burns, C. S.; Aronoff-Spencer, E.; Legname, G.; Prusiner, S. B.; Antholine, W. E.; Gerfen, G. J.; Peisach, J.; Millhauser, G. L. Copper Coordination in the Full-Length, Recombinant Prion Protein. Biochemistry 2003, 42, 6794−6803. (9) Jackson, G. S.; Murray, I.; Hosszu, L. L. P.; Gibbs, N.; Waltho, J. P.; Clarke, A. R.; Collinge, J. Location and Properties of MetalBinding Sites on the Human Prion Protein. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 8531−8535. (10) Qin, K. F.; Yang, Y.; Mastrangelo, P.; Westaway, D. Mapping Cu(II) Binding Sites in Prion Proteins by Diethyl Pyrocarbonate Modification and MatrixAssisted Laser Desorption IonizationTime of Flight (MALDI-TOF) Mass Spectrometric Footprinting. J. Biol. Chem. 2002, 277, 1981−1990. (11) Klewpatinond, M.; Davies, P.; Bowen, S.; Brown, D. R.; Viles, J. H. Deconvoluting the Cu2+ Binding Modes of Full-Length Prion Protein. J. Biol. Chem. 2008, 283, 1870−1881. (12) Viles, J. H.; Klewpatinond, M.; Nadal, R. C. Copper and the Structural Biology of the Prion Protein. Biochem. Soc. Trans. 2008, 36, 1288−1292. (13) Chattopadhyay, M.; Walter, E. D.; Newell, D. J.; Jackson, P. J.; AronoffSpencer, E.; Peisach, J.; Gerfen, G. J.; Bennett, B.; Antholine, W. E.; Millhauser, G. L. The Octarepeat Domain of the Prion Protein Binds Cu(II) with Three Distinct Coordination Modes at pH 7.4. J. Am. Chem. Soc. 2005, 127, 12647−12656. (14) Burns, C. S.; Aronoff-Spencer, E.; Dunham, C. M.; Lario, P.; Avdievich, N. I.; Antholine, W. E.; Olmstead, M. M.; Vrielink, A.; Gerfen, G. J.; Peisach, J.; et al. Molecular Features of the Copper Binding Sites in the Octarepeat Domain of the Prion Protein. Biochemistry 2002, 41, 3991−4001. (15) Wells, M. A.; Jelinska, C.; Hosszu, L. L. P.; Craven, C. J.; Clarke, A. R.; Collinge, J.; Waltho, J. P.; Jackson, G. S. Multiple Forms of Copper (II) Coordination Occur Throughout the Disordered NTerminal Region of the PrionProtein at pH 7.4. Biochem. J. 2006, 400, 501−510.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b02457.

■

Article

Validity of simulation procedures, data of structure of our model system (Cu-HGGGW complex) in this work and the parameters derived for the force field, and the coordination bonds related to the clusters in transition network (PDF)

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. ORCID

Ke Chen: 0000-0003-0433-5580 Wenfei Li: 0000-0003-2679-4075 Jun Wang: 0000-0003-2954-2607 Wei Wang: 0000-0001-5441-0302 Notes

The authors declare no competing financial interest. 5226


Article


(36) Sakharov, D. V.; Lim, C. Zn Protein Simulations Including Charge Transfer and Local Polarization Effects. J. Am. Chem. Soc. 2005, 127, 4921−4929. (37) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A PointCharge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (38) Riihimäki, E.-S.; Kloo, L. Computational Comparison of Cation Coordination to Human Prion Peptide Models. Inorg. Chem. 2006, 45, 8509−8516. (39) Babu, C. S.; Lim, C. Empirical Force Fields for Biologically Active Divalent Metal Cations in Water. J. Phys. Chem. A 2006, 110, 691−699. (40) Hodak, M.; Chisnell, R.; Lu, W. C.; Bernholc, J. Functional Implications of Multistage Copper Binding to the Prion Protein. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 11576−11581. (41) Loewenthal, R.; Sancho, J.; Fersht, A. R. Histidine-Aromatic Interactions in Barnase. Elevation of Histidine pKa and Contribution to Protein Stability. J. Mol. Biol. 1992, 224, 759−70. (42) Fernandez-Recio, J.; Vazquez, A.; Civera, C.; Sevilla, P.; Sancho, J. The Tryptophan/Histidine Interaction in Alpha-Helices. J. Mol. Biol. 1997, 267, 184−97. (43) Bonomo, R. P.; Cucinotta, V.; Giuffrida, A.; Impellizzeri, G.; Magri, A.; Pappalardo, G.; Rizzarelli, E.; Santoro, A. M.; Tabbi, G.; Vagliasindi, L. I. A Reinvestigation of Copper Coordination in the Octa-repeats Region of the Prion Protein. Dalton Trans. 2005, 0, 150−158. (44) Miura, T.; Hori-i, A.; Mototani, H.; Takeuchi, H. Raman Spectroscopic Study on the Copper(II) Binding Mode of Prion Octapeptide and Its pH Dependence. Biochemistry 1999, 38, 11560− 11569. (45) Nadler, W.; Hansmann, U. H. E. Optimized Explicit-Solvent Replica Exchange Molecular Dynamics from Scratch. J. Phys. Chem. B 2008, 112, 10386−10387. (46) Nadler, W.; Hansmann, U. H. E. Dynamics and Optimal Number of Replicas in Parallel Tempering Simulations. Phys. Rev. E 2007, 76, No. 065701. (47) Nadler, W.; Hansmann, U. H. E. Optimizing Replica Exchange Moves for Molecular Dynamics. Phys. Rev. E 2007, 76, No. 057102. (48) Sindhikara, D. J. Modular Reweighting Software for Statistical Mechanical Analysis of Biased Equilibrium Data. Comput. Phys. Commun. 2011, 182, 2227−2231. (49) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. J. Comput. Chem. 1992, 13, 1011−1021. (50) The PyMOL Molecular Graphics System, version 1.8; Schrödinger, LLC: Cambridge, MA, 2015. (51) Schwenk, C. F.; Rode, B. M. Extended Ab Initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations of Hydrated Cu2+. J. Chem. Phys. 2003, 119, 9523−9531. (52) Schwenk, C. F.; Rode, B. M. New Insights into the Jahn-Teller Effect through Ab Initio Quantum-Mechanical/MolecularMechanical Molecular Dynamics Simulations of CuII in Water. ChemPhysChem 2003, 4, 931−943. (53) Ohtaki, H.; Radnai, T. Structure and Dynamics of Hydrated Ions. Chem. Rev. 1993, 93, 1157−1204. (54) Persson, I.; Persson, P.; Sandström, M.; Ullström, A.-S. Structure of Jahn-Teller Distorted Solvated Copper(II) Ions in Solution, and in Solids with Apparently Regular Octahedral Coordination Geometry. J. Chem. Soc., Dalton Trans. 2002, 0, 1256−1265. (55) Giachin, G.; Mai, P. T.; Tran, T. H.; Salzano, G.; Benetti, F.; Migliorati, V.; Arcovito, A.; Longa, S. D.; Mancini, G.; D’Angelo, P.; Legname, G. The Non-octarepeat Copper Binding Site of the Prion Protein Is a Key Regulator of Prion Conversion. Sci. Rep. 2015, 5, 15253.

(16) Yen, C.-F.; Harischandra, D. S.; Kanthasamy, A.; Sivasankar, S. Copper-Induced Structural Conversion Templates Prion Protein Oligomerization and Neurotoxicity. Sci. Adv. 2016, 2, No. e1600014. (17) Younan, N. D.; Klewpatinond, M.; Davies, P.; Ruban, A. V.; Brown, D. R.; Viles, J. H. Copper(II)-Induced Secondary Structure Changes and Reduced Folding Stability of the Prion Protein. J. Mol. Biol. 2011, 410, 369−382. (18) Leal, S. S.; Botelho, H. M.; Gomes, C. M. Metal Ions as Modulators of Protein Conformation and Misfolding in Neurodegeneration. Coord. Chem. Rev. 2012, 256, 2253−2270. (19) Bocharova, O. V.; Breydo, L.; Salnikov, V. V.; Baskakov, I. V. Copper(II) Inhibits in Vitro Conversion of Prion Protein into Amyloid Fibrils. Biochemistry 2005, 44, 6776−6787. (20) Quaglio, E.; Chiesa, R.; Harris, D. A. Copper Converts the Cellular Prion Protein into a Protease-Resistant Species That Is Distinct from the Scrapie Isoform. J. Biol. Chem. 2001, 276, 11432− 11438. (21) Pushie, M. J.; Rauk, A.; Jirik, F. R.; Vogel, H. J. Can Copper Binding to the Prion Protein Generate a Misfolded Form of the Protein? BioMetals 2009, 22, 159−175. (22) Pushie, M. J.; Vogel, H. J. A Potential Mechanism for Cu2+ Reduction, β-Cleavage, and β-Sheet Initiation Within the N-Terminal Domain of the Prion Protein: Insights from Density Functional Theory and Molecular Dynamics Calculations. J. Toxicol. Environ. Health, Part A 2009, 72, 1040−1059. (23) Brown, L. R.; Harris, D. A. Copper and Zinc Cause Delivery of the Prion Protein from the Plasma Membrane to a Subset of Early Endosomes and the Golgi. J. Neurochem. 2003, 87, 353−363. (24) Pauly, P. C.; Harris, D. A. Copper Stimulates Endocytosis of the Prion Protein. J. Biol. Chem. 1998, 273, 33107−33110. (25) Wu, B.; McDonald, A. J.; Markham, K.; Rich, C. B.; McHugh, K. P.; Tatzelt, J.; Colby, D. W.; Millhauser, G. L.; Harris, D. A. The NTerminus of the Prion Protein Is a Toxic Effector Regulated by the CTerminus. eLife 2017, 6, No. e23473. (26) Evans, E. G. B.; Millhauser, G. L. Copper and Zinc-Promoted Interdomain Structure in the Prion Protein: A Mechanism for Autoinhibition of the Neurotoxic N-Terminus. Prog. Mol. Biol. Transl. Sci. 2017, 150, 35−56. (27) Hansmann, U. H. E. Parallel Tempering Algorithm for Conformational Studies of Biological Molecules. Chem. Phys. Lett. 1997, 281, 140−150. (28) Sugita, Y.; Okamoto, Y. Replica Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141−151. (29) del Pino, P.; Weiss, A.; Bertsch, U.; Renner, C.; Mentler, M.; Grantner, K.; Fiorino, F.; Meyer-Klaucke, W.; Moroder, L.; Kretzschmar, H. A.; et al. The Configuration of the Cu2+ Binding Region in Full Length Human Prion Protein. Eur. Biophys. J. 2007, 36, 239−252. (30) Weiss, A.; del Pino, P.; Bertsch, U.; Renner, C.; Mentler, M.; Grantner, K.; Moroder, L.; Kretzschmar, H. A.; Parak, F. G. The Configuration of the Cu2+ Binding Region in Full-Length Human Prion Protein Compared with the Isolated Octapeptide. Vet. Microbiol. 2007, 123, 358−366. (31) Li, P.; Song, L. F.; Merz, K. M. Parameterization of Highly Charged Metal Ions Using the 12-6-4 LJ-Type Nonbonded Model in Explicit Water. J. Phys. Chem. B 2015, 119, 883−895. (32) Rick, S. W.; Stuart, S. J. Potentials and Algorithms for Incorporating Polarizability in Computer Simulations. Rev. Comput. Chem. 2002, 18, 89−146. (33) Li, P.; Merz, K. M. Metal Ion Modeling Using Classical Mechanics. Chem. Rev. 2017, 117, 1564−1686. (34) Li, W.; Wang, J.; Zhang, J.; Wang, W. Molecular Simulations of Metal-Coupled Protein Folding. Curr. Opin. Struct. Biol. 2015, 30, 25−31. (35) Li, W. F.; Zhang, J.; Wang, J.; Wang, W. Metal-Coupled Folding of Cys2His2 Zinc-Finger. J. Am. Chem. Soc. 2008, 130, 892− 900. 5227


Article

The Journal of Physical Chemistry B (56) Viles, J. H. Metal Ions and Amyloid Fiber Formation in Neurodegenerative Diseases. Copper, Zinc and Iron in Alzheimer’s, Parkinson’s and Prion Diseases. Coord. Chem. Rev. 2012, 256, 2271− 2284. (57) Telpoukhovskaia, M. A.; Orvig, C. Werner Coordination Chemistry and Neurodegeneration. Chem. Soc. Rev. 2013, 42, 1836− 1846. (58) Ng, E. K. K.; Fu, A. W.-c.; Wong, R. C.-W. Projective Clustering by Histograms. IEEE Trans. Knowl. Data Eng. 2005, 17, 369−383. (59) Procopiuc, C. M. In Encyclopedia of Machine Learning; Sammut, C., Webb, G. I., Eds.; Springer US: Boston, MA, 2010; pp 806−811. (60) Tripathi, S.; Wang, Q.; Zhang, P.; Hoffman, L.; Waxham, M. N.; Cheung, M. S. Conformational Frustration in Calmodulin-Target Recognition. J. Mol. Recognit. 2015, 28, 74−86. (61) Ren, W.; Li, W.; Wang, J.; Zhang, J.; Wang, W. Consequences of Energetic Frustration on the Ligand-Coupled Folding/Dimerization Dynamics of Allosteric Protein S100A12. J. Phys. Chem. B 2017, 121, 9799−9806. (62) Li, W.; Wolynes, P. G.; Takada, S. Frustration, Specific Sequence Dependence, and Nonlinearity in Large-Amplitude Fluctuations of Allosteric Proteins. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 3504−3509. (63) Shoemaker, B. A.; Portman, J. J.; Wolynes, P. G. Speeding Molecular Recognition by Using the Folding Funnel: The Fly-Casting Mechanism. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 8868−8873. (64) Chodera, J. D.; Singhal, N.; Pande, V. S.; Dill, K. A.; Swope, W. C. Automatic Discovery of Metastable States for the Construction of Markov Models of Macromolecular Conformational Dynamics. J. Chem. Phys. 2007, 126, 155101. (65) Tan, C. Metal Cation Coupled Multiscale Conformational Motions of Protein. Ph.D. Dissertation, Nanjing University, Nanjing, China, 2014. (66) Wilson, C. J.; Apiyo, D.; WittungStafshede, P. Role of Cofactors in Metalloprotein Folding. Q. Rev. Biophys. 2004, 37, 285−314. (67) Young, T. R.; Kirchner, A.; Wedd, A. G.; Xiao, Z. An Integrated Study of the Affinities of the Aβ16 Peptide for Cu(I) andCu(II): Implications for the Catalytic Production of Reactive Oxygen Species. Metallomics 2014, 6, 505−517. (68) Barry, S. D.; Rickard, G. A.; Pushie, M. J.; Rauk, A. The Affinity of HGGG, GHGG, GGHG, and GGGH Peptides for Copper(II) and the Structures of Their Complexes  an Ab Initio Study. Can. J. Chem. 2009, 87, 942−953. (69) Alred, E. J.; Nguyen, M.; Martin, M.; Hansmann, U. H. E. Molecular Dynamics Simulations of Early Steps in RNA-Mediated Conversion of Prions. Protein Sci. 2017, 26, 1524−1534.

5228


Binding of Copper Ions with Octapeptide Region in Prion Protein

Recommend Documents