RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

A Study on Fe2+ – ␣-Helical-Rich Keratin Complex Formation Using Isothermal Titration Calorimetry and Molecular Dynamics Simulation YANYAN ZHAO,1 JAN K. MARZINEK,2,3 PETER J. BOND,4,5 LONGJIAN CHEN,1 QIONG LI,1 ATHANASIOS MANTALARIS,2 EFSTRATIOS N. PISTIKOPOULOS,2 MASSIMO G. NORO,6 LUJIA HAN,1 GUOPING LIAN3 1

College of Engineering, China Agricultural University, Beijing 100083, People’s Republic of China Centre for Process Systems Engineering (CPSE), Department of Chemical Engineering, Imperial College London, London SW7 2BY, UK 3 Unilever Discover, Colworth Park, Sharnbrook, Bedford MK44 1LQ, UK 4 The Unilever Centre for Molecular Science Informatics, Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge CB2 1EW, UK 5 Bioinformatics Institute (A*STAR), Singapore 138671, Singapore 6 Physical and Chemical Insights Group, Unilever R&D, Port Sunlight, Wirral CH63 3JW, UK 2

Received 12 October 2013; revised 7 January 2014; accepted 17 January 2014 Published online 11 February 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI 10.1002/jps.23895 ABSTRACT: Iron binding to protein is common in biological processes of dioxygen transport, electron transfer as well as in stabilizing drug–protein complexes. ␣-Helix is the most prevalent secondary structure of proteins. In this study, Fe2+ binding to ␣-helix has been studied by isothermal titration calorimetry (ITC) and explicitly solvated molecular dynamics (MD) simulation. Ferrous gluconate and ␣helix-rich keratin are used for the ITC study and the results revealed followed one set of identical sites binding model. The MD simulations further revealed that only the acidic side-chain functional groups and ␩2 (O,O) coordination modes are involved in the binding of Fe2+ to ␣-helix. The ITC results also showed that the binding of ferrous gluconate to keratin was entropy driven and the higher the temperature, the stronger the binding free energy. The favorable entropy of Fe2+ binding to keratin was attributed to the displacement of water molecules on the ␣-helix surface, and was confirmed via MD simulations. The most stable coordination states of Fe2+ and ␣-helix were identified via simulation: Fe2+ stacks between two glutamic acid side chain carboxylate groups, displacing water molecules. The binding free energies C 2014 Wiley Periodicals, calculated using MD simulation and the theoretical values were in excellent agreement with the ITC results.  Inc. and the American Pharmacists Association J Pharm Sci 103:1224–1232, 2014 Keywords: ferrous ion; protein binding; ITC; molecular dynamics; ␣-helix; keratin; thermodynamics; complexation

INTRODUCTION Iron is essential in a number of fundamental biological processes, such as dioxygen transport, electron transfer, oxidation as well as stabilizing drug–protein complexes.1–3 A number of proteins are involved in the uptake, transport, storage, and export of iron.4–7 "-Helix is the most prevalent protein secondary structural type, accounting on average for 30% of globular proteins.8,9 Keratin is the major structural fibrous protein of skin, hair, and nails, and is mainly composed of "-helices hence used as an example in this study.10 The relative bioavailability of iron compounds is also determined by their solubility.4 Iron compounds of high water solubility, such as ferrous sulfate and ferrous gluconate, were proved to have higher bioavailAbbreviations used: ITC, isothermal titration calorimetry; MD, molecular dynamics; US, umbrella sampling; MW, molecular weight; SD, steepest descent; PMF, potential of mean force; SMD, steered molecular dynamics; Ka , binding constant; H, enthalpy change; S, entropy change; G, ree energy change. Correspondence to: Lujia Han (Telephone: +86-010-6273-6313; Fax: +86-10627-36-778; E-mail: [email protected]); Guoping Lian (Telephone: +44-01234222741; Fax: +44-1234-248010; E-mail: [email protected]) Yanyan Zhao and Jan K. Marzinek are the equal contributors. This article contains supplementary material available from the authors upon request or via the Internet at http://onlinelibrary.wiley.com/. Journal of Pharmaceutical Sciences, Vol. 103, 1224–1232 (2014)  C 2014 Wiley Periodicals, Inc. and the American Pharmacists Association

1224

ability than those of poor water solubility.4 Therefore, ferrous gluconate was chosen here to study the binding affinity with "-helical-rich keratin.4 The atomic level study of structural behavior, that is, binding sites and the role of solvent effect as well as the equilibrium thermodynamic properties of iron cations binding to "-helical motif are essential to understand iron homeostasis, detoxification, and the bioavailability. Further investigations are necessitated by a detailed description of the thermodynamics, that is, equilibrium binding properties (free energy, enthalpic, and entropic contribution), thus the maximum reachable concentration as well as detailed atomic conformations, that is, binding sites and the role of water in the binding process. Isothermal titration calorimetry (ITC) has previously been applied for studying metal ion binding to amino acids, peptides, and proteins in solution.4 ITC measurements provide the direct heat changes associated with binding between metal ions and protein. By building thermodynamic models of the ITC measurements, the number of binding sites as well as free energies with associated enthalpy and entropy contribution can be directly determined.11,12 Reported studies of such ferrous ion binding include Escherichia coli ferritin,13 recombinant human H-Chain ferritin,14 and yeast frataxin.15 Because of the advantage of yielding atomically detailed data regarding the assessment of binding sites and associated

Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

conformational changes within proteins, molecular dynamics (MD) simulations of ligand–protein interaction have recently provided complementary information.16–19 A study20 using quantum mechanical–molecular mechanical approaches reported the possibility of three different interactions between ferrous ion and an "-helical peptide, and 02 (O,O) coordination modes were reported as the most stable type. However, this simulation was based upon an ideal protein environment and neglected solvent effects. Both the binding thermodynamic properties and structural binding site assessment between "-helix and ferrous ion remain unknown. Hence, in this study, using water as a solvent, the interaction between ferrous gluconate and keratin has been investigated. To address the above thermodynamic and structural challenges, we have employed the experimental ITC method and MD simulations. MD studies revealed that only one type of binding sites was stable—acidic side chain oxygens. This is in agreement with the ITC data, which fitted well with the one set of identical binding sites model. The most stable conformation obtained via MD simulation was used for the free energy calculations via umbrella sampling (US) method and agreed within 3% of error with experimental results at three different temperatures. High entropy observed in ITC experiment has been explained as a water displacement within MD simulation framework. Both the number of binding types and binding free energies, predicted from the MD simulation, agreed extremely well with the ITC experimental results providing new insight of Fe2+ binding the most common "-helical motif.

MATERIALS AND METHODS Experimental Materials Ferrous gluconate dihydrate was obtained from J&K Science, Beijing, China. The purity was 98%. Keratin samples were separated from sheep wool fibers and kindly provided by Professor Jian Lu’s group at Manchester University, UK. The keratin in water solution was a mixture of two main components with molecular weights (MWs) of about 46 and 65 kDa. For the analysis of ITC data, a MW of 50 kDa was suggested by the provider. Water used was tap water filtered with a Milli-Q system (Millipore, Beijing, China). Sample Preparation Keratin samples were dissolved in water and the pH was measured to be at 6.08. The concentration of keratin solution was 0.5 mg/mL. By taking a MW of the keratin sample at 50 kDa, the concentration of keratin was 1.0 × 10−5 M. Ferrous gluconate solution (6.0 × 10−4 M) was prepared by dissolving ferrous gluconate dihydrate powder into pure water and the pH was measured to be 4.5. Ferrous gluconate was regarded as completely dissociated at pH 4.5. ITC Measurements A nanoisothermal titration calorimeter (Nano-ITC low volume TA) was used. The solution of keratin was placed in a 190 :L sample cell and ferrous gluconate solution was loaded into the injection syringe. Both keratin and ferrous gluconate were degassed. The titration schedule included addition of 2 :L per injection of titrant with 25 injections. Each injection time was 5 s with 300 s intervals. The solution in the DOI 10.1002/jps.23895

1225

cell was mixed at 300 rpm by the syringe. Control experiments were also performed for the titration of ferrous gluconate solution into water. In order to remove the dilution and hydrolysis effect, the heat of the control experiments was subtracted from the titration experiment. Experiments of of ferrous gluconate titration to keratin were carried out at 25◦ C, 35◦ C, and 45 ◦ C. Data analysis of ITC mesurements was conducted using Nano Analyze software provided by TA Instruments (159 Lukens Drive, New Castle, DE 19720, USA), including baseline-adjusted experimental titration data, peak-integrated background-subtracted and proper binding model fitting. MD Simulations—Computational Details The structure of ferrous gluconate was built using Chimera 1.5.3.21 In order to create the topology, the CHARMM General Force Field22–25 for organic molecules (program version 0.9.6 beta) was employed. The topology was then processed and translated into the GROMACS format using in-house code. The Lennard-Jones parameters of Fe2+ predicted for the TIP3P26 water model were taken from Li et al.27 The keratin sample used for the ITC study was a mixture of proteins of 46 and 65 kDa (type I and II) extracted from sheep wool. A survey of known keratin sequences from sheep wool was conducted. Type II keratin microfibrillar component 5 (http://www.uniprot.org/uniprot/P25691) with MW of 55 kDa is considered to be close to the average MW of the keratin sample. Interactions of ferrous ion with protein are mainly governed by electrostatic interaction. In wool, acidic and basic residues accounts for approximately 13% and 21%, respectively.28,29 Hence, we have chosen helical segment 1A with a similar content of charged amino acids (11.4% of acidic and 20% of basic). Helical rod accounts for over 61% (307 out of 502 residues) of wool keratin. The head and tail domains have very low content of acidic residues (seven out of 502 residues, 1.4%). Segment 1A of keratin was built using the PyMol software.30 The structure was then processed through the H++ online tool (http://biophysics.cs.vt.edu/H±±)31–33 to predict likely pKa values at pH 4.5, in accordance with our experiments. Ionizable keratin side chains with respect to the pKa were assigned their default, charged state. The energy of the initial structures was then minimized in vacuum via the steepest descent (SD) algorithm. All MD simulations were performed using the GROMACS 4.5.534 package with the CHARMM22/CMAP35 force field and the TIP3P water model. Classical equations of motion were integrated with a Verlet leap-frog algorithm using a 2 fs time step. The LINCS36 algorithm was employed to constrain bond lengths. A 1.4 nm cutoff distance was used for the short-range neighbor list, updated every five steps (10 fs), and van der Waals interactions were cut off at 1.4 nm. The Particle Mesh Ewald method37,38 was used for the long-range electrostatics with a 1.2 nm real space cutoff. The velocity rescale thermostat with additional stochastic term39 and Parinello–Rahman barostat40 were used to maintain the temperature and pressure (at 1 bar). The initial velocities were set to follow a Maxwell distribution. Periodic boundary conditions were used in all directions. All simulations were performed on two Linux clusters at Imperial College London: Imperial College High Performance Computing Service, URL: http://www.imperial.ac.uk/ict/services/teachingandresearchser Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

1226

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

Figure 1. Examples of initial configurations used in MD simulations for the binding site assessment: type I—Fe2+ interaction with backbone carbonyl oxygen (ALA67) and amide nitrogen (SER68); type II—Fe2+ interacting with side chain hydroxyl oxygen (SER68) and backbone carbonyl oxygen (ASP64); type III—Fe2+ coordinated by GLU71 side chain carboxylate oxygens. Cartoon representation of helix—blue. Atoms are shown as red oxygens, cyan carbons, gray hydrogens, dark blue nitrogens, and yellow Fe2+ .

vices/highperformancecomputing and the Department of Chemical Engineering Computing Service, URL: http://imper ial.ac.uk/chemicalengineering: MD Simulations—Simulation Setup Initially, a keratin segment of "-helix was placed in the center of a cubic box of 512 nm.3 The termini were treated in their uncharged state being a part of the 55 kDa keratin. A previous study hypothesized three possibilities for Fe2+ interaction modes with "-helical polypeptide fragments.20 This involved: type I, in which ferrous ion interacted with amide nitrogen and carbonyl oxygen atoms of the polypeptide backbone; type II, in which ferrous ion interacted with backbone carbonyl oxygens together with the side chain atoms of uncharged polar residues (including serine glutamine, asparginine, or threonine); and type III, in which only charged side chain functional groups (of aspartate and glutamate residues) were responsible for forming interactions. Therefore, to assess each of these possible "-helix binding site modes, Fe2+ were manually placed within 0.2 nm of polar atoms of each of the three mentioned possible interaction sites. In particular, type II sites involved either serine, asparagine or glutamine side chain oxygens in combination with neighboring residue backbone carbonyl oxygens. In the case of type III, the most likely binding sites were negatively charged acidic residues, and hence the Fe2+ was placed in proximity of side chain functional oxygens of either aspartic or glutamic acid carboxylate groups. In the case of the type I model, Fe2+ was placed in proximity of backbone polar atoms of different residues (polar, nonpolar, basic, and acidic). Examples of initial configurations (type I, II, and III) for ferrous–keratin complexes are shown in Figure 1. We focused on the assessment of the stability of each binding site in order to identify the potential binding pocket responsible for mediating the interaction of keratin with other small molecules. Because of the fact that ferrous gluconate is fully dissociated in aqueous solution, gluconic acid moieties, and Fe2+ were placed separately and then TIP3P water molecules filled the remainder of the box. The initial mass concentration of ferrous gluconate was set to 3%, which resulted in 20 molecules. Twelve ferrous ions were placed within three different hypothZhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

esized binding sites and the remaining ions and gluconic acid moieties were placed randomly around the keratin "-helical segment. To neutralize the excess system charge resulting from the excess presence of acidic keratin "-helix residues, one Na+ ion was added. Using the SD algorithm, energy minimization was then conducted to relax any steric overlap between solute and solvent. The system was then equilibrated in the canonical ensemble (constant number of atoms, volume, and temperature) for 200 ps and then in the isothermal–isobaric ensemble (constant number of atoms, pressure, and temperature) for another 200 ps, allowing only water and neutralizing ions to move freely. Production simulations were then run for 100 ns at three different temperatures: 25◦ C, 35◦ C, and 45◦ C. Free energy calculations were performed using the US method.41–43 In this technique, an additional bias potential (harmonic) is introduced to obtain sufficient sampling along a defined reaction coordinate (e.g., distance between two molecules), and subsequent unbiasing provides a potential of mean force (PMF) curve—the binding free energy as a function of the distance. Initial coordinates for US windows were obtained by steered MD (SMD) simulations, in which Fe2+ was gradually pulled away from the "-helical segment of keratin surface. A constant velocity was applied to Fe2+ with a harmonic potential (“virtual spring”) providing initial configurations for free energy calculations along the distance between the two species. A systematic assessment of the influence of different binding constants, as well as pulling velocities, on the force versus time profiles in the SMD simulations was conducted, and is presented with the results.

RESULTS ITC Studies on Fe2+ Binding The ITC results of ferrous gluconate titration at 25◦ C are shown in Figure 2 as a plot of heat flux against time on the top, and on the bottom as peak-integrated, background-subtracted, concentration-normalized molar heat flow per aliquot versus injection numbers. The control experiment for injection of ligands into buffer consisted of a series of equal heats of dilution. Ferrous gluconate dissociated after titration into water solution, DOI 10.1002/jps.23895

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

Figure 2. ITC measurement of ferrous gluconate (6.0 × 10−4 M, in pure water solution, 25◦ C) titrated into pure water.

and was an endothermic process. Figure 3 shows the results of titration of ferrous gluconate (6.0 × 10−4 M, dissolved in pure water solution at 25◦ C) into keratin solution (1.0 × 10−5 M, in pure water solution, 25◦ C). After subtracting the blank experiment of ferrous gluconate titration into water, the interaction of ferrous gluconate with keratin showed an exothermic process. The one set of identical binding sites model is used to fit the ITC experimental data, Q= ⎡

nMtHV0 2



1 Xt + − × ⎣1 + nMt nK a Mt

Figure 3. ITC measurement of ferrous gluconate (6.0 × 10 4 M, in pure water solution, 25◦ C) into keratin solution (1.0 × 10−5 M, in pure water solution, 25◦ C). —Observed enthalpy changes of ferrous gluconate and keratin after subtracting the control experiment.

⎤  2 Xt 1 4X t ⎦ 1+ + − (1) nMt nK a Mt nMt

 

dvi Q(i) − Q(i − 1) Q = Q(i) − Q(i − 1) + V 2

(2)

Here, Mt is the concentration of keratin in the cell, Xt is the concentration of ferrous gluconate in the cell, Q is the heat change of each injection, V0 is the cell volume (190 :L), and dvi is the injection volume for the i-th injection. The fitting parameters are enthalpy change (H), number of binding sites per mol of keratin (mol/mol) (n), and the binding constant (Ka ). Figure 4 shows this model fitted the data well. It also shows that increasing temperature resulted in an increase of the binding free energy. All the fitted parameters (n, Ka , and H) and calculated free energy change (G) and entropy change (S) by Eqs. 3 and 4 are shown in Table 1. G = −RT ln K a DOI 10.1002/jps.23895

1227

(3)

Figure 4. ITC measurement of ferrous gluconate titrated (6.0 × 10−4 M) into keratin solution (1.0 × 10−5 M) after subtracting the control experiment. Solid lines are the fitted enthalpy changes using Eqs. 1 and 2.

S =

H − G T

(4)

Standard errors of these parameters were estimated by the data fitting to Eqs. 1 and 2. Table 1 reveals that 1 mol keratin can bind approximately 9 (859 ± 0.04) mol of ferrous gluconate. Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

1228

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

Table 1.

Thermodynamic Parameters for Ferrous Gluconate Binding to Keratin at Different Temperatures of 25◦ C, 35◦ C, and 45◦ C

Temperature (◦ C) 25 35 45

n

H (kcal/mol)

TS (kcal/mol)

Ka

G (kcal/mol)

GMD (kcal/mol)

MAEa

8.59 ± 0.34 9.42 ± 0.24 9.22 ± 0.27

−0.86 ± 0.05 −1.05 ± 0.04 −1.12 ± 0.05

6.57 ± 0.14 6.64 ± 0.13 6.88 ± 0.18

2.81 ± 105 2.85 ± 105 3.17 ± 105

−7.43 ± 0.19 −7.69 ± 0.17 −8.00 ± 0.23

−7.20 ± 0.20 −7.60 ± 0.20 −7.90 ± 0.20

3.09 1.17 1.25

Data obtained by fitting the isothermal titration calorimetry data to Eqs. 1 and 2. Results from molecular dynamics simulations (GMD ) are also provided for comparison. Mean absolute error  (MAE) of the free energy between experiment and molecular simulations is below 3.1%.  G−GMD  a MAE =  × 100%. G

Figure 5. The most stable MD simulation of Fe2+ binding to "-helix at three suggested sites type I, II, and III. Type III binding site was shown to be the most stable, with the ion remaining bound to the acidic side chain carboxylate oxygen atoms for the whole 100 ns simulation.

The large positive entropy values obtained indicate that Fe2+ binding to keratin is entropically driven. MD Simulations

not remain for significant proportions of the simulation time, and hence are not further discussed here. Free Energy Calculations

Assessment of Binding Sites 2+

The three types likely binding sites (I, II, and III) for Fe binding to protein reported in the literature20 were assessed based on 100 ns simulation trajectories. In Figure 5, the most stable examples of Fe2+ binding to "-helical segment of keratin at type I, type II, and type III sites as a distance over the simulation time are presented. This distance corresponds to that between the center of mass of polar atom(s) of a given binding site and ferrous ion. It can be noted that only the distance of iron from type III binding site remained constant. The most stable interaction of molecular coordination structure is shown in Figure 6a, which corresponds to the type III in Figure 5. In one case, Fe2+ was placed within 0.2 nm of both the side chain carboxylate oxygens of GLU57 and the neighboring side chain hydroxyl oxygen of THR58, which corresponded to type II bind site. The corresponding snapshots together with the distance of Fe2+ from the center of mass of both GLU57 and THR58 side chain oxygens are presented in Figure 7. It is observed that from initial type II coordinates Fe2+ drifted to the conformation of type III binding site. More examples of type I, II, and III stability, in terms of the distance from the interacting site versus time, are presented in the supporting information (SI) in Figure S1. The gluconic acid moieties of ferrous gluconate appear to have very small affinity toward wool keratin. Of 20 initially placed gluconic acid anions, only 1–2 (at 25◦ C) or 3–5 (at 35◦ C) became stacked onto the "-helix surface via ionic interactions and hydrogen bonds (Fig. 6b). The mediating atom of ionic interactions between "-helix and the gluconate tended to be the coordinated ferrous ion. However, those weak interactions did Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

From the MD simulation, the final conformations of the Fe2+ – "-helix complex at type III binding mode (ferrous ion bound to glutamic acid side chain functional groups—Fig. 6a) were extracted. Initial coordinates for the US windows were obtained using SMD with a constant pulling rate. First, an assessment of the effect of different force constants at the same velocity was performed. For this purpose, Fe2+ was pulled 2 nm away from keratin "-helical segment, with a velocity of 10 nm/ns, with harmonic spring constants of: 960, 480, 240, 120, and 48 kcal mol−1 nm−2 . The force versus time profiles of these SMD simulations are presented in SI in Figure S2. Similarly, to assess the influence of different pulling rates, velocities of 1, 5, and 10 nm/ns were applied, with a spring constant of 480 kcal mol−1 nm−2 (Fig. S3). The SMD trajectory was used to extract initial coordinates for free energy calculations by US method. Since the system is largely governed by ionic interactions, extensive sampling was required in order to obtain a reasonable histogram overlap. Thus, a mean spacing of 0.1 nm with 100 ns sampling per US window was used along this coordinate. As a result, 22 US windows of 100 ns each were calculated at each of three temperatures. The harmonic spring constant in the US windows ranged between 120 and 480 kcal mol−1 nm−2 . The weighted histogram analysis method WHAM44,45 was then used to combine all windows and yield the PMF curve (free energy along the reaction coordinate). The first 10 ns of each window were omitted for system equilibration. In Figure 8, the PMF curves are shown at three different temperatures for Fe2+ binding to "-helix together with an example of US histograms, showing excellent overlap. DOI 10.1002/jps.23895

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

1229

DISCUSSION

Figure 6. Final conformations of ferrous gluconate with "-helix: (a) Fe2+ stacked in between two glutamic acid (GLU79 and GLU82) side chain carboxylate groups, coordinated by water molecules; (b) Fe2+ and water molecules bound to a GLU71 side chain carboxylate oxygen. Iron mediates ionic interactions with a gluconic moiety which also forms a hydrogen bond (red dotted line) with glutamic acid; cartoon representation of helix—blue; Fe2+ —yellow sphere; gluconate—green; protein atoms: oxygens (red), carbons (cyan), hydrogen (gray), and nitrogen (blue).

The free energy is the difference between the initial and the final state at plateau. The obtained free energies were G25◦ C = −7.2 kcal/mol, G35◦ C = −7.6 kcal/mol, and G45◦ C = −7.9 kcal/mol, yielding a striking agreement with the experimental ITC data (Table 1). The error was estimated to be ±0.2 kcal/mol and was assessed by dividing the whole sampling trajectory into five blocks of 20 ns and calculating the individual PMF curves. In order to show the reduced entropic cost for Fe2+ binding due to water displacement from the acidic type III binding site at increasing temperatures, the average number of hydrogen bonds between acidic side chains and water molecules as a function of US distance (Fig. 9) was calculated. The hydrogen bonds were collected based on a 0.35 nm cutoff distance between donor and acceptor atoms as well as a hydrogen-donor–acceptor angle of 30◦ , and then averaged over the 100 ns simulation time of each window. DOI 10.1002/jps.23895

As previously mentioned, a recent study20 suggested three possible binding types between Fe2+ and protein amino-acid residues of ferritin (as illustrated in Fig. 1). However, the fact that our ITC results are best fitted by a model involving one set of identical binding sites indicates only one type of interaction dominates the binding of ferrous ion to "-helix-rich keratin (Fig. 4). It is possible that the other two types of interaction are too weak to be measured under the dilute solution conditions by ITC.13 The results from the molecular simulation confirmed that only ferrous ions at type III side chain carboxylate oxygens of both glutamic and aspartic acid remained stable for the whole simulation time (Figs. 5 and S1). Ferrous ions at type I and II locations drifted away, suggesting weak interactions. This is further confirmed in Figure 7: Fe2+ remained only 7 ns in the type II configuration and then drifted toward the side chain oxygens of GLU57, remaining exclusively coordinated by the carboxylate group until the end of the 100 ns simulation. In this case, the distance of Fe2+ from the GLU57 carboxylate oxygen atoms remained constant for the whole simulation time. The final coordinates of the most stable conformation were extracted for free energy calculations. During the assessment of force constants and pulling rates (Fig. S2), the weaker force constants (120 and 48 kcal mol−1 nm−2 ) were insufficient to extract Fe2+ from the keratin type III acidic pocket, and a force breaking point was not reached. With the stronger force constants (240, 480, and 960 kcal mol−1 nm−2 ), a maximum pulling force corresponding to the breaking point was obtained. This rupture force corresponds to the point when the ionic bond between the ion and the acidic carboxylate group was broken. On the basis of Figure S2, it is worth noting that with an increase of the strength of the spring, the breaking points were reached within shorter times, but the results were otherwise comparable. Thus, an intermediate spring constant of 480 kcal mol−1 nm−2 was chosen for subsequent pulling simulations. Application of different pulling velocities produced identical force versus time profiles (Fig. S3). Hence, the fastest pulling rate of 10 nm/ns (with a spring constant of 480 kcal mol−1 nm−2 ) was chosen to produce the initial coordinates for US windows. The results from free energy calculations by US are represented by PMF curves in Figure 8 at three different temperatures. The first PMF minimum corresponds to the configurations in which Fe2+ strongly interacted with acidic carboxylate oxygens of "-helix. The second local minimum at approximately 0.4 nm represents the first coordination shell when ferrous ion was cross-linked to the acidic side chain via a water molecule, as confirmed by the step-wise increase in hydrogen-bonds between "-helix and water at this distance (Fig. 9). The PMF then increases with distance, and reaches a plateau at around ∼1– 1.5 nm. It is observed that increasing temperature resulted in an increase of the binding energy in both experiments and MD simulations (Table 1). A previous ITC study on Fe2+ binding to Human H-Chain Ferritin also demonstrated large positive entropy changes.14 The free energy obtained by their ITC measurments was nearly identical (−7.09 kcal/mol, 25◦ C, pH 6.5) in comparison with our result (−7.43 kcal/mol, 25◦ C, in pure water solution).14 The large entropy change was hypothesized to be due to the displacement of water molecules bound to protein acidic residues46 and was further confirmed by the MD simulation. With the increase of the kinetic energy (temperature), there was a greater Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

1230

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

Figure 7. Fe2+ binding to "-helix at type II binding site. After 7 ns, Fe2+ moved away from type II binding site, initially located between GLU57 and THR58, to type III binding site occupying only carboxylate oxygens of GLU57, and remained there until the end of the 100 ns simulation. The distance with the carboxylate oxygen atoms remained stable throughout. Atoms are shown as red oxygens, cyan carbons, gray hydrogens, dark blue nitrogens, and yellow Fe2+ .

Figure 8. (a) Potential of mean force (PMF) curves for Fe2+ –"-helix polypeptide segment of keratin interaction at 25◦ C, 35◦ C, and 45◦ C, together with (b) an example of US histograms at 25◦ C.

propensity to release the water molecules from the "-helix surface, leading to the higher entropic contribution and thus the greater magnitude in the binding free energies. The average number of hydrogen bonds between acidic pocket and water molecules within US windows in Figure 9 confirms our entropic findings. It is observed that within 0.4 nm, the acidic side chains

Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

were occupied by ferrous ion, thus blocking access to water. The release of ferrous ion from the "-helix surface at distances greater than 0.4 nm was associated with an increasing number of hydrogen bonds between the acidic carboxylate oxygens and water molecules. Once the distance between ferrous ion and "-helix surface became greater than 0.7 nm, the average DOI 10.1002/jps.23895

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

1231

Figure 9. The average number of hydrogen bonds between acidic side chain oxygens and water molecules at 25◦ C, 35◦ C, and 45◦ C as function of the distance of Fe2+ and keratin "-helix acidic binding site. Hydrogen bonds were calculated from US windows based on a 0.35 nm cutoff distance between donor and acceptor atoms as well as a hydrogen-donor–acceptor angle of 30◦ , and then averaged over the 100 ns simulation time of each window.

number of hydrogen bonds remained constant. As higher temperatures resulted in a greater mobility of solvent molecules, the lower stability of associated hydrogen bonds with protein was observed. This explains the greater entropic contribution associated with the displacement (release) of water during the binding of ferrous ion at higher temperatures.

CONCLUSIONS The Fe2+ and "-helix-rich keratin binding in water environment is studied by both the ITC experiments and MD simulations. The results suggested that ferrous gluconate interaction with "-helix-rich keratin is dominated by the binding of ferrous ion to acidic side chain carboxylate group oxygens. ITC experiments also revealed that the interaction is exothermic and 1 mol of keratin can bind 9 mol of ferrous gluconate. Higher temperatures facilitated the binding process, suggesting that the interaction is entropy driven. This is because of the increased mobility of solvent molecules and hence lower stability of hydrogen bonds of ferrous ion with the acidic pocket. A greater displacement of water from the "-helix surface upon Fe2+ binding was observed at lower temperatures, confirming the importance of the entropic contribution to the binding process indicated by both the ITC experiments. The most stable coordination state of ferrous ion binding to "-helix observed during MD simulations was obtained when Fe2+ stacked between pairs of glutamic acid side chain carboxylate groups, with displaced water molecules bound around the cation, blocking their direct access to the "helix surface. The free energies calculated by MD simulations are in excellent agreement with the ITC experimental results, in further support of the predicted binding mode.

ACKNOWLEDGMENTS This research is financially supported by Unilever R&D Colworth, National Natural Science Foundation of China (Project No. 21006124), Program for New Century Excellent TalDOI 10.1002/jps.23895

ents in University (Project No. NCET-11–0477), Program for Changjiang Scholars and Innovative Research Team in University (Project No. IRT1293), and the EC’s Seventh Framework Programme (FP7/2007–2013) under Grant Agreement No. 238013, which are greatly acknowledged. Yanyan Zhao wishes to thank Professor Jian Lu from Manchester for supplying keratin samples. Jan Marzinek wishes to acknowledge Dr. Robert Farr and Matt Harvey for stimulating discussions as well High Performance Computing service at Imperial College London.

REFERENCES 1. Smith SJ, Du K, Radford RJ, Tezcan FA. 2013. Functional, metalbased crosslinkers for [small alpha]-helix induction in short peptides. Chem Sci 4(9):3740–3747. 2. La Mendola D, Magr`ı A, Campagna T, Campitiello MA, Raiola L, ¨ Bonomo RP, Rizzarelli E. 2010. A doppel "-helix Isernia C, Hansson O, peptide fragment mimics the copper(II) interactions with the whole protein. Chem Eur J 16(21):6212–6223. 3. Xu J. 2005. An insight into interaction of Fe2+ with glycylglycine: A DFT study. J Mol Struct: THEOCHEM 757(1–3):171–174. 4. Wilcox DE. 2008. Isothermal titration calorimetry of metal ions binding to proteins: An overview of recent studies. Inorg Chim Acta 361(4):857–867. 5. Perron N, Brumaghim J. 2009. A review of the antioxidant mechanisms of polyphenol compounds related to iron binding. Cell Biochem Biophys 53(2):75–100. 6. Sacco C, Skowronsky R, Gade S, Kenney J, Spuches A. 2012. Calorimetric investigation of copper(II) binding to A$ peptides: Thermodynamics of coordination plasticity. J Biol Inorg Chem 17(4):531–541. 7. Grossoehme NE, Akilesh S, Guerinot ML, Wilcox DE. 2006. Metalbinding thermodynamics of the histidine-rich sequence from the metal-transport protein IRT1 of Arabidopsis thaliana. Inorg Chem 45(21):8500–8508. 8. Baker EN, Hubbard RE. 1984. Hydrogen bonding in globular proteins. Prog Biophys Mol Biol 44(2):97–179. 9. Nick Pace C, Martin Scholtz J. 1998. A helix propensity scale based on experimental studies of peptides and proteins. Biophys J 75(1):422– 427. Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

1232

RESEARCH ARTICLE – Pharmaceutics, Drug Delivery and Pharmaceutical Technology

10. Vasconcelos A, Freddi G, Cavaco-Paulo A. 2008. Biodegradable materials based on silk fibroin and keratin. Biomacromolecules 9(4):1299– 1305. 11. Zhao Y, Chen L, Yakubov G, Aminiafshar T, Han L, Lian G. 2012. Experimental and theoretical studies on the binding of epigallocatechin gallate to purified porcine gastric mucin. J Phys Chem B 116(43):13010–13016. 12. Zhao Y, Marzinek JK, Chen L, Mantalaris A, Pistikopoulos EN, Han L, Lian G, Bond PJ, Noro MG. 2013. Molecular and thermodynamic basis for EGCG-keratin interaction—Part II: Experimental investigation. AIChE J 59(12):4824–4827. ´ 13. Bou-Abdallah F, Woodhall MR, Velazquez-Campoy A, Andrews SC, Chasteen ND. 2005. Thermodynamic analysis of ferrous ion binding to Escherichia coli ferritin EcFtnA†. Biochemistry 44(42):13837– 13846. 14. Bou-Abdallah F, Arosio P, Santambrogio P, Yang X, JanusChandler C, Chasteen ND. 2002. Ferrous ion binding to recombinant human H-chain ferritin. An isothermal titration calorimetry study. Biochemistry 41(37):11184–11191. 15. Cook JD, Bencze KZ, Jankovic AD, Crater AK, Busch CN, Bradley PB, Stemmler AJ, Spaller MR, Stemmler TL. 2006. Monomeric yeast frataxin is an iron-binding protein. Biochemistry 45(25):7767–7777. 16. Agnieszka KB. 2011. Thermodynamics of Ligand-Protein Interactions: Implications for Molecular Design, Thermodynamics - Interaction Studies - Solids, Liquids and Gases, Dr. Juan Carlos Moreno ˜ (Ed.). InTech, DOI: 10.5772/19447. PirajA¡n 17. Marzinek JK, Zhao Y, Lian G, Mantalaris A, Pistikopoulos EN, Bond PJ, Noro MG, Han L, Chen L. 2013. Molecular and thermodynamic basis for EGCG-keratin interaction—Part I: Molecular dynamics simulations. AIChE J 59(12):4816–4823. 18. Feenstra P, Brunsteiner M, Khinast J. 2012. Prediction of drugpackaging interactions via molecular dynamics (MD) simulations. Int J Pharm 431(1–2):26–32. 19. Woo H-J, Roux B. 2005. Calculation of absolute protein–ligand binding free energy from computer simulations. Proc Natl Acad Sci USA 102(19):6825–6830. 20. Jurinovich S, Degano I, Mennucci B. 2012. A strategy for the study of the interactions between metal–dyes and proteins with QM/MM approaches: The case of iron–gall dye. J Phys Chem B 116(45):13344– 13352. 21. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, Ferrin TE. 2004. UCSF chimera—A visualization system for exploratory research and analysis. J Comput Chem 25(13):1605–1612. 22. ParamChem Interface. Accessed 2013, at: https://www.paramchem. org. 23. Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, Shim J, Darian E, Guvench O, Lopes P, Vorobyov I, Mackerell AD. 2010. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem 31(4):671–690. 24. Vanommeslaeghe K, MacKerell AD. 2012. Automation of the CHARMM General Force Field (CGenFF) I: Bond perception and atom typing. J Chem Inf Model 52(12):3144–3154. 25. Vanommeslaeghe K, Raman EP, MacKerell AD. 2012. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of bonded parameters and partial atomic charges. J Chem Inf Model 52(12):3155–3168. 26. Berendsen HJC PJ, van Gunsteren WF, Hermans J. 1981. Interaction models for water in relation to protein hydration. Intermol Forces B14:331.

Zhao et al., JOURNAL OF PHARMACEUTICAL SCIENCES 103:1224–1232, 2014

27. Li P, Roberts BP, Chakravorty DK, Merz KM. 2013. Rational design of particle mesh ewald compatible Lennard-Jones parameters for +2 metal cations in explicit solvent. J Chem Theory Comput 9(6):2733– 2748. 28. Corfield MC, Robson A. 1955. The amino acid composition of wool. Biochem J 59(1):62–68. 29. Sahoo A SN. 2011. Nutrition for wool production. WebmedCentral Nutrition 2(10):WMC002384. 30. The PyMOL molecular graphics system, Version 1.5.0.4. Schr¨odinger, LLC. http://www.pymol.org/citing 31. Anandakrishnan R, Aguilar B, Onufriev AV. 2012. H++ 3.0: Automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations. Nucleic Acids Res 40(W1):W537-W541. 32. Myers J, Grothaus G, Narayanan S, Onufriev A. 2006. A simple clustering algorithm can be accurate enough for use in calculations of pKs in macromolecules. Proteins 63(4):928–938. 33. Gordon JC, Myers JB, Folta T, Shoja V, Heath LS, Onufriev A. 2005. H++: A server for estimating pKas and adding missing hydrogens to macromolecules. Nucleic Acids Res 33(suppl 2):W368-W371. 34. Hess B, Kutzner C, van der Spoel D, Lindahl E. 2008. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J Chem Theory Comput 4(3):435–447. 35. Bjelkmar Pr, Larsson P, Cuendet MA, Hess B, Lindahl E. 2010. Implementation of the CHARMM Force Field in GROMACS: Analysis of protein stability effects from correction maps, virtual interaction sites, and water models. J Chem Theory Comput 6(2):459–466. 36. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM. 1997. LINCS: A linear constraint solver for molecular simulations. J Comput Chem 18(12):1463–1472. 37. Darden T, Darrin Y, Pedersen L. 1993. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. J Chem Phys 98(12):10089–10092. 38. Essmann U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG. 1995. A smooth particle mesh Ewald method. J Chem Phys 103(19):8577–8593. 39. Bussi G, Donadio D, Parrinello M. 2007. Canonical sampling through velocity rescaling. J Cheml Phys 126(1):014101–014107. 40. Parrinello M RA. 1980. Phys Rev Lett 45:1196–1199. 41. Torrie GM, Valleau JP. 1974. Monte Carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical LennardJones fluid. Chem Phys Lett 28(4):578–581. 42. Torrie GM, Valleau JP. 1977. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J Comput Phys 23(2):187–199. 43. Buch I, Sadiq SK, De Fabritiis G. 2011. Optimized potential of mean force calculations for standard binding free energies. J Chem Theory Comput 7(6):1765–1772. 44. Hub JS, de Groot BL, van der Spoel D. 2010. g wham—A free weighted histogram analysis implementation including robust error and autocorrelation estimates. J Chem Theory Comput 6(12):3713– 3720. 45. Kumar S, Rosenberg JM, Bouzida D, Swendsen RH, Kollman PA. 1992. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J Comput Chem 13(8):1011–1021. 46. Sundstr¨om M, Hall´en D, Svensson A, Schad E, Dohlsten M, Abrahms´en L. 1996. The co-crystal structure of staphylococcal entero˚ resolution: Imolications for major histoxin type A with Zn2+ at 2.7 A tocompatibility complex class II binding. J Biol Chem 271(50):32212– 32216.

DOI 10.1002/jps.23895