Discontinuous Force Compression Curve of Single Bovine Carbonic Anhydrase Molecule Originated from Atomistic Slip

The compression process of a single bovine carbonic anhydrase II ( BCA II ) molecule by an atomic force microscope (AFM) tip in ultra high vacuum (UHV) condition is studied by the molecular dynamics (MD) simulations with the all-atom empirical force field model. The temperature is assumed to be 0 and 300 K, and the force curves are calculated with both the tip and surface modeled as rigid planar graphite sheets. At T = 0 K the normal force of the tip increases gently after the tip starts to compress the protein molecule, as the external force by the tip is consumed to relax the turn or coil structures on the outer part of the protein. However, once this region is deformed to take a flat structure, the normal force increases rapidly. The general trend in the force curves at T = 300 K is similar to that observed at T = 0 K, although the relaxation inside the protein is promoted thermally to reduce the normal force observed. Furthermore, the saw-tooth peaks observed in the force curves are found to originate from the sudden collective sliding motion of the residues, which is triggered by the slip of the locally contacting turn or coil structures at the atomic level. [DOI: 10.1380/ejssnt.2006.552]


I. INTRODUCTION
In the research field of single molecule force spectroscopy of biomolecules, the atomic force microscope (AFM) tip has been used to investigate their mechanical properties such as force extension curves [1][2][3][4][5][6][7][8][9]. In these experiments, one end of the target biomolecule is tied or immobilized to the substrate, while the other end connected to the AFM tip is moved apart from the substrate. The force extension curves obtained in this way have shown highly nonlinear features and saw-tooth patterns in many cases, which reveals the heterogeneous interactions inside the biomolecules. On the other hand, the measurements of force compression curves, which give information on hardness like Young's modulus, have been performed only for a few biomolecules [10][11][12].
Correspondingly, theoretical simulations for the force extension curves have been performed for several biomolecules [13][14][15][16] by adopting the scheme called the steered molecular dynamics (SMD) [17,18]. In this method, the virtual harmonic springs are attached to the two ends of the target biomolecule and they are moved apart in the opposite directions to mimic the movement of the AFM tip in experiment. These calculations have contributed to understanding of the unfolding process of biomolecules at the atomic level. However, atomistic simulations of compression of biomolecules have not been performed yet, and it has not been totally understood what occurs under the AFM tip and what is actually measured in experiment.
Our purpose in this work is to develop a simulation * Corresponding author: tagami@cms.nano.waseda.ac.jp method of the compression process of biomolecules under an AFM tip. In this article, we choose a case study of a bovine carbonic anhydrase II ( BCA II ) molecule adsorbed on the graphite surface in ultra high vacuum ( UHV ) condition, and examine the developed simulation method. The carbonic anhydrase (CA) works as a catalyst to promote the reversible hydration reaction, CO 2 + H 2 O = HCO − 3 + H + , with a high efficiency. This enzyme is fundamental in many physiological processes such as respiration, renal tubular acidification, bone resorption, and so on [19]. Among its many isozymes, CA II has been known to be the most catalytically active, and to be the most widely distributed in tissues. In the case of human, in particular, its deficiency causes diseases such as osteopetrosis [20,21], renal tubular acidosis [22], and so on. On the other hand, these diseases have been found to have a relationship to the folding patterns of CA II [23]. The BCA II is a homologue of human CA II, and has been also studied from the viewpoint of the protein folding process [24,25].
The substrate surface is assumed to be the highly oriented pyrolytic graphite (HOPG) surface, and is modeled by a graphite sheet of mono-layer thickness. The tip is also modeled by a graphite sheet to mimic the bottom of the microbead which has an extremely large radius of curvature and is frequently adopted in the AFM measurement of biomaterials. ( In experiments, the microbead is attached to the cantilever and in the vicinity of the tip, and thus it is different from the ordinal tip. However, in the viewpoint of the compression tool of materials, they play similar roles and their difference is only the radius of curvature. Thus, in this article, we use the word tip in a wide meaning instead of using the word microbead. ) The BCA II has a kind of multi-layered structure, i.e., the central core of a β-sheet surrounded by α-helix, turn, and coil structures. The dimensions of the molecule are approximately 46 × 51 × 51Å 3 . At T = 0 K the force curves are calculated by the energy minimization method, while at T = 300 K they are calculated by the molecular dynamics simulation method. The whole calculations are performed using the program called NAMD 2.5 [26] with the all-atom empirical force field called CHARMM 22. We found the following characteristic features by the simulations; At T = 0 K the normal force of the tip increases gently after the tip starts to compress the protein, since the mechanical energy injected by the tip is consumed to relax the turn or coil structures in the outer part of the protein. However, after this region is deformed to take a flat structure, the normal force increases rapidly. The general trend in the force curves at T = 300 K is similar to that observed at T = 0 K, although the thermal motion promotes the relaxation inside the protein and reduces the normal force. Furthermore, the saw-tooth peaks observed in the force curves are found to originate from the sudden collective sliding motion of the residues, which is triggered by the slip of the locally contacting turn or coil structures at the atomic level.
The rest of the paper is organized as follows. In §II the computational model and simulation method will be described. In §III the results of our calculations will be presented and the mechanism of the jump in the force curves will be analyzed at the atomic level. In §IV the paper will be summarized.

II. MODEL AND METHOD
The initial structure of the BCA II molecule is prepared from 1V9E.pdb by removing the extra water molecules and Zn 2+ ion, then by adding the missing hydrogen atoms to each of its amino acid residues [27]. The effect of the existence of Zn 2+ ion on the compression process will be studied in future work. The molecule is composed of 259 amino acid residues, and takes a sphere-like structure whose dimensions are about 46 × 51 × 51Å 3 . The substrate surface is modeled as a mono-layer graphite sheet whose lateral dimensions are 101.0 × 101.8Å 2 . The molecule is deposited on the surface so that its C-terminal residue faces the surface. The AFM tip is modeled by the monolayer graphite sheet whose size is the same as that of the substrate surface, and is set to be parallel to the substrate surface. The number of atoms in each component is 4,045 (protein), 4,032 (surface), and 4,032 (tip), respectively, which amounts to 12,109 atoms in total. The interatomic interaction is described by the all-atom empirical force field called CHARMM 22. For the graphite sheets, we adopt the potential parameter for the atom type "CA". The van der Waals interactions are cut off beyond 12Å. The system is negatively charged by -3 e, and no boundary conditions are imposed. The whole calculations in the followings are performed using the program NAMD 2.5 [26]. The relaxation of the surface and tip is not considered in this work. Figure 1 shows the structure of the BCA II molecule adsorbed on the graphite surface, which is obtained by the molecular dynamics (MD) simulation for 4.2 ns with the time step of 2 fs at T = 300 K. The N-and C-ter- The residues at N-and C-terminus are illustrated by the van der Waals representation with yellow and cyan color, respectively. The α-helices and β-strands are labelled by H1-H7 and S0-S12, respectively. The width (x-direction) and depth (y-direction) of the BCA II are 46.2 and 47.0Å, respectively.
minal residues ( Ser1 and Lys259 ) are drawn by the van der Waals representation with yellow and cyan colors, respectively. The residues with red and green/yellow colors correspond to the α-helix and β-strand structures, respectively. The green colored β-strands survive at the smallest tip-surface distance h = 30Å at T = 0 K, as will be presented later, while the yellow colored β-strands do not survive. The blue colored residues correspond to the turn or coil structures which take the α-helix structures in the crystal phase, but change their forms by thermal motion and adsorption on the graphite surface. These α-helices and β-strands are labeled by H1 to H7 and S0 to S12, respectively, whose definitions are listed in Table I. Hereafter, the axes with red, green, and blue colors are defined to be x, y, and z axis, respectively. Initially, the tip height is set to h = 60Å measured from the substrate surface, and is lowered in a stepwise manner. At each tip height, the normal force of the tip is calculated at the two system temperatures, 0 and 300 K. At both temperatures, only the protein molecules are allowed to relax. To be specific, in the simulation at T = 0 K, we adopt the energy minimization method. Namely, once the tip height is determined, the protein structure sandwiched between the substrate and tip is optimized for 50,000 steps with the conjugate gradient method. After this, the normal force of the tip is calculated by summing the atomic forces over the whole tip atoms. Then, all I: Amino acid residues in the secondary structure defined in Fig. 1. In the simulation at T = 300 K, we adopt the molecular dynamics method. The temperature is controlled by a LANGEVIN damping parameter of γ = 5 ps −1 . The time step is set to 2 fs and the bond length constraint is applied to hydrogen atoms. The atomic forces are monitored over 200,000 steps ( 0.4 ns ). The mean normal force is calculated by summing these force components over the whole tip atoms at each MD step, and then by averaging over these 200,000 MD steps. Next, the tip height is lowered by 0.5Å with its structure unchanged. These procedures are performed repeatedly. Since the thermal equilibration is not checked in this scheme, it may be called the quasi-equilibrium MD. We also present the results by the equilibrium MD in which, prior to monitoring the atomic forces, the MD simulation is performed for 2-4 ns at every tip height and the convergence of the root mean square displacement (RMSD) as well as the system energy are checked.

III. RESULTS AND DISCUSSION
A. General trend: two-step compression Figure 2 shows the calculated normal force at T = 0 K plotted against the tip-surface distance h. The general trend is described as in the following. First, an attrac- tive normal force sets in at h ∼ 54Å, and its maximum appears at h ∼ 50.2Å. It turns to be repulsive at h ∼ 48Å, and increases almost linearly until h = 43Å. Then it takes an oscillatory behavior with a succession of sawtooth peaks in the region of 35 < h < 43Å. After these, the slope of the force curves gets steeper, and finally exceeds 12 nN at h = 30Å. α-helix H4 are large. Noticeably the upper part of the protein shows an almost flat structure. At h = 35Å the lower part of the outer protein structure is found to be further compressed, which results in lowering the heights of the β-strands (S2, S3, S8) and α-helices (H4, H5, H7). In addition, the β-strand S1 changes into a coil structure. On the other hand, the residues Asp31 and Ile32 newly form a parallel β-sheet with Thr107 and Val108, although it is not clearly seen in Fig. 3(c). Finally, at h = 30Å, both the upper and lower parts of the outer structure are found to be compressed further. The β-strands at the left and right sides ( S0, S2, S3, S8, S12 ) change into the coil or turn structures, and the three α-helices denoted above decrease their heights. In this way, once the protein makes a close contact with the tip by making its upper part of the outer structure flat, it starts to become stiffer suddenly, which shows up as a rapid increase in the normal force at h < 40Å ( see Fig. 2 ). Figure 4 shows the force curves calculated at T = 300 K. The solid and broken curves correspond to the results obtained by the quasi-equilibrium and equilibrium MD simulations, respectively. In the quasi-equilibrium MD, first, the normal force takes attractive values whose maximum appears at h ∼ 49.5Å. It turns to be repulsive at h ∼ 46.8 A, and increases gradually until h = 40Å. Then it drops to almost zero at h ∼ 38.5Å, and increases again with a steeper slope. Finally, the normal force takes 5.5 nN at h = 30Å. Basically, the force curve in the equilibrium MD shows a similar trend to that in the quasi-equilibrium MD, although the locations of the local maxima or minima differ slightly. In the followings, however, we focus on the results obtained with the latter method, as in the former method the protein moves from the center of the substrate surface to the vicinity of its edge during the long simulation time for equilibration repeated at each tip height. This might give an artificial contribution to the values of the normal force as the periodic boundary conditions are not applied. Figure 5 shows snapshots at T = 300 K at the four , which are taken at the last MD step at each tip height. As a general trend, the feature of structural deformation is similar to that found at T = 0 K ( see Fig.  3 ). Namely, first, the upper outer part of the BCA II molecule is compressed, then the lower outer part is also compressed. At h = 40Å, the upper part of the protein becomes almost flat, which indicates that the contact with the tip is formed not at a single point but by a planar area, just as in the same way to the results at T = 0 K. It shows up as a rapid increase of the normal force at h < 39Å. However, the three remarkable differences can be observed between the snapshots at T = 0 and 300 K. At T = 300 K, (1) the α-helix H4 is pushed up and becomes close to the tip at h = 30Å. (2) The heights of the α-helices H5 and H7 are lowered. In particular, the axis of H5 is inclined to be almost parallel to the surface at h = 30Å. (3) The heights of the β-strands S5 and S6 are lowered. These features are commonly observed in the equilibrium MD simulation. Such structural relaxation inside the protein promoted by thermal motion reduces the normal force of the tip which takes 5.5 nN at h = 30Å at T = 300 K (  Fig. 4 ), lower than 12.0 nN at T = 0 K ( Fig. 2 ).

B. Origin of saw-tooth peaks
As denoted in the previous subsection, the force curve at T = 0 K shows saw-tooth peaks at several tip heights. The four saw-tooth peaks with largest top-to-bottom heights are marked with arrows in Fig. 2. Although even at T = 300K a similar feature is observed, it is difficult to distinguish the purely mechanical effects from the thermal fluctuation at this stage. Thus, in this subsection, we will discuss what occurs at such saw-tooth peaks observed at T = 0 K. Hereafter, the atomic displacement at a certain saw-tooth peak is defined to be the atomic displacement from the top to the bottom of a selected peak. At the peak a in Fig. 2, it turns out that the coils and turns at the upper part of the BCA II, including H1 and H2, slide along the lateral direction. The maximum displacement is observed to be 5.49Å at His14 which belongs to the turn structure H1. In addition, this seems to be triggered by the structural change at Tyr6, Asn242, and Trp243. The locations of the three residues are illustrated with the van der Waals representation in Fig.  6 (a), which are viewed from the [100] direction. Here, Tyr6 lies below His14, and Asn242 and Trp243 form a part of a coil structure. Figures 6 (b) and 6 (c) show the closeup of these residues at the top and bottom of the saw-tooth peak, respectively. As clearly seen, Tyr6 slips on Asn242 and Trp243, by which the stress imposed by the tip is released. It shows up as a sudden jump in the force compression curve.
At the saw-tooth peak marked b in Fig. 2, the maximum displacement of 3.22Å occurs at Leu187 ( not shown ) which belongs to a coil structure. As a general feature, the left side of the protein ( see Fig. 3(b) ) moves downward while the right side moves upward, by which the height difference on the top surface of the molecule is reduced. At the intermediate region which is sandwiched by the above two domains moving in the opposite directions, a structural change like a dislocation is observed. For example, Phe19 moves downward while Phe129 moves upward ( not shown ).
At the saw-tooth peak c in Fig. 2, the lower part of the protein is found to be stretched along the y axis. The maximum displacement of 6.29Å is found at Tyr50 which belongs to a coil structure. The phenyl group in this residue abruptly rotates to face the graphite substrate surface, synchronized with the rotation of the phenyl group of Phe177 which lies just above it. Figure 7(a) shows the relevant residues located nearby. At the same time with the structural change as denoted above, the OH group of Phe177 moves downward, which pushes out Lys157 toward the [010] direction. On the other hand, Ala76 slips on Gly51 ( see Figs. 7(b) and 7(c) ) which is the adjacent residue to Tyr50, and comes close to the graphite surface. This brings about the sliding of the turn and coil structure which contains Ala76 toward the [010] direction. In this way, the BCA II gets stretched along the y axis.
Finally, we comment on the saw-tooth peak marked d in Fig. 2. The maximum displacement is found at the residue Ile21 with 4.64Å. As a general feature, 1) the α-helices H5 and H7 do not move significantly while the α-helix H4 moves downward, and 2) the turn and coil structures on the upper part moves along the lateral directions. Figure 8(a) shows the three residues relevant to such a structural change occurring in a collective manner. Thr197 and Pro12 belong to different turn structures while Val120 belongs to the β-strand S6. From the top to bottom of the saw-tooth peak, these residues change their positions as illustrated in Fig. 8(b) to Fig. 8(c). Namely, Pro12 moves toward the [010] direction while Thr197 moves toward the [010] direction. Accompanied with this, the turn structure which contains Thr197 and Leu202 is pushed out toward the [010] direction. Then, the latter residue touches and pushes down Ala133 which lies at the end of the α-helix H4.
In this way, we find that the saw-tooth peaks have close relationship with the collective sliding motion of the residues which is triggered by the slip of the turn or coil structures at the atomic level. We should point out here that similar saw-tooth peaks have been obtained in the elongation of the BCA II molecule [16] and other biomolecules. However, the origin of such peaks has been attributed to the breakdown of the hydrogen bonds between the residues, typically formed between the βstrands, and is considered to be different from the mechanism presented in this article. We expect that the sawtooth peaks based on the atomistic slips will be observed in the compression measurement of the BCA II molecule as well as the other protein molecules.

IV. SUMMARY
The compression process of a single protein molecule ( BCA II ) by an AFM tip is studied in the UHV condition by the molecular dynamics (MD) simulations based on the all-atom empirical force field model. The temper-ature is assumed to be 0 and 300 K, and the force curves are calculated with both the tip and surface modeled as rigid planar graphite sheets. At T = 0 K the normal force of the tip increases gently after the tip starts to compress the protein molecule, as the external force by the tip is consumed to relax the turn or coil structures on the outer part of the protein. However, once this region is deformed to take a flat structure, the normal force increases rapidly. The general trend in the force curves at T = 300 K is similar to that observed at T = 0 K, although the thermal motion promotes the relaxation inside the protein and reduces the normal force observed. Furthermore, the saw-tooth peaks observed in the force curves are found to originate from the sudden collective sliding motion of the residues, which is triggered by the slip of the locally contacting turn or coil structures at the atomic level.
However, the following problems remain to be solved. The surrounding water molecules and Zn 2+ ion are significant for the activity of this protein molecule in its native condition. The effects of these factors on the mechanical deforming process should be clarified by further calculations. Furthermore, we do not achieve the Young's modulus of the BCA II from a limited number of force curves obtained in this work, because such a quantity should be evaluated through a statistical average over a large number of force curves. In order to obtain the Young's modulus, more simulations should be required by changing the initial conditions, orientations of the molecule, and so on.