The structure and dynamics of chitin nanofibrils in an aqueous environment revealed by molecular dynamics simulations

CEITEC – Central European Institute of Tec 625 00 Brno, Czech Republic. E-mail: jarosl National Centre for Biomolecular Research Kamenice 5, 625 00 Brno, Czech Republic Max-Planck-Institut für Eisenforschung Düsseldorf, Germany Institute of Physics of Materials, Academy o 22, 616 00 Brno, Czech Republic † Electronic supplementary information (E See DOI: 10.1039/c6ra00107f Cite this: RSC Adv., 2016, 6, 30710


Introduction
In contrast to most man-made structural materials, biological materials are inherently structurally heterogeneous.Most of them consist of a matrix of structural biopolymers and various other organic components that are produced by the organisms themselves which control their formation and, either directly or indirectly, the formation of subassemblies which lead to the hierarchical organization of these materials.Among the known structural biopolymers, chitin is the second most abundant biopolymer in nature, aer cellulose, 1 and forms the organic matrix in the cell walls of fungi and the exoskeletons of Arthropoda. 2 Arthropod exoskeletons are formed from the cuticle, a hierarchically structured bre-based composite, where chitin is responsible for both shaping and reinforcing the material.The polysaccharide molecule chitin is a linear polymer of b-1,4 linked N-acetylglucosamine.3][4] Two other forms of chitin found in nature are b-chitin and g-chitin.b-Chitin can be found for example in squid pens. 5ts abundance is very low due to the parallel arrangement of chitin chains, which does not provide as strong a hydrogen bond pattern as in a-chitin.g-Chitin has been described in insect cocoons or other special structures formed by insects, 4 and it contains a repeating motif where two parallel chains are complemented by one in the opposite orientation.
Typical a-chitin crystallites occurring in the cuticle have a polygonal section prole and consist of 18-25 chitin molecules with a diameter between 2 and 5 nm and lengths of about 300 nm. 67][8][9] Potential chitinorganizing proteins currently under discussion are a class of evolutionarily conserved cuticle proteins featuring the so-called Rebers and Riddiford (R&R) domain 9 such as resilin, and the obstructor proteins 10 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. 11owever, the exact functions of these proteins have not been identied yet.Chitin-protein nanobrils 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 a-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 nanobrils.During the process of cuticle formation, 16 chitin molecules are produced by the membrane inserted enzyme chitin synthase. 17The enzyme catalyses the b-1-4 linkage between N-acetylglucosamines supplied from the cytoplasm.9][20][21] Here, the chitin and chitin-binding proteins present in the extracellular space organize into the typical nanobrils that form the exoskeleton, most likely in a self-assembly process.This process involves the assembly of chitin molecules into larger brils 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-specic structural differentiation in the nal tissue observed in many species of arthropods.
Even though crystalline a-chitin has been intensively studied in past years, 12,14,22,23 its structural properties remain unclear.Likewise experimental evidence concerning nanobril composition, structural parameters and hierarchical structure were never further expanded to explain the self-organization of chitin nanobrils.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 a-chitin nanobrils, 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 environment 24 and during adhesion to montmorillonite clay. 1 More complex studies try to describe adsorption of chitosan to chitin nanocrystallites, 25 thermodynamics of a-chitin de-crystallization, 26 and properties of chitin-protein composites. 27Recently, a coarse grain model of crystalline a-chitin was derived from experimental data and atomistic molecular dynamics simulations. 28owever, to our best knowledge, none of these studies did characterize self-assembly of a-chitin nanobrils and structures of possible intermediates in a comprehensive way.Due to methodological constraints, the self-assembly process of a-chitin nanobrils 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 nanobrils of various dimensions are.
Here, we evaluate the geometry of single chitin molecules and chitin nanobrils 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 nanobrils in terms of free energy changes.The studied nanobril models were chosen such as to reect successive intermediate bril congurations that presumably occur during the self-assembly process leading to the nanobril size reported in previous studies of arthropod cuticle. 29,30This 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.

Materials and methods
Chitin nanobril models were constructed in silico from the experimentally determined structure of an a-chitin crystal (Fig. 2). 12Each model is composed of oligosaccharide chains of constant length.Each chain contains twenty units of N-acetylglucosamine (GlcNAc) linked by b-1,4 glycosidic bonds capped with hydroxyl groups on both chain terminals.The chains are then packed in the a and b lattice directions of the source achitin crystal, which preserves native contacts among individual chains, including both hydrophobic interactions as well as hydrogen bonds.These lead to either a single chain (1000) or brils in an antiparallel arrangement as combinations of 3 (1110), 4 (2200), 6 (2220, 3300, 6000), 8 (4400), 9 (3330), 12 (6600), 19 (6760) and 34 (8998) chains.The number of non-zero digits in the name is the number of chitin chains in the ab plane, whereas the numeral at a given position is the amount of chains in the corresponding ac plane.Detailed information about the simulated systems is provided in Table S1 in ESI (ESI †).
Molecular dynamics simulations performed with the Amber 12 package 31 were employed to evaluate the structural and dynamic behavior of the prepared models.The chitin nano-brils were described with a GLYCAM 06g force eld. 32To mimic the native conditions in the chitin nanobril assembly zone, each structure was put in a truncated octahedral box, which was lled with water molecules described by the TIP3P water model. 33Long-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 xed 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 xed, enabling an integration time step of 2 fs to be used.Temperature was kept xed 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. 35The length of each production run was 200 ns, conguration snapshots were taken each 10 ps.
Trajectory analyses were performed in the ptraj and cpptraj modules 36 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. 38he number of contacts was calculated using a Fermi switch function: where d ij is the distance between the hydrogen H and acceptor A atoms of analyzed contact.The b and d 0 parameters were chosen as 10.0 ÅÀ1 and 2.5 Å.These parameters are able to discriminate hydrogen bonds from long-range non-specic interactions for a varied range of hydrogen and acceptor atoms, as can be seen from a comparison of the switch function and histograms of the distances d ij depicted in Fig. S1.† In our analysis, the directionality of a hydrogen bond is not taken explicitly into account.Thus in some cases, the obtained results might be contaminated by contacts that cannot be considered as strong hydrogen bonds, but rather electrostatic interactions.These situations are however rare as it is demonstrated in Fig. S2 † showing distribution of angle q(D-H/A) for selected hydrogen bonds.Majority of considered contacts have the angle higher than 110 , which is recommended minimum value for the hydrogen bond. 38We also compared our approach with the traditional one, which is usually based on simple true/false consideration based on fulllment of distance and angle criteria.Depending on the used criteria, our approach overestimates or underestimates the number of hydrogen bonds up to 25% or 5% for strict and weaker criteria, respectively (for details see Table S2 †).The axial chirality was described as the pseudo-dihedral angle a 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 nanobrils.The rst 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 j(C5 n , C4 n , O4 n , C1 n+1 ) and f(C4 n , O4 n , C1 n+1 , O5 n+1 ) dihedral angles.Rotamers of terminal CH-CH 2 OH groups were monitored using the dihedral angle d(C4, C5, C6, O6) and orientations of terminal CH 3 CONH-CH groups were monitored using the dihedral angle 3(C1, C2, N2, C2N).Ring conformations were analyzed by measuring the f and q ringpuckering coordinates for sugar rings (C1, C2, C3, C4, C5, and O5) and pseudo-rings (O4 n , C4 n , C3 n , O3 n , O5 n+1 , C1 n+1 ). 39he binding strength between chitin chains was evaluated by means of the MM/PBSA [40][41][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 nanobril).The absolute free energy G is composed of two major contributions: (a) the free energy of the solute in vacuum G vac for a solute conguration ensemble generated in the solvent and (b) the solvation free energy of the solute G sol , representing the vertical transfer from vacuum to solvent (thus this energy does not include any deformation energy).The rst term is calculated from the solute enthalpy H, which is taken as a time average of the total solute energies E MM calculated by means of molecular mechanics (MM), and entropy S. The calculation of the entropy term will be discussed later in more detail.
The solvation free energy is composed of two contributions: electrostatic G el sol and non-electrostatic G ne sol .The electrostatic part was evaluated by employing a Poisson-Boltzmann (PB) method 41 to describe the solvent implicitly.The PB method requires the denition of the atomic radii separating the interior and exterior of the system.In this work modied Bondi parameters 40 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 nonelectrostatic contribution was calculated by employing a surface tension model using the solvent accessible surface area 43,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: where G is the absolute free energy of the system containing n chitin chains, and G 1000 is the absolute free energy of the 1000 system.The same equation was employed for the calculation of individual contributions to the total reaction free energy.
The essential dynamics analysis (EDA) method 45,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 lter out thermal noise from analyzed trajectories.In the rst step, the covariance matrix was calculated from uctuations 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 uctuation 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 motionessential motion, whereas modes with smaller eigenvalues represent rapid motionthermal motion.New trajectories containing only correlated atomic motions individually projected from the source trajectory on ve eigenvectors with the largest eigenvalues were calculated together with projection factors p.
For further analysis, the trajectories were coarse-grained.These simplied 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 S t and rotational S r 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 S v was calculated for the same temperature using a quasi-harmonic approach. 47To 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 innite time. 48igures showing 3D structures of the studied systems were prepared in PyMol 49 and VMD. 50

Results and discussion
In the natural system, the pH and the ionic composition of the aqueous uid present in the ecdysial gap where chitin nano-brils assemble is controlled by the cells. 51Within our model environment, the experimentally determined low concentrations of Na + , K + , Mg 2+ , and Ca 2+ ions in the moulting uid of the crustacean Porcellio scaber during initial cuticle formation would have resulted in very few ions present in the simulation box.Since preliminary simulations in isotonic solution yielded similar results, all models were simulated in a water environment without the presence of any ions.Due to computational demands, the sizes of the models had to be kept reasonably small.The length of chitin chains was limited to twenty N-acetylglucosamine units, which corresponds to approximately 10.6 nm-long chains.This is about 30 times smaller than the reported length of chitin nanobrils occurring in nature.However, the models are long enough to avoid articial interactions between terminal sides in the axial direction.The average structures of our nanobril model systems obtained from the entire production MD trajectories are shown in Fig. 3.It is immediately obvious that an increasing number of chitin chains in the a and b directions have a strong impact on the resulting overall structure of individual nanobrils.
The overall exibility of the different nanobril models was monitored by root-mean square deviation (RMSD, Fig. 4 and S3 †).This property quanties 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 uctuation ranging from 3.3 to 15.1 Å, with a maximum at 5.3 Å.This indicates that the structure is rather exible.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 exible due to the lack of surrounding partners that would stabilize the structure.Surprisingly, it was found that some of the larger nanobrils, namely the 1110, 2200, 6000, 6600 systems, deviate even more from the crystalline conguration, indicating an even higher exibility.To explain this observation, the hydrogen bond network in nanobrils was analyzed, as this is thought to be one of the key interactions that inuence their shape, stability and mechanical properties.
The hydrogen bond analysis was performed for specic and non-specic pairs of hydrogen bond donor and acceptor atoms.The specic hydrogen bonds represent geometrically necessary interactions occurring in the crystalline a-chitin while the non-specic hydrogen bonds represent possible interactions that might coincidentally occur during its evolution or its deformation.The specic 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.
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 brils, 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 pseudorings, 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 singlechitin 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 quantied by the ring-puckering f and q angles (Fig. S4 and S5 †).The polar q angle is close to 180.0 for both kinds of rings, which correspond to the observed chair conformation.The azimuthal f angle shows different distributions for every simulated system, but as the polar q angle is always close to 180.0 , these changes actually reect only very small geometrical changes.The only signicant deviation of the polar q 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 signicantly 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 studies 52 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 specic 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 nanobril size and shape.One possible explanation for this might be that these hydrogen bonds are weaker than the previous ones and could thus be inuenced by the water environment at the surface of the nanobrils.To conrm 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 nanobrill w ) or when they are located inside the nanobril (l i ).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: hHBi ¼ n w l w + n i l i where n w and n i are the number of hydrogen bonds per given stacked layer edge that need to be determined for layers exposed to water or the nanobril internal environment, respectively.This simple model does not consider all contacts with water properly, for example on the ab and ac facets of the nanobrils.Nevertheless it performed very well, as the square of the Pearson correlation coefficient (R 2 ) between the number of hydrogen bonds obtained from the analyzed trajectories and that predicted by the model, found by the least-square tting approach, is in both cases higher than 0.99 (Table S3, Fig. S6 †).For the O2N-H6O hydrogen bond, the found parameters are 1.04 AE 0.35 and 3.92 AE 0.18 for n w and n i , respectively.For the O6-H60 hydrogen bond, they are 1.64 AE 0.61 and 5.57 AE 0.31, respectively.Indeed, these values conrm that the O2N-H6O and O6-H60 hydrogen bonds are weakened by the presence of the solvent, because in both cases the parameter n w is about a third of n i .
We also focused on the analysis of other possible hydrogen bonds that do not commonly occur in crystalline chitin.The aim was to nd possible correlations between the presence of such non-specic hydrogen bonds and the observed structural shape of the simulated nanobrils.However, the analysis of the O6-H3O, O6-H2N, O3-H6O, O3-H2N, O2N-H3O atom pairs did not conrm 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 nanobril.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 noncovalent (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 nanobril 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 specic hydrogen bonds.The evolution of the fourth specic 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-specic hydrogen bonds formed by the O6-H2N and O2N-H3O atom pairs was also monitored.Two changes occurred during the simulation: the rst one happened in the interval from the 100 th to 150 th ns and the second change, which was more rapid, at about the 425 th ns.These two changes split the simulations into three regions, for which representative structural snapshots are shown in Fig. S7.† During the rst 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-specic ones.From the specic 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 150 th to 400 th 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 specic 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-specic hydrogen bonds should be close to zero.At the end of the 1110 simulation, the nal 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 rst third.The further rearrangement must involve the structure refolding due to misaligned stacking.This was the reason why we did not continue with the simulation, because a signicantly longer simulation time would be required.
The other simulated systems did not undergo any signicant rearrangements like the one observed for the 1110 system.However, some of them still deviated signicantly 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 congurational chirality introduced by the stereogenic centers in the building saccharide units).As a measure of axial chirality, the pseudo-dihedral angle a was calculated and its distribution is shown in Fig. 6.Its average values, including the corresponding uctuations, are summarized in Table 2.
The factor that inuences the degree of axial chirality most signicantly is the number of stacked layers of chitin chains interacting in the b direction.The axial chirality decreases when a nanobril 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 a with different structural descriptors (Fig. S8-S10 and Tables S5 and S6 †).
The signicant correlation between the pseudo-dihedral angle a and structural descriptors was found for dihedral angles describing the geometry of the glycosidic bonds.The  dihedral angles j and f exhibit single Gaussian distributions (Fig. S9 †).The peak maxima (average values) (Table 2) then correlated well with the pseudo-dihedral angle a with R 2 of 0.97 and 0.99, respectively.The obtained linear models (Fig. S10 †) enable the extrapolation of j and f to a ¼ 0, which should correspond to the situation in crystalline a-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 j and f were not determined for the real structure of crystalline a-chitin, but extrapolated from the behavior of smaller chitin nanobrils, 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 eld employed for the description of chitin nanobrils.
We also checked if the observed agreement between our simulations, experimental, and ab initio data cannot be consequence of specic force-eld parameters that would force the system to adopt conguration with values of the j and f 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.
4][55][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 nanobril.Moreover it was shown that the phenomena is force-eld dependent. 54On timescale of our simulations ($200 ns) we did not observe such disappearance of the twist.Splitting of distribution of the pseudo-dihedral angle a 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 nanobrils 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 nanobril.Nevertheless they were not able to recognize that magnitude of the twist might exhibit anisotropy due to size of nanobrils along the a and b directions because their nano-brils were symmetrical.Also the molecular origin of the twist will be probably different as the authors identied 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 a and average values of the ring-puckering coordinate q 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 †).R 2 is 0.95 and 0.94, respectively.
The correlations found indicate that the twist distortion of the simulated nanobrils is a consequence of the internal axial chirality of individual chitin chains.The internal axial chirality is caused by deviations of the j and f 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 nanobrils 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 q puckering coordinate describing such pseudo-ring conformations.
The axial chirality of chitin nanobrils is also reected in key so-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 soness (e.g.how the energy of the system changes with the particular movement).Projections of the trajectories to the ve essential modes with highest (soest) 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 nonuniform distribution is caused by slow axial distortion, which is also visible in the distribution of the pseudo-dihedral angle a (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 nanobrils observed in arthropod exoskeletons, these properties might be important for the interaction between chitin nanobrils 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 so in the other directions, which might facilitate small structural adjustments that enable proper alignment to the protein binding sites in living systems.
The fact that the chitin brils in mature and developing exoskeleton are always conned by various chitin-binding proteins makes it difficult to obtain direct experimental evidence for axial chirality in the a-chitin nanobrils, since the interactions between chitin and proteins are likely to inuence the structural conformation of the chitin brils.However, synchrotron X-ray micro-diffraction experiments have shown that in arthropod exoskeletons the crystallographic axes a and b of the chitin brils 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 nanobrils which is retained despite the addition of the protein coating.
The self-assembly of chitin chains into chitin nanobrils 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 nanobril formation.
The changes in reaction enthalpy in the gas phase DH r /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 DH r /n decreases with an increasing number of chains to the limiting value (extrapolated for an innite number of chitin chains), which we estimated to be about À651.1 AE 58.4 kcal mol À1 .This behavior was expected, because with an increasing number of chitin chains the nanobril 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 nanobrils (negative enthalpy) we found that the solvent destabilizes their mutual binding.This is caused by the relatively strong interaction between nanobrils and solvent molecules, and the fact that the contact of individual disassembled chitin chains with the solvent is higher than in an assembled nanobril.The contribution of the solvent to the binding is expressed by the value of hDG sol,r i/n.The value of hDG sol,r i/n, which is the sum of hDG el,r i/n and 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.hDG ne,r i/n, increases up to the limiting value (Table S8, † Fig. 8), which we estimated to be about 394.7 AE 34.7 kcal mol À1 .This value is smaller than the limiting value for the reaction enthalpy in the gas phase, conrming 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 AE 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 nano-brils 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 nanobril 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 nanobril surface (the proportion of the chitin chains that remain solvated).The value of hDG sol,r i/n must then converge to the opposite value of the solvation energy of a single chitin chain.Indeed, the value of hDG sol,r i/n converges to 394.7 AE 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 nanobrils occurring in the natural chitin nanocomposites, which usually consist of 18-25 chitin molecules. 6It can be assumed that the evolution of chitin nanobrils with this geometry was at least partially driven by the fact that they are the smallest nanobrils 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 simplied 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 (DS r /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. 59Despite 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.

Conclusions
Employing a theoretical approach, we elucidate the structure, dynamics, energy characteristics and aspects of self-assembly of a single chitin chain and 10 chitin nanobrils composed of an increasing number of molecules in an aqueous environment.The molecular dynamics simulations of the smallest nanobril model, which includes three chitin chains, indicate that nano-bril formation is driven by a variety of weak non-specic hydrogen bonds during the early stages of nanobril formation.These non-specic hydrogen bonds, which are formed especially by the O6-H2N and O2N-H3O atom pairs, do not normally appear in crystalline chitin or assembled nanobrils.We can speculate that their presence might be important for the efficiency of nanobril self-assembly.Simulations of the other nine nanobrils show relatively stable structures.The most signicant discovery is axial chirality, which is exhibited by all nanobrils.The similar behavior was already described for cellulose nanobrils 55,56 but it seems that the molecular origin of this phenomenon in these two crystalline polysaccharides might be different.The largest axial distortion was observed for the 2220 and 6000 nanobril systems, where the twist along the axis of the 20 unit long chains is nearly 240 and 180 degrees, respectively.This axial chirality is a consequence of the intrinsic axial chirality of single chitin chains, which is altered by interactions between chains.The stacking of chains on top of each other in the a direction preserves their axial chirality, while the presence of additional nanobrils in the b direction signicantly decreases chirality.This can be explained by dipoledipole interactions between the acetamide groups from adjacent stacked chitin layers.These alter the mutual orientation of adjacent sugar rings along glycosidic bonds and thus inuence the overall axial chirality.Essential dynamics analysis revealed that axial twist distortion is the soest dynamic mode.This might be important for the formation of protein-coated chitin nanobrils in living organisms, where it could play a role in the binding of proteins.The calculated binding free energies reveal that the major driving force for binding is enthalpy, which is suppressed by desolvation and entropy.The progress of the total free reaction energy per single chitin chain reaches an energy plateau at about 20 chitin chains and above.This indicates that the 6760 and 8998 chitin nanobrils exhibit binding interactions that resemble those in crystalline chitin.The resulting excellent mechanical properties together with the rather high exibility of nanobrils in the a and b directions may very well be one reason why evolution favored chitin nanobrils with this geometry as basic components for biocomposites such as the exoskeleton of Arthropoda.

Fig. 1
Fig.1The eleven structural models of chitin nanofibrils analyzed in this study.Their denomination is based on the number of horizontally aligned chitin molecules (up to four), the numerical value of each digit indicates the number of vertically stacked molecules in the structure.

Fig. 2 (
Fig. 2 (Top) Schematic structure of repeating N-acetylglucosamine dimers with atom names highlighted in blue.(Bottom) 3D structure of repeating N-acetylglucosamine dimers in a chitin chain, axes shown correspond to crystallographic axes.

Fig. 3
Fig. 3 Average structures of simulated nanofibril systems.Building blocks are color-coded to highlight structural features.

Fig. 5
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.

Fig. 6
Fig. 6 Distribution of pseudo-dihedral angle a quantifying axial chirality of nanofibril chains.

Fig. 8
Fig. 8 Dependence of change in reaction enthalpy DH r , electrostatic DG el,r and non-electrostatic DG ne,r contributions to solvation free reaction energies and total free reaction energy DG r per single chitin chain on the number of chitin chains n (excluding the entropy contributions).Lines are meant to only highlight the progress of individual energy contributions.
project CEITEC 2020 (LQ1601) and the project CZ.1.05/1.1.00/02.0068 nanced from the European Regional Development Fund.The computational resources were provided by the CES-NET LM2015042 and the CERIT Scientic Cloud LM2015085, provided under the programme "Projects of Large Research, Development, and Innovations Infrastructures".This work was also supported by the IT4Innovations Centre of Excellence project (CZ.1.05/1.1.00/02.0070,individual projects IT4I-9-11, OPEN-4-2, OPEN-5-10, OPEN-6-10 and OPEN-6-23), funded by the European Regional Development Fund and the national budget of the Czech Republic via the Research and Development for Innovations Operational Program, and the Ministry of Education, Youth and Sports of the Czech Republic via the programme "Projects of Large Research, Development, and Innovations Infrastructures" (LM2011033).The parts of this study that were conducted at the Max-Planck-Institut für Eisenforschung, Düsseldorf, Germany benetted from nancial support from the German Research Foundation (DFG) within the priority programme SPP 1420.M. F. acknowledges nancial support from the Academy of Sciences of the Czech Republic through the Fellowship of Jan Evangelista Purkynĕ.Z. S. is grateful for nancial support from the Max Planck Society during her stay at the Max Planck Institute for Iron Research in Düsseldorf, Germany.

Table 1
Statistics of specific hydrogen bonds a hHBi s(HB) hHBi/n hHBi s(HB) hHBi/la nnumber of chitin chains; lnumber of stacking layers; hHBiaverage number of hydrogen bonds; s(HB)variance (uctuation) of number of hydrogen bonds; hHBi/naverage number of hydrogen bonds related to a single chitin chain; hHBi/laverage number of hydrogen bonds related to a single stacked layer of chitins.

Table 2
Average values of dihedral angles a, j and f and their root mean square fluctuations s

Table 3
Comparison of dihedral angles j and f