{'fraction_non_alphanumeric': 0.03118020958598179, 'fraction_numerical': 0.023621370898471054, 'mean_word_length': 4.979969183359014, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Molecular Dynamics Study to Assess the Positive Ion Distribution and the Effects of Protonation on the N-terminus Region of a Xylanase

Xylanase enzyme families play an important role in many industries due to its ability to digest the hemicellulose within plant cell walls. Developing the enzyme that can withstand a highly acidic environment in a stomach requires an extended information in molecular scale for further protein engineering. In this study, an atomistic molecular dynamics simulation was performed to assess the positive ion distribution pattern on the N-terminus region of a Xylanase enzyme molecule. Then, another series of simulations were performed for the enzyme protonated at the aspartate residues 11 and/or 32 identified by the previous calculation. Conformational analysis showed that ASP32 protonation caused a significant disruption at the N-terminus region, as the nearby arginine residue 31 (ARG31) became free from the Coulomb interaction with the negatively charged ASP32 sidechain and cause the N-terminus region to misfold.


Introduction
In modern animal feedstock industry, enzyme technologies are used during the processing of food ingredients [1]. One of the most popular enzyme family is Xylanase. Xylanase possess its function on hemicellulose digestion, which is an important step prior to the total cellulose digestion for the plant cell walls [2]. To implement xylanase into the feedstock, the enzyme is mixed with the other ingredients for the primary digestion before undergoing a granulation process by an extruding machine. However, to further improve the efficiency of the enzymes, secondary digestion of the enzymes within the stomachs has been proposed. To reach this goal, enzyme molecules with high resistivity to high temperature of the extruding machine and to low pH in the stomachs are needed [3]. A number of enzyme engineering attempts have been done both experimentally and computationally to improve thermostability of xylanase and other families of cellulase and hemicellulase enzymes [4]. Many studies on the molecular biology and molecular biophysics of xylanase enzymes from the GH11 family (Xyn11A) refer the threedimensional structures of an Xyn11A molecule as a 'jellyroll' or a 'right hand' (see figure 1.) [5,6]. Atomistic molecular dynamics simulations, a physics-based computational approach, have long been used in the studies of protein structures and functions due to their abilities to provide atomistic details of both equilibrium and non-equilibrium processes occurs within protein structures under some extreme environments. Our recent work on improving the thermostability of a xylanase enzyme of the from Bacillus sp. showed that introducing a pair of disulfide bond connecting the helix and the 'back-hand' beta hairpin-loop greatly improve the thermostability [7]. However, another weak-spot on the protein structure was found at the 'Pinky-finger' region at the N-terminus. Thus, one of the ongoing studies is focusing on the addition of a hydrophobic residue to stabilize the beta-strands around the 'Pinky-finger' region.
In this work, the distribution of positive ions and the stability of protonated 'Pinky-finger' region were investigated. In the first part, an atomistic molecular dynamics simulation of a wildtype Xyn11A recombinase from Bacillus sp. was performed to see the distribution of positive ions about each residue at the 'Pinky-finger' region (amino acid residues 1-50). High salt concentration was chosen in order to accelerate the sufficient statistics on the ion distribution. After that, two residues were chosen to create protonated enzyme structures for further series of MD simulation runs. Conformational analysis addressed the molecular mechanisms of the local structural disruptions due to the protonation, which could serve as information for further protein engineering towards the acid-stable Xyn11A. The 'Pinky finger' region was indicated in Pink, while the beta strands B1, A1, A2 and B2 are labelled, and residues ASP11 and ASP32 are represented by balls and sticks.

Methods
Atomistic structure of Xyn11A was obtained from a homology modelling and were parameterized by the AMBER03 force-field [8]. The protein was solvated in TIP3P explicit solvent [9] (66.5×66.5×66.5)Å 3 water box, with Na+ or Cl-counter-ions to neutralize the charged amino acid residues. Additional Na+ and Cl-ions were added to emulate 0.5 M salt concentration. Energy minimization using the steepest descent method and 0.5 ns MD equilibration, where each structure was gradually heated from 100K to 333K, was performed. After that, 100 ns production MD run was carried out by GROMACS 5.1.2 at 1 bar and at 333 K [10]. Berendsen and velocity rescale algorithms were employed to regulate pressure and temperature in the NPT ensemble [11]. The particle mesh Ewald (PME) method [12] and a cutoff distance of 10 Å were used. LINCS constraint algorithm [13] was applied to covalent bonds associated with hydrogen, so that the timestep of 2.0 fs was allowed. After simulations finished, all the trajectory snapshots were fitted to their starting structures by using the least square fitting method from the GROMACS module trjconv. Then, to observe the ion distribution and determine the sites for protonation, minimum distances between Na+ ions and Pinky-finger region were analyzed. After the protonation sites were decided, protonated starting structures were built and followed similar simulation protocol as the non-protonated structure. Root mean square deviation (RMSD) of each Xyn11A structure, excluding 'Thumb' region, relative to the starting structure was then calculated for each 10 ps of the simulation trajectories at all protonated conditions to measure the ensemble average of local conformational changes along a 100-ns trajectory. RMSD at a given time t is defined by while mi and ri are the mass and the position of atom i, M is the total mass and N is the total number of atoms. Also, the Pinky-finger region (residue 1-50) was monitored. Then, root mean square fluctuation (RMSF) per residue was calculated to measure the local fluctuation affected by amino acid protonation. RMSF of an atom i is defined by while ri is the position of atom i and T is the total period of time, Visual molecular dynamics (VMD) program [14] was used to visualize the molecular structures.

MD Simulations to Study Ion Distribution Around a Xylanase Molecule
In order to determine the positive ion binding sites within the Pinky-finger region, probability of Na+ ions with minimum distances from amino acids less than 4 Å was analyzed from the MD trajectory of Xyn11A wildtype structure at 333K. From the probability of ion binding at each residue shown in figure  2, the two most ion-populated residues were ASP11 and ASP32 (indicated by red arrows). It could be explained that aspartate sidechains had low pKa value and was negatively charged at neutral pH environment. In an acidic environment in which protons were abundant, aspartate sidechains were more likely to be protonated. Therefore, ASP11 and ASP32 were chosen as the protonated residues in the latter simulations in order to monitor the conformational changes and structural disruptions due to the amino acid protonation.

Figure 2.
Probability of Na+ ions to be in close proximity (minimum distance d < 4 Å) with each amino acid. Probability at the residues ASP11 and ASP32 were indicated by arrows. Figure 3a and 3b showed the Root Mean Square Deviation (RMSD) values calculated from the MD trajectories of Xyn11A structures with different protonation states . RMSD was firstly calculated for the whole protein structure excluding the thumb region. Figure 3a showed that the RMSD at the 'nothumb' region was found the highest for the structure with ASP32 protonation after 90 ns MD simulation (4.2 ± 0.1 Å). The rest of the structures had lower RMSD values (3.0 ± 0.2 Å) than the ASP32-protonated structure. Similar calculations were done for the 'Pinky-finger' region and the ASP32-protonated structure also displayed the highest RMSD (4.4 ± 0.2 Å), followed by the ASP11/ASP32-protonated structure (3.4 ± 0.2 Å). These results were surprising because the structure with a single protonation at ASP32 underwent significantly larger conformational changes than that with double protonation. To investigate the local flexibility at each amino acid residue, Root Mean Square Fluctuation (RMSF) was calculated for all four structures. The dashed lines in figure 3c indicated the positions of ASP11 and ASP32 protonation. The RMSF results showed that, for all three protonated states, there were significant increases in the fluctuations of B1/A1 loops, including the ASP11 residue, compared to the nonprotonated state. However, the ASP32 residue was unaffected even when ASP32 itself was protonated.   Figure 4a shows the conformation of the Pinky-finger region (residue 1-50) from the wildtype, non-protonated Xyn11A. A hydrogen bonding network formed between two hairpin-loops, B1/A1 (residues 11-15) and A2/B2 (residues 31-35), was similar to that found between beta-strands, but with some mispairing due to the flexibility of the loop. Also, additional hydrogen bond were found between one of the terminal oxygen atoms of ASP11 and a secondary amine (N-H) group at the backbone of nearby residue. ASP11 protonation, shown in figure 4b, reduced the probability of h-bonding between ASP11 and the backbone, resulted in the increased flexibility of the B1/A1 loop, consistent with RMSF result. The most significant changes of RMSD were observed when Xyn11A was protonated at the ASP32 residue and the charge of the protonated-ASP32 sidechain became neutral. As a result, Coulomb interaction between ARG31 and ASP32 sidechains was removed. ARG31 then created a non-native hydrogen bonding network, which exposed the B1/A1 loop to further misfolding created by the extended hydrogen bonding network with ASP11. Meanwhile, the ASP11/ASP32 protonated structure shown in figure 4d was less disrupted, as one of the two terminal oxygen atoms became protonated and was less likely to form the extended hbonding network seen in the ASP32 protonated structure. (c) ASP32 (d) ASP11/ASP32 Figure 4. sample conformational snapshots of the 'Pinky finger' region (residue 1-50) taken from the simulations of Xyn11A with a) no protonation, b) ASP11, c) ASP32 and d) ASP11/ASP32 protonation. Residues 11 and 32 were represented by balls and sticks the other hairpin-loop residues (B1/A1: residues 11-15 and A2/B2: residues 31-35) were represented by licorice, in which C, H, O and N atoms were shown in grey, white, red and blue, respectively. The protonated sites were highlighted in purple and the hydrogen bonds were represented in dashed lines under the conditions that donor-acceptor distances < 3.5 Å and donor-hydrogen-acceptor angle < 45°.

Conclusions
A series of atomistic molecular dynamics simulations were performed to predict the positive ion distribution on the 'Pinky-finger' region of a GH11 xylanase enzyme, and to investigate the conformational changes due to the protonation at the predicted sites. Conformational analysis showed that ASP32 protonation caused ARG31 to form a dynamic hydrogen bonding networks with the B1/A1 hairpin loop (residues 11-15). This non-native local stabilization caused the 'Pinky-finger' region to be misfolded and became the nucleation point for total destabilization. Simulation results provided information for further protein engineering to increase the enzyme acido-tolerability, as the removal of the misfolding source -residue ARG31 -was suggested.