Zora
Střelcová
ab,
Petr
Kulhánek
ab,
Martin
Friák
acd,
Helge-Otto
Fabritius
c,
Michal
Petrov
c,
Jörg
Neugebauer
c and
Jaroslav
Koča
*ab
aCEITEC – Central European Institute of Technology, Masaryk University, Kamenice 5, 625 00 Brno, Czech Republic. E-mail: jaroslav.koca@ceitec.muni.cz
bNational Centre for Biomolecular Research, Faculty of Science, Masaryk University, Kamenice 5, 625 00 Brno, Czech Republic
cMax-Planck-Institut für Eisenforschung GmbH, Max-Planck-Straße 1, 40237 Düsseldorf, Germany
dInstitute of Physics of Materials, Academy of Sciences of the Czech Republic, Žižkova 22, 616 00 Brno, Czech Republic
First published on 18th March 2016
Chitin is one of the most abundant structural biomolecules in nature, where it occurs in the form of nanofibrils that are the smallest building blocks for many biological structural materials, such as the exoskeleton of Arthropoda. Despite this fact, little is known about the structural properties of these nanofibrils. Here, we present a theoretical study of a single chitin molecule and 10 α-chitin nanofibrils with different numbers of chains in an aqueous environment that mimics the conditions in natural systems during self-assembly. Our extensive analysis of the molecular dynamics trajectories, including free energy calculations for every model system, reveals not only the structural properties of the nanofibrils, but also provides insight into the principles of nanofibril formation. We identified the fundamental phenomena occurring in the chitin nanofibrils such as their hydrogen bonding pattern and resulting helical shape. With increasing size, the nanofibrils become increasingly stable and their structural properties approach those of crystalline α-chitin if they consist of more than 20 chains. Interestingly, this coincides with the typical size of chitin nanofibrils observed in natural systems, suggesting that their evolutionary success was at least partially driven by these specific structure–property relations.
Typical α-chitin crystallites occurring in the cuticle have a polygonal section profile and consist of 18–25 chitin molecules with a diameter between 2 and 5 nm and lengths of about 300 nm.6 They are always surrounded by a sheath of proteins, resulting in rod-shaped nanofibrils that have been identified as the basic structural unit of the cuticle in a large variety of species from all groups of Arthropoda.6–9 Potential chitin-organizing proteins currently under discussion are a class of evolutionarily conserved cuticle proteins featuring the so-called Rebers and Riddiford (R&R) domain9 such as resilin, and the obstructor proteins10 that possess tachycitin-like chitin-binding domains. The long-known sterical compatibility between the amino groups of the acetylglucosamine and amino acid residues of the proteins have been proposed to facilitate the formation of hydrogen bonds between the molecules.11 However, the exact functions of these proteins have not been identified yet. Chitin–protein nanofibrils further organize into stacks of planar arrays that form the so-called twisted plywood structure commonly found in all exoskeletal elements of Arthropoda. While the ground-state structure of α-chitin,7,8,12,13 its mechanical properties and the impact on the overall mechanical performance of the composite cuticle have been extensively studied both theoretically and experimentally,14,15 little is known about how it organizes to form the chitin–protein nanofibrils. During the process of cuticle formation,16 chitin molecules are produced by the membrane inserted enzyme chitin synthase.17 The enzyme catalyses the β-1-4 linkage between N-acetylglucosamines supplied from the cytoplasm. The formed chitin molecules are extruded across the plasma membrane into the extracellular space.18–21 Here, the chitin and chitin-binding proteins present in the extracellular space organize into the typical nanofibrils that form the exoskeleton, most likely in a self-assembly process. This process involves the assembly of chitin molecules into larger fibrils and their interaction with chitin-binding proteins, which are the fundamental steps not only for the formation of the exoskeleton itself, but also establish the fundamentals that lead to the elaborate function-specific structural differentiation in the final tissue observed in many species of arthropods.
Even though crystalline α-chitin has been intensively studied in past years,12,14,22,23 its structural properties remain unclear. Likewise experimental evidence concerning nanofibril composition, structural parameters and hierarchical structure were never further expanded to explain the self-organization of chitin nanofibrils. Since suitable in vitro and even more so in vivo experimental approaches are barely practicable in terms of elucidating the basic principles of the self-assembly of single chitin molecules to α-chitin nanofibrils, this study employs theoretical approaches using molecular dynamics in fully solvated systems to shed light on these processes. There are only few studies described in the past that employ molecular dynamics for description of chitin properties. For example, properties of a single chitin chain was studied by molecular dynamics approach in aqueous environment24 and during adhesion to montmorillonite clay.1 More complex studies try to describe adsorption of chitosan to chitin nanocrystallites,25 thermodynamics of α-chitin de-crystallization,26 and properties of chitin–protein composites.27 Recently, a coarse grain model of crystalline α-chitin was derived from experimental data and atomistic molecular dynamics simulations.28
However, to our best knowledge, none of these studies did characterize self-assembly of α-chitin nanofibrils and structures of possible intermediates in a comprehensive way. Due to methodological constraints, the self-assembly process of α-chitin nanofibrils in vivo as well as in vitro cannot be fully covered in an integral modeling experiment due to the size of simulated models, which are simply beyond current technical possibilities. Therefore, the process is studied in reverse by observing what the properties of already assembled chitin nanofibrils of various dimensions are.
Here, we evaluate the geometry of single chitin molecules and chitin nanofibrils of various sizes (Fig. 1), correlate the observed global structural features with local geometrical parameters, and analyze the driving forces responsible for the self-assembly of nanofibrils in terms of free energy changes. The studied nanofibril models were chosen such as to reflect successive intermediate fibril configurations that presumably occur during the self-assembly process leading to the nanofibril size reported in previous studies of arthropod cuticle.29,30 This corresponds to the 6760 system (Fig. 1), which should be approximately 3 nm in diameter. Starting from the 6760 system, we considered one bigger and nine smaller systems. The shapes of smaller systems were selected in such a way that they reasonably span all important chain rearrangements in the both a and b crystallographic directions while avoiding simulations of too similar systems.
Molecular dynamics simulations performed with the Amber 12 package31 were employed to evaluate the structural and dynamic behavior of the prepared models. The chitin nanofibrils were described with a GLYCAM 06g force field.32 To mimic the native conditions in the chitin nanofibril assembly zone, each structure was put in a truncated octahedral box, which was filled with water molecules described by the TIP3P water model.33 Long-range interactions under the periodic boundary conditions were treated by the particle-mesh-Ewald method,34 with a direct summation cut-off of 9.0 Å. The same cut-off distance was used for the treatment of van der Waals interactions. Each system was equilibrated prior to the production of simulations. The system was slowly heated from 5 K to 300 K. During the equilibration, the solute was kept fixed by positional restraints with decreasing strength from 50.0 to 0.0 kcal mol−1 A−2, with the initial structure used as a reference structure. The barostat set to 100 kPa was enabled at later stages of the equilibration to reach the proper density of the water environment.
Production simulations were performed at a constant temperature of 300 K and a pressure of 100 kPa. The lengths of all bonds containing hydrogen atoms were kept fixed, enabling an integration time step of 2 fs to be used. Temperature was kept fixed by Langevin thermostat employing a collision frequency of 1.0 ps−1. Pressure was maintained by a weak coupling barostat with a feedback time constant of 1.2 ps. The production dynamics were run on GPU accelerators.35 The length of each production run was 200 ns, configuration snapshots were taken each 10 ps.
Trajectory analyses were performed in the ptraj and cpptraj modules36 of the Amber package or using in-house developed scripts. The root-mean square deviation (RMSD)37 was calculated for each snapshot along the trajectory to the reference structure, which was the initial structure (e.g. crystal-like structure).
Hydrogen bonds were analyzed by evaluating number of contacts between the H and A atoms in patterns of the following form D–H⋯A–Z that might be considered as a stabilizing interaction, where the three dots denote the hydrogen bond.38 The number of contacts was calculated using a Fermi switch function:
The axial chirality was described as the pseudo-dihedral angle α between four points. The second and third points are the centers of mass of the sugar rings aligned in the ab terminal planes. These two points represent the helical axis of nanofibrils. The first and fourth points are the centers of mass of the top-edge sugar rings (in the b direction). The geometry of glycosidic bonds was evaluated using the ψ(C5n, C4n, O4n, C1n+1) and ϕ(C4n, O4n, C1n+1, O5n+1) dihedral angles. Rotamers of terminal CH–CH2OH groups were monitored using the dihedral angle δ(C4, C5, C6, O6) and orientations of terminal CH3CONH–CH groups were monitored using the dihedral angle ε(C1, C2, N2, C2N). Ring conformations were analyzed by measuring the ϕ and θ ring-puckering coordinates for sugar rings (C1, C2, C3, C4, C5, and O5) and pseudo-rings (O4n, C4n, C3n, O3n, O5n+1, C1n+1).39
The binding strength between chitin chains was evaluated by means of the MM/PBSA40–42 as implemented in Amber 12. The method is based on post-processing the trajectories obtained during production molecular dynamics to assess the absolute free energy of the solute (chitin nanofibril). The absolute free energy G is composed of two major contributions: (a) the free energy of the solute in vacuum Gvac for a solute configuration ensemble generated in the solvent and (b) the solvation free energy of the solute Gsol, representing the vertical transfer from vacuum to solvent (thus this energy does not include any deformation energy). The first term is calculated from the solute enthalpy H, which is taken as a time average of the total solute energies EMM calculated by means of molecular mechanics (MM), and entropy S. The calculation of the entropy term will be discussed later in more detail.
G = Gvac + Gsol = H − TS + Gsol = 〈EMM〉 − TS + 〈Gelsol〉 + 〈Gnesol〉 |
The solvation free energy is composed of two contributions: electrostatic Gelsol and non-electrostatic Gnesol. The electrostatic part was evaluated by employing a Poisson–Boltzmann (PB) method41 to describe the solvent implicitly. The PB method requires the definition of the atomic radii separating the interior and exterior of the system. In this work modified Bondi parameters40 were used, together with an interior relative dielectric constant set to 1.0. The exterior relative dielectric constant was set to 80.0, which corresponds to water. The non-electrostatic contribution was calculated by employing a surface tension model using the solvent accessible surface area43,44 with a default tension constant for the given atomic radii sets (0.0072 kcal mol−1 A−2). Due to the computational complexity of the MM/PBSA, snapshots taken every 1 ns were included in the analysis of the last 150 ns.
The change in the total reaction free energy related to a single chitin chain (the structure labeled 1000) was calculated via the following equation:
The essential dynamics analysis (EDA) method45,46 was used to identify correlated atomic motions occurring during simulations. It is a post-processing method for molecular dynamics trajectories that can also be used to filter out thermal noise from analyzed trajectories. In the first step, the covariance matrix was calculated from fluctuations of the C1, C3, and C5 atoms around their average positions, employing trajectory snapshots taken every 10 ps. Then, eigenvalues and corresponding eigenvectors of the covariance matrix were evaluated. The eigenvalues represent the total mean square fluctuation of the system along the corresponding eigenvectors. Eigenvectors were employed to divide overall motion along the trajectory into separate modes. The modes with greater eigenvalues represent slow motion – essential motion, whereas modes with smaller eigenvalues represent rapid motion – thermal motion. New trajectories containing only correlated atomic motions individually projected from the source trajectory on five eigenvectors with the largest eigenvalues were calculated together with projection factors p.
For further analysis, the trajectories were coarse-grained. These simplified trajectories were employed in the vibrational entropy calculations. Each N-acetylglucosamine unit was represented by a point (CG bead) with a mass equal to the total mass of the N-acetylglucosamine unit, and its position was taken as the center of mass of all atoms in the given unit. The intention behind this procedure was to decrease the number of degrees of freedom while keeping the overall motions of systems close to representations with atomic resolution.
The entropy contributions were calculated separately for translation, rotation and vibrational motions using the cpptraj module from the Amber package. The translational St and rotational Sr contributions were calculated using statistical mechanical formulas for an ideal gas model at a temperature of 300 K and a pressure of 100 kPa. The vibrational entropy Sv was calculated for the same temperature using a quasi-harmonic approach.47 To eliminate the impact of trajectory length on results, the vibrational entropy was calculated from varying blocks of the trajectory with lengths of 50 ns, 75 ns, and 150 ns and then extrapolated to infinite time.48
Figures showing 3D structures of the studied systems were prepared in PyMol49 and VMD.50
![]() | ||
Fig. 3 Average structures of simulated nanofibril systems. Building blocks are color-coded to highlight structural features. |
The overall flexibility of the different nanofibril models was monitored by root-mean square deviation (RMSD, Fig. 4 and S3†). This property quantifies the structural deviation of the simulated models from ideal, crystal-like structures. The larger this value is at a given time, the more the model structure deviates from the initial structure. The single chitin chain exhibits a broad fluctuation ranging from 3.3 to 15.1 Å, with a maximum at 5.3 Å. This indicates that the structure is rather flexible. Despite this fact, the structure maintains a straight shape without tendencies to fold into more compact structures. This was expected, as the single chitin chain would be the most flexible due to the lack of surrounding partners that would stabilize the structure. Surprisingly, it was found that some of the larger nanofibrils, namely the 1110, 2200, 6000, 6600 systems, deviate even more from the crystalline configuration, indicating an even higher flexibility. To explain this observation, the hydrogen bond network in nanofibrils was analyzed, as this is thought to be one of the key interactions that influence their shape, stability and mechanical properties.
The hydrogen bond analysis was performed for specific and non-specific pairs of hydrogen bond donor and acceptor atoms. The specific hydrogen bonds represent geometrically necessary interactions occurring in the crystalline α-chitin while the non-specific hydrogen bonds represent possible interactions that might coincidentally occur during its evolution or its deformation. The specific hydrogen bonds are represented by the O2N–H2N, O5–H3O, O2N–H6O, and O6–H60 atom pairs, and their statistics are summarized in Table 1.
System | n | l | O5–H3O | O2N–H2N | ||||
---|---|---|---|---|---|---|---|---|
〈HB〉 | σ(HB) | 〈HB〉/n | 〈HB〉 | σ(HB) | 〈HB〉/l | |||
a n – number of chitin chains; l – number of stacking layers; 〈HB〉 – average number of hydrogen bonds; σ(HB) – variance (fluctuation) of number of hydrogen bonds; 〈HB〉/n – average number of hydrogen bonds related to a single chitin chain; 〈HB〉/l – average number of hydrogen bonds related to a single stacked layer of chitins. | ||||||||
1000 | 1 | 0 | 14.04 | 1.74 | 14.04 | 0.00 | 0.03 | — |
1110 | 3 | 0 | 44.01 | 2.76 | 14.67 | 2.04 | 1.12 | — |
2200 | 4 | 2 | 71.76 | 1.76 | 17.94 | 36.11 | 1.62 | 18.06 |
2220 | 6 | 3 | 110.00 | 1.64 | 18.33 | 56.18 | 1.61 | 18.73 |
3300 | 6 | 4 | 110.70 | 1.50 | 18.45 | 72.87 | 1.92 | 18.22 |
3330 | 9 | 6 | 166.40 | 1.72 | 18.49 | 113.40 | 2.19 | 18.90 |
4400 | 8 | 6 | 146.80 | 1.90 | 18.35 | 108.50 | 2.70 | 18.08 |
6000 | 6 | 5 | 109.00 | 1.87 | 18.17 | 89.91 | 2.64 | 17.98 |
6600 | 12 | 10 | 221.30 | 2.13 | 18.44 | 182.10 | 3.66 | 18.21 |
6760 | 19 | 16 | 353.40 | 2.03 | 18.60 | 308.30 | 3.07 | 19.27 |
8998 | 34 | 30 | 634.00 | 2.52 | 18.65 | 567.30 | 5.22 | 18.91 |
System | n | l | O2N–H6O | O6–H60 | ||||
---|---|---|---|---|---|---|---|---|
〈HB〉 | σ(HB) | 〈HB〉/n | 〈HB〉 | σ(HB) | 〈HB〉/l | |||
1000 | 1 | 0 | 0.01 | 0.08 | — | 0.05 | 0.18 | — |
1110 | 3 | 0 | 9.01 | 4.12 | — | 3.24 | 1.52 | — |
2200 | 4 | 2 | 9.13 | 2.65 | 4.56 | 10.89 | 4.29 | 5.45 |
2220 | 6 | 3 | 13.48 | 2.92 | 4.49 | 13.48 | 2.92 | 4.49 |
3300 | 6 | 4 | 18.66 | 3.32 | 4.67 | 31.79 | 2.95 | 7.95 |
3330 | 9 | 6 | 28.63 | 4.00 | 4.77 | 59.01 | 3.91 | 9.84 |
4400 | 8 | 6 | 26.09 | 3.99 | 4.35 | 44.27 | 4.49 | 7.38 |
6000 | 6 | 5 | 16.19 | 3.78 | 3.24 | 12.53 | 2.77 | 2.51 |
6600 | 12 | 10 | 45.63 | 5.85 | 4.56 | 67.58 | 5.90 | 6.76 |
6760 | 19 | 16 | 89.73 | 5.92 | 5.61 | 143.50 | 7.07 | 8.97 |
8998 | 34 | 30 | 192.70 | 9.57 | 6.42 | 266.30 | 10.19 | 8.88 |
The largest number of hydrogen bonds per chain was detected for the O5–H3O atom pair. This type of hydrogen bond is responsible for stabilization along the c-axis of a single chitin chain due to an interaction between the O5 atom and H3O–O3 hydroxyl group from adjacent sugar rings. Except for the 1000 and 1110 fibrils, all systems exhibit values that range from 17.9 to 18.6. These values are close to the theoretical limit, which is 19.0 for twenty-unit-long chains. The 1000 and 1110 systems exhibit only slightly smaller values (about 14). These might be caused by the absence of stabilizing adjacent chitin chains in the ac plane, and the resulting complete exposure of both systems to the water environment. These hydrogen bonds are also responsible for the formation of seven-membered pseudo-rings, which are roughly in the chair conformation. Together with the sugar rings, which are also in the chair conformation, they form a thermodynamically stable belt-like structure of mutually interconnected rings. This arrangement makes single-chitin chains resistant to axial rotations along glycosidic bonds, which prevents a short single chain collapsing into a more compact structure. This is in agreement with the behavior observed in the simulation of the 1000 system. The conformation of rings was further quantified by the ring-puckering ϕ and θ angles (Fig. S4 and S5†). The polar θ angle is close to 180.0° for both kinds of rings, which correspond to the observed chair conformation. The azimuthal ϕ angle shows different distributions for every simulated system, but as the polar θ angle is always close to 180.0°, these changes actually reflect only very small geometrical changes. The only significant deviation of the polar θ angle from its ideal value of 180.0° was observed for the 1000, 1110, 2200, and 6000 systems. In these systems, the ring enclosures formed by the O5 and H3O atoms are not buried in the structure, but located on the surface and thus fully exposed to water, making the enclosures significantly weaker.
The second most common type of hydrogen bond is formed between stacked chitin chains (in the a direction) by interaction between the carbonyl O2N and amide H2N atoms. The number of hydrogen bonds per number of stacked layers is between 17.98 and 19.27, which is close to the theoretical limit of 20.0. The value slightly increases with an increase in the number of adjacent stacked layers in the b direction. For example, the 2200 system exhibits a value of 18.06, while the 2220 system exhibits a value of 18.73. However, this effect is less pronounced if the number of layers stacked in the a direction increases (see the 3300, 3330 and 6000, 6600 systems). From previous studies52 it follows that this type of hydrogen bond is responsible for maintaining the structural integrity of chitin crystallites in the a direction, which is in agreement with the behavior observed in all our simulated systems that contain stacked chains.
The remaining specific hydrogen bonds are formed by the O2N–H6O and O6–H60 atom pairs. Their values, related to either the number of chains or the number of stacked layers, are not as constant as in the two previous cases (Table 1), but strongly depend on the nanofibril size and shape. One possible explanation for this might be that these hydrogen bonds are weaker than the previous ones and could thus be influenced by the water environment at the surface of the nanofibrils. To confirm this hypothesis, we tried to rationalize this observation by including the effect of the solvent into the model. Instead of entering the number of hydrogen bonds purely due to the total number of stacked layers, the stacked layers were counted separately for the situation where they are exposed to the solvent (they are close to the surface of the nanofibril – lw) or when they are located inside the nanofibril (li). To make the model robust, each edge of stacked layers was counted independently. The relationship was then considered as a simple linear model of the following form:
〈HB〉 = nwlw + nili |
We also focused on the analysis of other possible hydrogen bonds that do not commonly occur in crystalline chitin. The aim was to find possible correlations between the presence of such non-specific hydrogen bonds and the observed structural shape of the simulated nanofibrils. However, the analysis of the O6–H3O, O6–H2N, O3–H6O, O3–H2N, O2N–H3O atom pairs did not confirm such correlations, as the number of these hydrogen bonds related to a single chain is negligible (smaller than one in all cases except for the 1110 system, see Table S4†).
During the simulation of the 1110 model system, we observed a structural reorganization of the nanofibril. This is characterized by a very large change in the RMSD, the largest observed deviation among all the tested systems. The RMSD shows a broader distribution in the range from 5.7 to 15.7 Å with a maximum at 8.9 Å (Fig. 4) for the 200 ns-long simulation. Detailed analysis of the trajectory showed that the structure of the model underwent a rearrangement, in which the non-covalent (mostly dipole–dipole) interaction between adjacent chains in the b direction is disrupted. This occurred during the equilibration stage of the simulation, which indicates that this kind of interaction is very weak and was the reason for their replacement with hydrogen bond interactions. To quantify this process, the evolution of emerging hydrogen bonds (HB) over time was calculated (Fig. 5). During the 200 ns-long simulation, we observed a sudden change in the structure of the nanofibril chain, which motivated us to extend the simulation to 500 ns to observe possible further reorganization. The hydrogen bonds were monitored for the O2N–H2N, O2N–H6O, and O6–H6O atom pairs as representatives of those specific hydrogen bonds. The evolution of the fourth specific hydrogen bond for the O5–H3O atom pair is not shown, since the count was constant throughout the whole simulation (around 44.1, see Table 1). The evolution of non-specific hydrogen bonds formed by the O6–H2N and O2N–H3O atom pairs was also monitored.
![]() | ||
Fig. 5 Number of hydrogen bonds emerging during simulation of the 1110 model. Lines represent running averages over 5 ns time intervals and vertical bars are variances in the same time interval. |
Two changes occurred during the simulation: the first one happened in the interval from the 100th to 150th ns and the second change, which was more rapid, at about the 425th ns. These two changes split the simulations into three regions, for which representative structural snapshots are shown in Fig. S7.† During the first change, the structure that still resembled the starting structure collapsed and became more compact, which resulted in an increase in the number of hydrogen bonds, especially non-specific ones. From the specific hydrogen bonds, only the O2N–H6O type prevails. However, a very slow increase in two other types is also visible in the region from the 150th to 400th ns. The average increase is 0.29 and 0.28 hydrogen bonds per 100 ns for the O2N–H2N and O6–H6O atom pairs. The second change is accompanied by an increase in the number of O2N–H2N hydrogen bonds. This resulted in the formation of regions in the structure where the chains change their position to the stacked conformation (Fig. S7,† highlighted with ovals). This stacking is not perfect yet, as the chains are not properly aligned vertically in some parts (they are slightly displaced relative to each other in the c direction). In addition, the chains even cross each other and stack in an inappropriate order in one part of the structure. Taken together, the observed changes suggest that the 1110 system slowly transforms into the 3000 system, where the stacking is the most dominant interaction due to specific hydrogen bonds. From the analysis of hydrogen bonds in the other simulated models, we can predict the numbers of hydrogen bonds for the 3000 arrangement to be: 36.44 (O2N–H2N), 4.16 (O2N–H6O), and 6.56 (O6–H6O). The number of non-specific hydrogen bonds should be close to zero. At the end of the 1110 simulation, the final values were: 12.29 (O2N–H2N), 7.25 (O2N–H6O), 3.83 (O6–H6O), and 7.33 (O6–H2N). From the comparison of the predicted and current value of the leading structural descriptor (O2N–H2N), we can assume that the transformation is approximately in the first third. The further rearrangement must involve the structure re-folding due to misaligned stacking. This was the reason why we did not continue with the simulation, because a significantly longer simulation time would be required.
The other simulated systems did not undergo any significant rearrangements like the one observed for the 1110 system. However, some of them still deviated significantly from the crystal-like geometries. This is especially noticeable for the 2200, 6000, and 6600 systems. The cause of this behavior is a twist distortion along the c-axis of the structure (see Fig. 3). This distortion makes all the systems axially chiral (in addition to the configurational chirality introduced by the stereogenic centers in the building saccharide units). As a measure of axial chirality, the pseudo-dihedral angle α was calculated and its distribution is shown in Fig. 6. Its average values, including the corresponding fluctuations, are summarized in Table 2.
System | 〈α〉 | σ(α) | 〈ψ〉 | σ(ψ) | 〈ϕ〉 | σ(ϕ) |
---|---|---|---|---|---|---|
1000 | — | — | −129.0 | 20.5 | −78.0 | 11.1 |
1110 | — | — | −130.6 | 22.4 | −79.1 | 10.8 |
2200 | 236.2 | 48.9 | −139.7 | 10.5 | −83.4 | 8.2 |
2220 | 100.3 | 22.5 | −145.0 | 9.7 | −86.8 | 7.8 |
3300 | 117.7 | 11.0 | −145.5 | 9.1 | −86.2 | 7.5 |
3330 | 59.2 | 19.9 | −148.0 | 8.4 | −88.0 | 7.2 |
4400 | 112.4 | 10.8 | −146.4 | 8.9 | −86.2 | 7.1 |
6000 | 169.4 | 11.4 | −142.7 | 9.6 | −85.3 | 7.5 |
6600 | 99.9 | 9.0 | −146.9 | 8.5 | −86.6 | 7.1 |
6760 | 15.4 | 5.8 | −149.9 | 7.3 | −89.1 | 6.7 |
8998 | 11.3 | 6.3 | −150.9 | 7.8 | −89.2 | 7.1 |
The factor that influences the degree of axial chirality most significantly is the number of stacked layers of chitin chains interacting in the b direction. The axial chirality decreases when a nanofibril is growing in this direction (compare the 2200 → 2220, 3300 → 3330, and 6000 → 6600 systems). The highest axial chirality among the simulated systems is exhibited by the 2200 and 6000 systems. The helical twist of 180° along the entire length is particularly obvious in the average structure of the 6000 system shown in Fig. 3. In an attempt to explain this behavior at the molecular level, we focused on possible correlations of the pseudo-dihedral angle α with different structural descriptors (Fig. S8–S10 and Tables S5 and S6†).
The significant correlation between the pseudo-dihedral angle α and structural descriptors was found for dihedral angles describing the geometry of the glycosidic bonds. The dihedral angles ψ and ϕ exhibit single Gaussian distributions (Fig. S9†). The peak maxima (average values) (Table 2) then correlated well with the pseudo-dihedral angle α with R2 of 0.97 and 0.99, respectively. The obtained linear models (Fig. S10†) enable the extrapolation of ψ and ϕ to α = 0, which should correspond to the situation in crystalline α-chitin (no axial chirality is observed in the crystal). The comparison of the obtained values with the data extracted from two available crystal structures and ab initio calculations is summarized in Table 3. Although the values ψ and ϕ were not determined for the real structure of crystalline α-chitin, but extrapolated from the behavior of smaller chitin nanofibrils, their correspondence with the experimental as well as ab initio data is very good. They are approximately in the middle between two values determined by X-ray crystallography and are also close to the values determined by ab initio calculations. This can be taken as a performance validation of the force field employed for the description of chitin nanofibrils.
We also checked if the observed agreement between our simulations, experimental, and ab initio data cannot be consequence of specific force-field parameters that would force the system to adopt configuration with values of the ψ and ϕ dihedral angles equal to values found by the extrapolation. We found that minima on the potential energy surface describing rotation on both glycosidic bonds do not coincide with the extrapolated values of dihedral angles (Fig. S11†). This indicates that observed behavior is a result of orchestration of all type of interactions among them dihedral, electrostatic and van der Waals interactions playing predominant role.
The twist was also observed in simulations of cellulose nanofibrils.53–56 However it was shown that on longer timescales (∼800 ns) the twist can disappear due to reorganization of hydrogen bonds that starts at the solvent/solute boundary and slowly propagates into the body of cellulose nanofibril. Moreover it was shown that the phenomena is force-field dependent.54 On timescale of our simulations (∼200 ns) we did not observe such disappearance of the twist. Splitting of distribution of the pseudo-dihedral angle α in the 3330 as well as 8998 systems (Fig. 6) is rather caused by an increase of the twist for a short period during the simulations. To further support this different observation, we extended length of simulation for 6000, 4400, and 6760 systems up to 300 ns (Fig. S12†). Once again, no sudden change or decrease in the twist was observed.
A recent study from the same authors performed on cellulose nanofibrils of various sizes and on shorter timescales (∼200 ns)55 shows results that are very similar to our work. They found that the twist is decreasing with increasing size of nanofibril. Nevertheless they were not able to recognize that magnitude of the twist might exhibit anisotropy due to size of nanofibrils along the a and b directions because their nanofibrils were symmetrical. Also the molecular origin of the twist will be probably different as the authors identified HO2 group as a key player, which is however absent in chitin as this position is occupied by the N-acetylamide group.
Another correlation was found between the pseudo-dihedral angle α and average values of the ring-puckering coordinate θ for the pseudo-rings formed by the hydrogen bonds between the O5–H3O atom pairs. As a consequence, the pseudo-dihedral angle also correlates with the number of O5–H3O hydrogen bonds per chain (Table S6†). R2 is 0.95 and 0.94, respectively.
The correlations found indicate that the twist distortion of the simulated nanofibrils is a consequence of the internal axial chirality of individual chitin chains. The internal axial chirality is caused by deviations of the ψ and ϕ dihedral angles from their ideal values observed in the crystal structure. In our model, these extrapolated ideal average values are −151.0° and −89.4°, while in the single chitin chain (the 1000 system) they are only −129.0° and −78.0°, respectively. The internal axial chirality of single chains is suppressed in larger chitin nanofibrils by non-covalent interactions between adjacent stacked layers of chitin chains. These probably include dipole–dipole interactions between acetamide groups in the b direction. As a consequence of distortion along glycosidic bonds, the number of O5–H3O hydrogen bonds that stabilize a belt-like structure of chitin chains was also found to correlate well with the overall axial chirality, together with the θ puckering coordinate describing such pseudo-ring conformations.
The axial chirality of chitin nanofibrils is also reflected in key soft-mode movements detected by the essential dynamics analysis of coarse-grained trajectories employing the C1, C3, and C5 atoms. This analysis separates the overall dynamic motion into linearly independent motions scored by their softness (e.g. how the energy of the system changes with the particular movement). Projections of the trajectories to the five essential modes with highest (softest) eigenvalues are summarized in Fig. S13.† Non-uniform (Gaussian) distributions of projections were found for several systems (most notably for 1110, 2200, 3330, 6600, and 8998). This indicates that the motion is too slow to be properly sampled on the timescale of the analyzed simulations. With the 1110 system, the non-uniform distributions of several modes are caused by the complex structural rearrangement, which was already discussed. With the 2200, 3330, 6600, and 8998 systems, the non-uniform distribution is caused by slow axial distortion, which is also visible in the distribution of the pseudo-dihedral angle α (Fig. 6, most notably the split distribution for the 3330 system).
The slowest essential modes show changes in the axial chirality for all systems (apart from 1000 and 1110), indicating that a twist distortion is energetically the most favorable deformation (Fig. 7, M1). The other modes show various types of distortion along the c direction, with the 6760 system they are variants of shear distortions (Fig. 7, M2). Since the 6760 system most closely resembles the geometry of chitin nanofibrils observed in arthropod exoskeletons, these properties might be important for the interaction between chitin nanofibrils and chitin-binding proteins in natural systems. While chitin is very stiff in the c direction due to the covalent nature of glycosidic bonds supported by the hydrogen bond network, it is rather soft in the other directions, which might facilitate small structural adjustments that enable proper alignment to the protein binding sites in living systems.
![]() | ||
Fig. 7 Movements along two essential modes for the 6760 system. The M1 mode shows a twist distortion while the M2 mode shows a shear distortion. |
The fact that the chitin fibrils in mature and developing exoskeleton are always confined by various chitin-binding proteins makes it difficult to obtain direct experimental evidence for axial chirality in the α-chitin nanofibrils, since the interactions between chitin and proteins are likely to influence the structural conformation of the chitin fibrils. However, synchrotron X-ray micro-diffraction experiments have shown that in arthropod exoskeletons the crystallographic axes a and b of the chitin fibrils rotate around the longitudinal axis c.58 This supports our results, since this rotation might very well originate from an initial axial distortion of the newly assembled chitin nanofibrils which is retained despite the addition of the protein coating.
The self-assembly of chitin chains into chitin nanofibrils is controlled by the kinetics and thermodynamics of individual assembly steps. The kinetics of self-assembly is very hard to describe by means of computational techniques, as the process is composed of several reaction pathways which are difficult to capture and characterize. An example of such a process partially occurred in the simulation of the 1110 model. In contrast, thermodynamics can be easily described by means of various computational techniques. In this work, we selected the MM/PBSA method. Using this method, we calculated the reaction free energies released during nanofibril formation.
The changes in reaction enthalpy in the gas phase ΔHr/n show that the process is strongly exothermic. As we use molecular mechanics (pair-wise additive approach) to describe the system, it is possible to further decompose the enthalpy into individual energy terms (Table S7,†Fig. 8 and S14†). From the decomposition, it follows that the major driving forces in the binding are equally contributing electrostatic and van der Waals forces. This is very surprising, as we expected that the electrostatic contribution would be the major player due to the emerging hydrogen bond network and dipole–dipole interactions of acetamide groups. The value of ΔHr/n decreases with an increasing number of chains to the limiting value (extrapolated for an infinite number of chitin chains), which we estimated to be about −651.1 ± 58.4 kcal mol−1. This behavior was expected, because with an increasing number of chitin chains the nanofibril structures approach the structure of crystalline chitin (the ideal structure) where the overall interaction network is maximized.
While the interaction between chitin chains is strongly favorable towards the assembly of chitin nanofibrils (negative enthalpy) we found that the solvent destabilizes their mutual binding. This is caused by the relatively strong interaction between nanofibrils and solvent molecules, and the fact that the contact of individual disassembled chitin chains with the solvent is higher than in an assembled nanofibril. The contribution of the solvent to the binding is expressed by the value of 〈ΔGsol,r〉/n. The value of 〈ΔGsol,r〉/n, which is the sum of 〈ΔGel,r〉/n and 〈ΔGne,r〉/n, increases up to the limiting value (Table S8,†Fig. 8), which we estimated to be about 394.7 ± 34.7 kcal mol−1. This value is smaller than the limiting value for the reaction enthalpy in the gas phase, confirming that the destabilization by the solvent is weaker than the mutual interaction between chitin chains.
The change in total reaction free energy per chitin chain is provided in Fig. 8. Its extrapolated limiting value is about −256.6 ± 24.2 kcal mol−1. The initial very steep change in the total binding energy reaches an energy plateau with a very slow energy decay towards the limiting energy starting from approximately 20 chains (n = 20). From this perspective, the 676 and 8998 nanofibrils reached dimensions at which the binding enthalpy and desolvation energies are already maximized. From the enthalpy point of view, this means that the shape of the nanofibril and formed internal interaction network reached a situation very close to the crystalline chitin. From a desolvation point of view, as the number of chitin chains increases, the total surface of chitin chains that lose contact with the solvent (the proportion of the chitin chains that are desolvated upon binding) becomes predominant compared to the nanofibril surface (the proportion of the chitin chains that remain solvated). The value of 〈ΔGsol,r〉/n must then converge to the opposite value of the solvation energy of a single chitin chain. Indeed, the value of 〈ΔGsol,r〉/n converges to 394.7 ± 34.7 kcal mol−1 while the solvation energy of the 1000 model is −399.1 kcal mol−1.
The beginning of the energy plateau coincides with the size of the nanofibrils occurring in the natural chitin nanocomposites, which usually consist of 18–25 chitin molecules.6 It can be assumed that the evolution of chitin nanofibrils with this geometry was at least partially driven by the fact that they are the smallest nanofibrils with crystalline chitin-like properties, which combine small size for structural versatility with excellent chemical and mechanical properties.
We also use our data to predict changes in reaction entropy. Due to the length of the trajectories, we only obtained qualitative data on the simplified coarse-grained models of the simulated systems (Table S9 and Fig. S15†). The obtained data show a decrease in translational, rotational and vibrational contributions to the total reaction entropy per single chitin chain (ΔSr/n) with an increasing number of chains, which could be expected as the structures become more compact and rigid. Due to coarse-graining of the analyzed trajectories, the model does not correctly account for all entropy changes, such as for example the decrease in hydroxyl group rotational entropy upon the formation of hydrogen bonds. Thus in future studies it will be necessary to employ more advanced techniques for entropy calculations, such as the two-phase thermodynamic (2PT) model.59 Despite bottlenecks in the methodology used, entropy does not change the progress of the reaction free energy discussed above. The entropy converges towards about −275 cal mol−1 K−1, (a very crude estimate due to the irregular progress of entropy), which corresponds to a free energy contribution of about +82.5 kcal mol−1. This, together with the reaction free energy can be used to predict the total contribution to binding per single unit of chitin chain (e.g. N-acetylglucosamine), which is about −8.7 kcal mol−1. Since this value is in very good agreement with typical binding energies including saccharide moieties,60 our methodology is suitable for providing consistent predictions of binding energies and affinities at least at a qualitative level.
Footnote |
† Electronic supplementary information (ESI) available: Tables S1–S9, Fig. S1–S15. See DOI: 10.1039/c6ra00107f |
This journal is © The Royal Society of Chemistry 2016 |