1. Introduction

Minerals and proteins exist in close proximity and at nanoscale dimensions in biology. Interactions at these interfaces are vital to the functions of almost all structural materials in nature such as teeth, seashells, and bone. Knowledge of these interfaces, in particular, is useful in understanding the mechanical and physiochemical behavior of natural biocomposites such as bone, teeth, nacre, and so forth consisting of mineral and organic phases. The percentage constituent of mineral and organics is different in the various biological materials. In bones and teeth, the mineral content is about 60% and 90%, respectively [1], whereas, in nacre it varies from 95% to 98% [2]. The major component of the organic phase in nacre is proteins. Although organics are present in small proportion, they significantly alter the mechanical behavior of biocomposites. Nacre, for example, exhibits fracture toughness about 3 000 times more than the pure aragonite [2]. The proteins present in the natural biocomposite nacre show properties significantly different from any other bulk proteins. Our finite element modeling study, also later verified by experiments, shows that proteins present in nacre have modulus of about 15–20 GPa [3–7]. These simulations were done based on a parametric study of varying values of elastic modulus of the organic from 5 MPa to 100 GPa. This is about twenty times higher than the modulus of bulk proteins. The protein in nacre exists in a 20 nm space, between aragonite platelets of 200–250 nm thickness in an interlocked “brick and mortar” architecture [6]. Two possible factors which could cause the nacre protein to possess such high modulus are (i) its confinement in 20 nm space, and (ii) the mineral-protein interactions. Our current research focuses on the effect of mineral-protein interactions. The effect of confinement will be the focus of our future work. We have found in our previous work that the proximity of mineral influences the mechanical behavior of proteins [8]. We have reported that more energy is required to pull a protein at mineral proximity than in the absence of mineral. Recently, we have also reported the large influence of mineral on mechanical behavior of protein through simulations of collagen-hydroxyapatite in bone [9]. Here, we present the detailed mechanisms responsible for large changes in mechanical response during protein unfolding in the proximity of the mineral. In this molecular modelling simulation, the GS domain [10, 11] of the Lustrin A, protein which is found close to mineral surface in nacre, is pulled under two conditions: (i) in the proximity of mineral aragonite, (ii) in the absence of any mineral proximity. All other conditions of simulation remain identical for the two cases. As the GS domain is pulled keeping one end constrained, the protein domain unfolds. The mechanisms involved in the unfolding pathways are significantly different under the two conditions of pulling, resulting in large differences in energy required for unfolding at mineral proximity.

Molecular dynamics is a useful technique for studying the folding/unfolding behavior of proteins [12, 13]. The reversible unfolding of individual immunoglobulin domains was successfully investigated using a combination of steered molecular dynamics (SMD) and atomic force microscopy (AFM) [14–16]. The SMD technique has been successful in reproducing the stretching events for individual solvated I27 domain [17]. This technique has also been used for studying the mechanical properties of clay minerals at nanometer length scale [18–21]. In our previous work, [21] we have used steered molecular dynamics to understand the deformational behavior of the beta barrel and beta planar sheet structures in proteins. Also, in literature, the helix-coil transition of peptide deca-alanine in vacuum was studied using molecular dynamics simulation [22]. Researchers [23] were able to investigate the different transition states of biopeptides at different temperatures.

The model protein structure used in this study is the nacre protein Lustrin A. The primary reason for this choice is that this protein is located close to the aragonite surface and thus can potentially interact with the mineral. The protein Lustrin A is made of 1428 amino acid residues and has molecular weight of about 112 KDa. The complete structure of this protein is not known. Hence, at present, this study is confined to the response of a single-domain Glycine-Serine (GS) of the protein Lustrin A. Molecular modeling techniques such as coarse-grained modeling [24] can be applied in the future to study the response of the full protein when the structure is known. Replica-exchange molecular dynamics method has been successfully used in literature to study the folding unfolding behavior of proteins and peptides [25]. There are various experimental techniques such as laser tweezers [26], florescence quenching [27], and AFM [28], which are found to be useful in characterizing the protein folding and refolding response. Most of these techniques are fairly successful in the study of single-molecule protein folding. These studies have primarily looked at the molecular mechanisms in the unfolding of proteins, however recently we have reported measurement of adhesion forces between aragonite and nacre proteins by pulling proteins off the aragonite surface using an AFM tip [29]. Experiments were conducted on freshly cleaved samples of nacre, and protein molecules were pulled in the presence of aragonite at various velocities using an AFM tip. The results indicate that protein molecules can withstand large pulling forces of the order of 6 nN before separating from the aragonite or the AFM tip [29].

2. Model Construction and Simulation Details

One of the most widely studied nacre proteins is Lustrin A. It is a domain-based protein, and is located close to the interface formed by the inorganic (aragonite) and organics (primarily proteins). Lustrin A is populated by ten cysteine rich domains with nine protein rich domains sitting in between. It has a domain made of 275 amino acids, and significantly rich in glycine and serine. It is located close to the C-terminal and named as GS domain in this work. For the sake of these simulations, GS domain is the organic entity. The details about the possible role of each of the domains in Lustrin A are discussed in our previous work [8]. From the known primary sequence of GS-domain random, structures were generated, minimized, and then solvated by “SOLVATE” module of VMD. The aragonite is built by extending its unit cell [30] to 16 units each in - and -direction and 2 units in -direction. The solvated GS domain built earlier is placed next to the aragonite to build the organic inorganic model used in this study. The details of model building are discussed in our earlier paper [8]. It should also be noted that, since water has significant role in the influence of mineral proximity, appropriate solvation of protein is very important. In our model [8], before pulling the GS domain is completely solvated, thus enclosing the entire domain in a water box.

Steered molecular dynamics (SMD) is used to study the influence of mineral proximity on the mechanical response of protein (GS domain). One of the alpha-carbon atoms present close to aragonite surface is pulled to unfold the GS domain while fixing the other end located opposite to it. Same pulling and fixed atoms are used when GS is pulled in the presence and absence of aragonite. Three different magnitudes of velocity 0.25 Å/ps, 0.50 Å/ps, and 1.00 Å/ps are applied in unfolding GS domain. Using each of these velocities, the GS domain is pulled for a period of 250 ps. The stiffness of SMD spring used is 5 Kcal/mol/Å. Molecular dynamics software NAMD [31] is used for all simulations, and VMD [32] is used for visual analysis. CHARMm forcefield is applied in this work. The parameters for GS domain are obtained from standard CHARMm protein parameter file [13]. The CHARMm bonded parameters for aragonite are used from the literature whereas the nonbonded parameters are derived from existing Buckingham potential. The CHARMm aragonite parameters used in this work are shown in Table 1.

Table 1: CHARMm forcefield parameters for aragonite.

3. Results and Discussion

It is been found from our previous work [8, 9] that the proximity of mineral influences the mechanical response of proteins as observed in both seashells and bone. More energy is required to unfold a protein when pulled in the presence of aragonite than when pulled in its absence. It is also observed that energy required to unfold at mineral proximity depends on the velocity of pulling. The specific reason for this velocity dependence is not understood yet and is a part of our future research. In the current work, studies are conducted to find mechanisms leading to large increase in energy required for unfolding of protein at mineral proximity and in the absence of the mineral.

Protein molecules exhibit a tertiary structure due to several cross links (hydrogen bond, disulphide bridge, etc.) between domains and within domains. Specific hydrogen bonds give rise to mathematically defined structures such as alpha-helix, beta-sheets. A protein has several turns (such as hairpin) between strands. When a protein is pulled (in the absence of mineral), all these bonds offer resistance to pulling. When the protein is pulled in presence of mineral the resistance comes from both structural features of the protein molecule as well as the resistance forces due to interaction of protein with mineral.

The load-displacement characteristics of the protein domain (GS domain) pulled at mineral proximity and without the mineral at a velocity of 1.00 Å/ps are shown in Figure 1. The ratio of the area under the load-displacement (L-D) curve in the presence of aragonite to the area under L-D curve in the absence of aragonite is defined as “Work Factor.” The work factors obtained at three different velocities are shown in Table 2. From this table, it is observed that when the protein is pulled at velocity of 0.25 Å/ps, about ten times more energy is necessary to unfold the GS domain by same magnitude of displacement in the presence of aragonite than in the absence of it. It is observed from Figure 1 that the L-D responses are significantly different under the two conditions. The primary factors leading to the larger area under the L-D curve pulled at mineral proximity are

Table 2: “Work factor” at different velocities.

Figure 1: The L-D characteristics of the GS domain, pulled at a velocity of 1.00 Å/ps: at proximity of aragonite and in the absence of aragonite.

(a)the presence of peaks of higher magnitude in L-D response curve for mineral proximity;(b)and the presence of higher value of load for a given displacement. In the following two sections, we present the results of our study about the mechanisms arising from the above two factors.

3.1. Section I

Peaks in the L-D curve are present in both plots as indicated in Figure 1. The peaks observed in the absence of mineral aragonite have a height of about 150 to 200 pN, whereas the peaks in the presence of aragonite are of larger magnitude and are in the range of 1500 to 2000 pN and have lower frequency of occurrence. The magnitude of peak heights in the absence of aragonite corresponds to the opening of loops and turns and breaking of hydrogen bonds. From the nature of peaks observed in two cases, it can be inferred that different mechanisms are involved in the formation of peaks when pulling protein at mineral proximity. From the trajectory analysis of unfolding, it is observed that peaks occur until water enters between the aragonite surface and GS strands. As pulling is continued and water enters in between, no additional peaks are observed in the L-D curve. This observation indicates that direct interaction between aragonite surface and GS strands may result in the formation of these peaks. We have investigated the peaks up to 50displacement. This displacement corresponds to the point after which significant water enters in between the aragonite surface and GS strand. From the trajectory analysis, it is observed that when the GS strand is pulled at mineral proximity, a mechanism similar to the physical phenomenon of “latching” is observed to take place between aragonite surface and GS strand in proximity. The “latch” observed here is not a physical latch but a result of strong nonbonded interactions between an atom or a group of atoms of amino acid residues (attached to GS strand) and aragonite surface. These latches are formed when the interacting group(s) come(s) close to the aragonite surface, but as pulling is continued they are observed to break. These latches thus go through a process of “forming" and “breaking" as the GS strand is pulled along the aragonite surface, giving rise to the peaks in the L-D curve.

Further, we zoomed into a small section, circled in Figure 2(a) where, a “latch” is observed. The “forming" and “breaking" phenomenon of a latch is represented in Figures 2(b) and 2(c). The aragonite layer close to the surface is shown in the figures. The L-D characteristics of one of the peaks are shown in Figure 2(d). As the strand is pulled in the direction of the arrow in the figure, a group of atoms (marked as blue balls) comes closer to the surface, builds a strong interaction, arrives between the carbonate groups of aragonite, and gives rise to what we called as a “latch." The formation of the “latch” corresponds to the base point of the peak as marked in Figures 2(b) and 2(d). Once the latch is formed, it offers resistance to pulling, causing the magnitude of load to rise sharply. The load curve keeps rising as pulling is continued until the load level where the latch can no longer sustain the pulling force results. At that point, the latch breaks, and thus the force decreases. The apex of a peak therefore corresponds to the point where a latch breaks (Figure 2(c)).

Figure 2: The mechanism of “latching” and formation of peaks: (a) latching site (b) formation of latch (c) breaking of latch (d) peaks resulting from latching mechanism.

Aragonite exhibits orthorhombic crystal symmetry. The surface formed by (001) is rich in negatively charged oxygen atoms belonging to carbonate group. The positively charged calcium atoms are located very close to the surface. The charges on the surface oxygen and calcium atoms are found to be −0.9995 and +2.000, respectively. The blue atom group constitutes mainly of hydrogen and/or oxygen atoms. The hydrogen and oxygen atoms of amino acid strongly interact with oxygen and calcium of aragonite, respectively, forming these latches. As polar water molecules enter between the aragonite surface and GS strand, these interacting groups can no longer come close enough to form the latches and thus no sharp peaks are observed. In Figure 1, no peaks are observed beyond a displacement of about 60 Å. The peaks within the same L-D curve (at mineral proximity) are different in terms of height, base width, shoulder characteristics, etc. This heterogeneity in peak characteristics depends on the number of latches formed at a given time, the atoms involved (blue group) in the “latching” action and on the presence of other interacting groups closer to the aragonite surface.

Thus, a peak in the L-D curve starts rising when a side chain group builds up a strong nonbonded interaction with aragonite surface. It falls when this interaction is broken upon further pulling. Thus, the energy stored in the rising half of the peak is due to the resistance offered by the latch and is thus expected to be approximately equal to the nonbonded energy holding these two interacting groups together. The nonbonded energy which is the summation of van der Waals and electrostatic energy between the interacting group of amino acid side chains, and the aragonite section (shown in Figure 2(b)) is determined. The variation of the electrostatic and van der Waals energy with time (frame) is shown in Figures 3(a) and 3(b). In Figure 2(d), it is observed that Peak 0 starts rising at time frame 85 (17 ps). The electrostatic and van der Waals energy corresponding to this time frame as observed in Figures 3(a) and 3(b) is 2000 pN-and 700 pN-Å, respectively. The area bound by the rising half of the Peak 0 (, and in Figure 2(d)) is 2500 pN-Å, which is close to the total nonbonded interaction (2700 pN-Å) energy between the interacting groups. From a similar analysis of Peak 1 (Figure 4(a)), where the peak starts at time frame of 135 (27 ps), the total attraction energy (900 pN-Å) is found to be very close to the area bound by the rising part of the curve (1100 pN-Å). The variation of the electrostatic and van der Waals energy with time (frame) for Peak 1 is shown in Figures 4(b) and 4(c).For other major peaks, these values are also in close agreement. Thus it is quantitatively shown that the large peaks observed in the proximity of the mineral are a result of direct nonbonded interactions between the side chain atoms and the complementary atoms on the mineral surface.

Figure 3: The nonbonded energy between the interacting groups giving rise to Peak 0: (a) electrostatic energy (b) van der Waals energy.

Figure 4: Peak 1 (a) load-displacement plot (b) electrostatic (c) van der Waals energy.

3.2. Section II

As observed in Figure 1 and also mentioned earlier, one of the factors which contribute to the work factor is the higher magnitude of load at any given displacement, when pulled at mineral proximity. The peak load observed in the figure is about 13 nN, a high magnitude, predominantly as a result of mineral protein interactions. This high magnitude of force resulting from constant velocity pulling simulations is possible without fracture of the single protein strand being pulled. Flow induced fracture of single polymer chains is reported to be between 2.5 to 13.4 nN [33, 34] and the fracture strength is dependent upon the polymer and rate of loading. In this section, the potential mechanisms involved are described. As the GS domain is pulled, the water molecules present closer to the aragonite surface orient and arrange themselves in a definite pattern and further undergo very little movement. On moving away from aragonite surface, the water molecules are structurally altered, and behave like “normal solvent water" attached to the protein.

A clear difference in behavior is observed between water molecules closer to the surface and away from the surface both in terms of structure and movement during pulling. We have therefore categorized the water into two types. The water molecules close to the aragonite surface are named here as Aragonite-Bound Water (A-BW), and water molecules moving with the protein are named as Protein-Bound Water (P-BW). The schematic representation of these two types of water molecules is shown in Figure 5. At low displacement of pulling, very few water molecules are observed to be present between aragonite and GS strands. However, this number increases as pulling is continued. The strong hydrogen-bonded interaction between carbonate surface of aragonite and polar water molecules causes some of the water molecules to lie close to the aragonite surface. The water molecules form a periodic pattern which is primarily due to the periodicity of the aragonite surface in the - plane. The P-BW shows no definite pattern, with the water molecules surrounding the protein and moving freely with the protein. As the GS is pulled, the A-BW remains tightly bound to the aragonite surface and exhibits almost no movement. The P-BW on the other hand moves freely with the GS strand. On careful observation, it is observed that there is a 3-4thick layer of water (about two water molecules thick) which behaves as a transition layer between A-BW and P-BW. In this layer, the slipping of P-BW is observed to happen with respect to the A-BW. The P-BW forms a sheath around the GS strands. When GS is pulled, the strands move with sheath of water around it. P-BW thus moves with GS strands and A-BW remains bound to the surface. A schematic representation of this mechanism is shown in Figure 6. This figure illustrates that when GS is pulled, the interaction contributing to resistance to displacement includes the following:

Figure 5: Two “types” of water molecules observed during pulling of GS domain (a) simulation (b) schematic representation. A-BW: aragonite bound water and P-BW: protein bound water.

Figure 6: Schematic representation of movement of the GS domain and protein-bound-water (P-BW) when GS domain is pulled.

(i)the net attractive interactions of GS-(A-BW) and GS-aragonite;(ii)the net attractive interaction of (P-BW)-(A-BW) and (P-BW)-aragonite and GS-(P-BW).

To study the contribution of these different interactions to the L-D characteristics of GS domain, aragonite sections of smaller - dimension and part of GS strand close to it are selected. We have used two such sections for the analysis of the interaction energy. One of them is shown in Figure 7. This section extends to 20in the -direction and 50in the -direction. Figures 7(a) and 7(b) show this section at 180 ps and 250 ps, timeframes, respectively. The white region between A-BW and P-BW, in reality, contains water molecules, however these are not shown here for clarity. This region is the layer of transition or slippage as mentioned earlier. Figure 7(b) shows the movement of P-BW (green atoms) as the strand is pulled. The interaction energy for each of the pairs is shown in Table 3. The electrostatic interaction energy is dominant compared to the van der Waals energy. The attractive nonbonded interaction energy between mineral aragonite and A-BW is pN-Å. Thus, the A-BW is tightly bound to the aragonite surface and does not move when the GS domain is pulled alongside the mineral.

Table 3: Values of interaction energies between different pairs of constituents of aragonite-water-GS domain.

Figure 7: Snapshots taken during the steered molecular dynamics simulation of the pulling of the GS domain in the presence of aragonite mineral proximity at (a) 180 picoseconds and (b) 250 picoseconds.

The attractive interaction energy between aragonite and P-BW is pN-Å, a very significant magnitude. The P-BW and the GS domain have an attractive interaction energy of pN-Å. In addition, aragonite and A-BW also have attractive interactions with the GS domain, although smaller magnitudes of pN-and pN-Å, respectively. The dominant resistance to pulling of the GS is the attractive interaction between GS domain and the P-BW, which in turn has significant interaction with the mineral. Although the direct interaction between the mineral and GS is small because of the distance, water plays a critical role in building a “bridge” between the mineral and the protein, facilitating development of a large resistance to pulling in the presence of mineral. As a result, the load experienced at the pulling end is high. This therefore results in higher magnitude of load at a given displacement in the L-D curve.

4. Conclusions

In our previous work, we have shown that the proximity of mineral influences the mechanical response of protein. More energy is required to unfold a protein in presence of aragonite than without it. The L-D curve in the presence of aragonite shows different responses from the one when pulled without mineral proximity. In the presence of aragonite, larger peaks are observed, and the magnitude of load at any given displacement is significantly higher. In this work, we have found quantitatively the mechanisms leading to the difference in load displacement response of protein at mineral proximity. The following hold.

(i)At the early stage of pulling, the peaks in the L-D plot at mineral proximity are quantitatively correlated to the interaction energy between the atoms involved in the latching phenomenon of amino acid side chain to aragonite surface.(ii)The role of water in mineral and protein interaction is very significant.(iii)The water closer to the mineral surface is highly oriented and does not move while the protein strand is pulled. Water layer around the strands moves with the strand as the protein is pulled.(iv)The high magnitude of load for a given displacement originates from attractive interactions between the protein, protein-bound water, and the mineral. Here, for the first time, quantitative description of the mechanics responsible for the large differences in the magnitude of force needed to unfold proteins in the proximity of the mineral in nacre is provided. This work provides clues as to possible reasons for extraordinary properties of organic phase such as high elastic modulus and large deformation before failure observed in nacre.

Acknowledgments

Teragrid allocation (TGDMR060001T) and NCSA supercomputing resources were used for majority of the simulations. NDSU center of supercomputing applications (CHPC) and Biomedical Research Infrastructure Network (BRIN) were used for model development. Hardware support for NAMD at NDSU was provided by Dr. Gregory Wettstein. P. Ghosh acknowledges support from ND EPSCoR for doctoral dissertation award at NDSU.