R. Sandya Rani and
Moumita Saharay*
Department of Physics, Osmania University, Hyderabad, India. E-mail: saharaymou@gmail.com
First published on 16th January 2019
The protein-mediated biomineralization of calcium carbonate (CaCO3) in living organisms is primarily governed by critical interactions between the charged amino acids of the protein, solvent, calcium (Ca2+) and carbonate (CO32−) ions. The present article investigates the molecular mechanism of lysozyme-mediated nucleation of amorphous calcium carbonate (ACC) using molecular dynamics and metadynamics simulations. The results reveal that, by acting as nucleation sites, the positively charged side chains of surface-exposed arginine residues form hydrogen bonds with carbonates and promote aggregation of ions around them leading to the formation and growth of ACC on the protein surface. The newly formed ACC patches were found to be less hydrated due to ion aggregation-induced expulsion of water from the nucleation sites. Despite favorable electrostatic interactions of the negatively charged side chains of aspartate and glutamate with calcium ions, these residues contribute minimally to the growth of ACC on protein surface. The activation barrier for the growth of partially hydrated ACC patches on lysozymes was determined from the free energy profiles obtained from metadynamics simulations.
In living organisms, the process of mineral formation is primarily governed by the interactions between the charged biopolymers,9–12 ionic constituents (Ca2+ and CO32−), and water. For instance, the negatively charged proteins in corals and sea urchin spicules are shown to stabilize the amorphous phase of CaCO3 even in in vitro experiments in the absence of any biological cells.13–15 A recent study on the formation of coral skeleton demonstrated that the ACC phase is indeed the metastable precursor to the crystalline CaCO3, and that the growth rate of aragonite crystals through the non-classical crystallization process11,12 involving aggregation of ACC particles is ∼100 times faster than the classical route of ion-by-ion formation of crystals.16 Earlier studies have also shown that the binding affinity of Ovocleidin-17 (OC-17), a hen eggshell protein, for ACC nanoparticles is much higher than that for the calcite surface17,18 suggesting that favorable binding of the protein to ACC is critical for biomineralization and eggshell formation. The lower calcite-protein binding affinity is leveraged for the detachment of the protein from the crystal after the ACC to crystal transition.18 Thus, OC-17 or a class of functionally similar proteins/biopolymers can be regarded as catalysts to facilitate the crystallization of CaCO3.
Although previous computer simulation studies have investigated the mechanism of the attachment of the protein19,20 to the ACC nanoparticles and CaCO3 crystal surfaces, little is known about the role of protein in promoting the formation of ACC from the ionic solution. Another key feature that is yet to be understood is the atomistic details of the mechanism and energetics of protein-ACC association/dissociation process. Some proteins/biopolymers hinder the crystallization of CaCO3 by stabilizing3,21,22 the amorphous phase, and the biomineralization is activated as and when the living organism requires it. The nucleation of CaCO3 in the presence of negatively charged biopolymers10 is thought to occur in two steps: (i) the binding of Ca2+ with the biopolymer to form Ca2+-bound globules, (ii) the diffusion and association of CO32− ions with the globules to form polymer-mediated ACC phase that then grows in size. These polymers, referred to as ‘ion sponge’, critically determine where and how the CaCO3 nucleation will take place and proceed further. However, the precise molecular details associated with these steps are still missing.
In vitro experiments have shown that the positively charged lysozyme increases the nucleation rate of CaCO3 crystals.23,24 The large negative ζ-potential (−14 mV) around this protein favors the formation of negatively charged double layers of CO32− around its surface.23 To understand the microscopic details of the onset of protein-mediated aggregation of Ca2+ and CO32− and to identify the vital molecular attributes and the key residues that govern the mineralization process,25,26 we have performed a series of atomistic molecular dynamics simulations of lysozyme, a hen eggshell protein. Our results reveal that the positively charged side chains of solvent-exposed arginine residues act as nucleation sites for the formation and growth of ACC on protein surface, while the negatively charged side chains of glutamate and aspartate residues contribute minimally to biomineralization.
System | CaCO3 pairs | Water molecules | Orthorhombic box dimension (Å3) |
---|---|---|---|
LW | — | 30284 | 70.0 × 70.0 × 70.0 |
S1LC | 235 | 14289 | 76.2 × 76.2 × 76.2 |
S2LC | 474 | 14238 | 76.5 × 76.5 × 76.5 |
hACC | 1620 | 1094 | 50.0 × 50.0 × 50.0 |
The LW model allows us to characterize the lysozyme–water interactions in the absence of Ca2+ and CO32− ions, whereas S1,2LC models are used to probe the influence of these ions on lysozyme. The initial coordinates of the hen egg-white lysozyme (PDB Id: 2cds) was taken from the X-ray crystal structure. To neutralize the net positive charge on lysozyme, 8 Cl− were added in LW system, while 4 extra CO32− were added in S1,2LC models. The initial configurations of S1,2LC were generated by randomly placing desired numbers of Ca2+ and CO32− ions in cubic water boxes of sizes specified in Table 1. The lysozyme was then solvated in this ionic solution by removing the ions and water molecules that were in hard contact with the protein using CHARMM-GUI27–29 package. CHARMM36 (ref. 30) force field (FF) parameters were used for the protein and ions (Ca2+, CO32−) together with the standard TIP3P water model.
For hACC, H2O:CaCO3 ratio was maintained at 0.67 to be consistent with our earlier investigation.31 The MD simulations were performed using NAMD32 package. The systems were equilibrated for 1 ns in NPT ensemble followed by 15 ns of MD simulations in NVT ensemble with the time step of 2 fs. The temperature and pressure were maintained at 303 K and 1 atm, respectively. The temperature was controlled by the Langevin thermostat using a damping coefficient of 1 ps−1. The electrostatic interactions were evaluated using particle mesh Ewald method with a grid spacing of 1 Å. Periodic boundary conditions were used and a distance cutoff of 12 Å was chosen to calculate the non-bonded interaction energies between different species.
To understand the flexibility of the near neighbor environment around the positively charged side chains, we studied the evolution of hydrogen bonds formed by the side chains with the nearest CO32− and water molecules. An arginine side chain was considered to be hydrogen bonded to CO32− if the CZ–Oc (refer Fig. 1 for atom labels) distance was less than 4.5 Å or to a water molecule if the CZ–Ow distance was less than 4.0 Å.33–35 For each type of h-bond, the time correlation function, Chb(t), which is defined as follows,
C(t) = 〈h(0)h(t)〉/〈h〉 |
Fig. 1 Atom names used for the descriptions of side chains in (a) arginine, (b) aspartate, and (c) glutamate. |
We performed well-tempered metadynamics36–38 simulations, using NAMD package, to estimate the free energy required for the association of ions and dehydration of hydrated Ca2+–CO32− clusters around selected side chains of lysozyme. The separation (rCV) between the tagged ions and the residue of interest was chosen as the collective variable (CV). The Gaussian hills were added at a regular interval of 100 fs, and the width of each Gaussian hill was 1 Å. The total run length for the metadynamics simulation was 100 ns.
Fig. 2 Arrangement of Ca2+ (red spheres) and CO32− (black) ions around arginine (blue), aspartate (green), and glutamate (green) residues in S1LC for (a) initial configuration and (b) after 15 ns of MD simulation. Water molecules are not shown for clarity. VMD39 software package was used for visualization. |
The clusters of ions are observed to grow in size during the course of MD simulation. The largest ion cluster is formed around ARG14, ARG125, and ARG128. Due to their close proximity, these three surface residues form an ‘attractive triangle’, which favors the rapid growth of amorphous calcium carbonate on lysozyme. The positively charged guanidinium (Gdm+) functional groups of arginine residues play a crucial role in the nucleation and growth25 of crystalline CaCO3 by favorably interacting with the negatively charged CO32− ions. The Gdm+–CO32− and Ca2+–CO32− interactions facilitate the ion association and further growth of amorphous CaCO3-like domains around these residues.
g(r) | S2LC | hACC | ||
---|---|---|---|---|
Atom pair | Peak position r (Å) | N(r) | Peak position r (Å) | N(r) |
Cc–Cc | 4.8 | 3.0 | 4.6 | 9.2 |
Ca–Cc | 2.6, 3.2 | 2.3 | 2.6, 3.2 | 4.6 |
Ca–Ow | 2.2 | 3.0 | 2.2 | 0.8 |
Ca–Ca | 3.3, 3.8 | 1.3 | 3.3, 3.7 | 4.4 |
Ca–Oc | 2.1 | 3.2 | 2.07 | 5.4 |
Cc–Ow | 3.3 | 12.8 | 3.1 | 4.04 |
Although the peak intensities in ion–ion RDFs are relatively higher in S2LC than in hACC, peak positions remain mostly unchanged. The ion–ion RDFs indicate that the ions are held together strongly to form a semi-rigid and well-structured coordination shells in S2LC. The present structural analysis reveals that local environment around lysozyme favors the aggregation of ions and the formation of amorphous calcium carbonate, although, the concentration of ions is much lower than the hydrated ACC (hACC) studied here.
In Fig. 4, the RDFs between NGdm (guanidinium-nitrogen) and Cc (carbonate carbon) atom pair shows three maxima at 3.34, 5.18, and 7 Å. The coordination number at the first minimum is 1.24. The first and second maxima are separated by a plateau of width ∼0.7 Å with negligible intensity, showing a crystal-like ordering of Cc atoms in this region. The number of Oc (carbonate oxygen) around NGdm is ∼1.4 at the first minimum (2.9 Å) of gNGdm–Oc(r) (Fig. 4). The underlying geometry could possibly indicate a configuration where one of the H atoms of a NGdm forms single h-bond with Oc. This preferential arrangement compares well with the geometry optimized isolated cluster25 of Gdm+–CaCO3–2H2O, where the distance between NGdm and Oc is ∼2.65 Å. The second and third maxima in gNGdm–Oc(r) appear at 3.4 and 4.6 Å, respectively. The coordination numbers at the corresponding minima in S1LC are 2.3 (2.0 in S2LC) and 4.7 (4.5 in S2LC), respectively.
The first maximum in gNGdm–Ca(r) (Fig. 4) shows two sub peaks at 4.2 and 5.5 Å in both the systems. This characteristic arrangement of Ca2+ around NGdm is somewhat similar to that of Ca2+-environment around Cc in hACC (Fig. 3), where the first peak of gCa–C(r) is also divided into two sub peaks. This indicates two preferred arrangements: for longer distance, Ca2+ is coordinated with only one NGdm or Oc making a ‘monodentate’ geometry whereas for shorter distance Ca2+ is equidistant from two NGdm or Oc making a ‘bidentate’ geometry. With increased ionic concentration in S2LC, the intensity of these peaks decreases which implies a less-structured coordination shell with same number of neighboring Ca2+.
The Cc coordination number around CZ (Fig. 4) is 2 within the first coordination shell for rCZ–Cc < 5.4 Å, which indicates that each pair of NGdm atoms (NH1, NH2; and NH2, NE) of Gdm+ interacts with a single CO32−. The arrangement and density of anions up to the second coordination shell is mostly similar in both the systems while the ion density increases significantly for S2LC beyond second coordination shell.
However, the near neighbor environment around the negatively charged functional groups (Fig. 5) of ASP and GLU residues are mostly populated by Ca2+. The RDF between carbonyl-oxygen (OD) of ASP and Ca2+ shows two sharp maxima at 2.2 and 4.2 Å, separated by a zero-intensity plateau of 1 Å width which signifies a crystal-like near neighbor ordered structure up to ∼5 Å from ASP. In addition, the gCG–Ca(r) (Fig. 5) also shows two high-intensity peaks at ∼2.2 and 3.2 Å signifying two preferred geometries of the cations with respect to the COO− group of ASP/GLU, bidentate and monodentate respectively. The radial distributions of Cc and Oc (Fig. 5) around COO− group of ASP/GLU are poorly structured with coordination numbers less than 1 up to the first minima at ∼5.0 Å for both the distributions.
The calculated atomic density distribution (Fig. 6), i.e. density map, of CO32− around Gdm+ shows two density peaks of Cc atoms within the CZ–Cc distance cutoff of 4.3 Å, the first minimum in gCZ–Cc(r) (Fig. 4). In contrast, Oc density map within the first minimum of gCZ–Oc(r) at 3.5 Å shows four high-density regions along the N–H bond directions.
In general, the non-polar compounds approach the arginine side chain from the directions normal to Gdm+ plane, while the polar compounds, e.g. water, carboxylates/carbonyls, form in-plane h-bonds40,41 with Gdm+. Our results reveal that in the lysozyme–CaCO3 systems, the CO32− ions replace the water molecules from the first coordination shell of arginine to form the TNTO coordination geometry. The second coordination shell is primarily formed by the Ca2+ ions with approximately three near-neighbor (NN) cations for each CO32− (Fig. 6). The observed distribution of NN Ca2+ and CO32− around arginine facilitates the formation of further layers of ions leading to their aggregation around the surface-exposed arginine side chains of lysozyme. In S1,2LC, the coordination number of water oxygen (Ow) around arginine is less than that of in LW. The aggregation of ions around Gdm+ in lysozyme–CaCO3 systems by expelling water molecules indicates that arginine residues act as nucleation sites for the formation of biominerals. The observed water expulsion at the protein-ACC interface is in-line with Freeman et al.18's work that showed that the water population at the ACC-protein interface is much less than that at the calcite-protein interface.
The binding energies of Ca2+–CO32− and Ca2+–H2O complexes are ∼520.8 kcal mol−1 (ref. 25) and 56.9 kcal mol−1,42 respectively. Among these two attractive interactions, the hydrophilic nature of calcium dominates which results in the separation of the CO32− from its Ca2+ counterpart via a water layer (Fig. 7) in between. These water molecules are ‘tightly bound’ by both types of ions. As the life time of Ca–H2O h-bond is large, the possibility of breaking this h-bond within the simulation time limit is slim. Thus, the removal of water from the nearest hydration layer and aggregation of ions around aspartate/glutamate seems unlikely.
A convoluted picture of the distribution of ions and the ion density with reference to a given side chain residue are presented here by charge density plots (Fig. 8). The average charge density (ρ) around arginine demonstrates sharp oscillations around the neutral axis, ρ = 0, up to ∼7 Å, indicating the formation of four alternating layers of cations and anions. However, the charge density around the COO− bearing groups (aspartate/glutamate) shows only one prominent positive peak with a much weaker next layer of anions before the distribution merges to the neutral axis. Thus, we can conclude that arginine may play the major role in the nucleation and growth of CaCO3 crystal in lysozyme, whereas the other two residues may facilitate the protein to bind to the positively charged surface of the CaCO3 crystal.
Each time correlation function was fitted with a multiexponential function to obtain average time constant or relaxation time (τhb, Table 3). The values of τhb for CO32− are orders of magnitude higher than those for water, which is indicative of higher residence times of CO32− ions in the vicinity of the nucleation sites. The shorter relaxation times of ‘loosely bound’ water molecules may be necessary to facilitate the dehydration of the ACC patches on the protein surface.
System | Function | τhb (ns) |
---|---|---|
S1LC | Cchb(t) | 10.10 |
Cwhb(t) | 0.69 | |
S2LC | Cchb(t) | 154.53 |
Cwhb(t) | 0.12 |
The calculated free energy profile estimates the energy cost to bring an anion or a Ca2+–CO32− ion pair that is located away from the protein surface to the vicinity of arginine side chain. The separation (rCV) between the Cc atom of a distant CO32− and the CZ atom of ARG61 was chosen as a collective variable (CV) to describe arginine–Ca2+–CO32− association. Fig. 10 shows the calculated free energy profile as a function of rCV in the range from 6.0 Å to 20.0 Å.
Fig. 10 (a) Free energy profile for the association of CO32− around the arginine side chain along the reaction coordinate (rCV). (b) Snapshot from the global minimum of the free energy profile showing the arrangement of ions and water. Color scheme: same as in Fig. 6. |
In the initial configuration, the chosen CO32− ion pairs with a neighboring Ca2+ to form an isolated CaCO3 unit surrounded by water. This particular Ca2+ ion was observed to be bound to the selected CO32− throughout the range of rCV. The free energy profile shows a shallow and broad minimum at rCV = 11.4 Å (Fig. 10a). The binding energy (ΔF) between the dissociated state at rCV = 20 Å and the global minimum is ∼8.4 kcal mol−1. This value compares well with the binding energy (∼10 kcal mol−1) between OC-17 (a hen eggshell protein) and calcite surface, where arginine residues serve as the anchors to the crystal surface.18 However, this barrier height (ΔF) is 4 kcal mol−1 higher than that estimated for a single CO32− around a Gdm+ functional group25 and that between an isolated arginine residue and ACC,26 which can be attributed to the inhomogeneous environment of the carbonate ion in the present study. The interactions between an anchoring agent (here, arginine) and the ions (Ca2+ and CO32−) play a crucial role in the initial stages of the ion aggregation, but further growth of ACC is generally governed mainly by the interactions between ions and water molecules. Thus, the barrier height reported here is an approximate estimate on the energetics for the growth of existing ion aggregates irrespective of the anchoring agent.
Further analysis of the metadynamics trajectory shows that the addition of this single CO32− in the arginine environment increases the ion-cluster size by at least three more CaCO3 pairs via near-neighbor electrostatic interactions (Fig. 10b).
Two model systems of lysozyme were generated with different ion concentrations, and both show essentially similar results, apart from a denser ACC phase in the solution with higher concentration of ions. The three dimensional density profile of ions around arginine shows that, on the average, each arginine is coordinated by at least two CO32− ions forming four h-bonds with the Gdm–NH2-groups making a twin-nitrogen-twin-oxygen (TNTO) geometry. This strong association with the nearest neighbors initializes the foundation for the successive alternating layers of positive and negative charges replacing the solvent water molecules. Moreover, the charge density plots also confirm this layering of charges suggesting that arginine can play a pivotal role in the nucleation of CaCO3 crystals. To elucidate the growth of ionic layers around a given arginine residue, we performed metadynamics calculations to estimate the required free energy for ion association in the higher coordination shell. The computed free energy profile shows a barrier height of ∼8 kcal mol−1 between the bound and unbound states of a CO32− ion. We also observed that the addition of a single CO32− may fetch few more Ca2+ and CO32− pairs via near-neighbor electrostatic interactions and thus increases the size of the ion cluster around arginine. In contrast, the local environment around glutamate and aspartate retain the strongly bound solvent water through hydrophilic interactions with Ca2+ that inhibit the formation of next layer of CO32− ions.
This journal is © The Royal Society of Chemistry 2019 |