Effect of Acetylation on the Nanofibril Formation of Chitosan

1. Introduction

Chitosan is a linear copolymer of D-glucosamine (GlcN) and N-D-acetylglucosamine (GlcNAc) linked by $\beta$-1,4-glycosidic bonds. Owing to the solubility in aqueous solution, pH responsiveness, and antimicrobial properties [1], chitosan has found vast applications in biotechnology and pharmaceutics [2,3]. Chitosan is produced by deacetylation of chitin, nature’s second most abundant biopolymer, which is found in the shells of crustaceans and fungal cell walls [1]. The deacetylation is never complete, resulting in chains with varying distributions of acetyl groups or degree of acetylation (DA) [4,5]. Additionally, depending on the method of deacetylation, the relative arrangement of GlcN and GlcNAc units or the pattern of acetylation (PA) may be different [4,5]. Generally, chitosan can be prepared chemically by treatment of chitin with a strong base which results in a random PA [4] or enzymatically via the action of chitin deacetylases which yields non-random PA [6,7].

In analogy to proteins, the DA and PA along with the degree of polymerisation (DP) are the primary parameters that describe chitosan’s sequence [5,8] (Figure 1). Recent experimental work demonstrated that the DP, DA, and PA influence the physicochemical properties and biological activities of chitosan [4,5,8,9]. However, the sequence-property or sequence-activity relationship of chitosan is poorly understood, which hinders the rational design and development of chitosan-based functional materials for novel applications [8]. Historically, this was a result of limited control in the synthesis of chitosan with specific DA, PA, and DP [5]. Improvements in synthesis and characterisation methods [5] have allowed for the control of DA and DP but control of PA remains a challenge [6,8].

Several experimental studies have examined the effect of DA on the physicochemical, materials, and biological properties of chitosan. Foster et al. found that increasing DA within the range of 15–28% is associated with the decreasing surface roughness, tensile strength, crystallinity of the chitosan films [10]. These results are consistent with an earlier study by Cao et al. [11], and may be related to the observation [12] that higher DA leads to lower cell adhesion and lower cell proliferation on the chitosan films. More recently, Moerschbacher and coworkers examined the effect of PA on chitosan’s solution properties and gelation behavior using a variety of experimental techniques including ¹³C-NMR, rheology measurements, and circular dichroism spectroscopy [8]. They found that, in stark contrast to the commercially available random-PA chitosan, the block-PA chitosan produced with chitin deacetylases had lower intrinsic viscosity in water, limited particle formation, and weak gelation behavior [8].

Despite the experimental efforts to investigate the sequence-property relationship of chitosan, very few computational or theoretical work has been published so far. Franca et al. [13] performed molecular dynamics (MD) simulations of a prebuilt chitosan nanoparticle (or nanofibril) comprised of 9 (3 by 3 matrix) antiparallel chains each having 10 monomer units in the block or alternating PA and different DAs. The results suggested that the chain solubility is inversely related to the DA and the chitosan nanofibril with the block PA is more stable [13]; however, these conclusions were based on an extremely short (20 ns production run) MD trajectory of a minimally sized nanofibril. Importantly, the simulations started from a nanoparticle, which precludes the investigation of chitosan’s self-assembly behavior. More recently, Tsereteli and Grafmüller [14] developed a coarse-grained model for titratable chitosan and validated the model using Metropolis Monte-Carlo simulations of a single chitosan chain with various DP (656–2100), DA (5–100%) and PA at different solution pH conditions. Based on the radius of gyration of the single chitosan chain, the study [14] suggested that the increased DA leads to increased compaction or chain flexibility and the chain with the alternating PA is more compact than the block-PA chain. Although these results are interesting and novel, the simulations did not account for chain-chain association, which is critical for understanding the sequence-property relationship of chitosan. Furthermore, due to the coarse-grained nature of the model, atomic details cannot be gleaned.

The objective of the present work is to understand at the atomic level the effects of the DA and PA on chitosan’s self-assembled microscopic structure. MD simulations rely on the accurate parameters in the molecular mechanics force field or energy function to describe the dynamics of molecules. Thus, we first refined the parameters in the current all-atom CHARMM molecular mechanics force field of chitosan [15,16]. We then conducted, to our best knowledge, the first de novo all-atom MD simulations of chitosan chains with different DAs and PAs from the initially randomly distributed configurations in solution. The simulations allowed us to examine, to our best knowledge, for the first time, how the DA and PA influence the formation and the structure of the nanofibril. Comparison with the relevant experimental data is discussed. The present work provides an important step towards the atomic-level understanding of how acetylation modulates chitosan’s physicochemical properties and contributes to the rational design of chitosan-based materials with tunable degree and pattern of acetylation. The improved force field parameters are disseminated to the community to facilitate future studies of chitosan-based materials.

2. Results and Discussion

2.1. Optimization of the molecular mechanics force field and summary of simulations

The CHARMM36 carbohydrate force field [15,16] already contains optimized parameters for a subset of $\beta$-linked carbohydrates, including $\beta$-acetylglucosamine (GlcNAc), the parent molecule of $\beta$-glucosamine (GlcN). These two molecules differ only in two torsion angle parameters related to the primary amino group (Figure 1). In our previous studies of chitosan [17,18], the aliphatic primary amine parameters were adopted according to the CGENFF force field [19]; however, we found that the potential energy profiles of the two torsional angles deviate significantly from the quantum mechanical (QM) calculations (Figure 2). Thus, following the protocol used in the development of the CHARMM36 carbohydrate force field [15,16], we manually adjusted the torsion parameters until they matched the QM data (Figure 2b and d). All subsequent simulations were performed with these optimized parameters which are made freely available (see Data Availability).

The chitosan de novo self-assembly simulations were initiated with 24 chains randomly distributed in a 85 Åx85 Åx85 Å cubic box filled with water and periodically replicated in three dimensions to give the weight percent concentrations of 8.6% and 9.2% for the 20% and 50% DA systems respectively. Each chain is comprised of 10 monomer GlcN or GlcNAc units, whereby the DA was 20% or 50% for all chains. At 50% DA, two different patterns of acetylation, alternating and block, were used. For each combination of DA or PA, triplicate runs were conducted starting from the same initial arrangements of the 24 chains but different initial velocities for atoms (velocity seeds). All trajectories were continued until the total number of interchain hydrogen bonds (h-bonds) and total solvent accessible surface area (SASA) plateaued, resulting in ∼2 $\mu$s for DA of 20% and ∼3-4.5 $\mu$s for DA of 50% runs (Table 1). Note, since the present work is focused on the effect of DA and PA on the self-assembled structure, the neutral glucosamine was used for the deacetylated units, i.e., high pH condition. Thus, the electrostatic repulsion between the protonated GlcN residues was not considered.

2.2. Chitosan solubility and self-assembly is dependent on the degree of acetylation

We first investigated the effect of DA using simulations of the 20% and 50% acetylated chitosan chains. We found that at 20% DA, no self-assembly occurred and the chitosan chains remained in solution after 2 $\mu$s, with only transient formation of intermolecular hydrogen bonds (h-bonds, Supplemental Figure S1). In contrast, at 50% DA, chain association via the formation of intermolecular h-bonds occurred rapidly after 1 $\mu$s of simulation (Supplemental Figure S1). This difference suggests that acetylation promotes chain-chain association, consistent with an early experimental study which showed that at above 50% DA even the protonated chitosan chains can aggregate [21]. Note, the observation that the 20%-DA chitosan chains did not self assemble within the simulation time reflects the short chain length (DP=10). With a long chain length (or increased DP), chitosan chains can form a larger number of intermolecular h-bonds (see later discussion), and thus self assembly can occur at lower DA, e.g., the nanofibrils of fully deacetylated 20-mer chitosan chains are stable (Romany and Shen, unpublished data). Thus, the effect of DA on self assembly may be amplified in our simulations.

To characterize solubility, we calculated the percent solvent exposure of the chitosan chains defined as the solvent accessible surface area (SASA) relative to that of a fully solvated chain. The data for chitosan chains with 20% and 50% DAs were compared. We also included a similar calculation of chitin, which is chitosan with 100% DA (taken from our previous study [22]). The probability distributions show that the % solvent exposure decreases with DA (peak of 90%, Figure 3a). At DA 50% and 100%, the percent solvent exposure per chain is below 60%, whereas at DA at 20% the solvent exposure per chain is near 100% (peak of 90%, Figure 3a). This can be explained, as water is excluded around the self-assembled chains with DA above 50%, whereas the unassembled chains with 20% DA are surrounded by solvent.

To understand the physical basis for the effect of DA on the solubility and self-assembly of chitosan, we examined the h-bond formation between water and GlcNAc or GlcN. We found that the primary amine nitrogen of GlcN is able to accept h-bonds from water, but the acetylated nitrogen of GlcNAc does not form h-bonds with water (Figure 3b, purple bar). The same behavior is observed for chitosan chains with both block and alternating PAs (Figure 3b, purple and blue bars). Since the lower-DA chitosan chains have higher amine nitrogen content as compared to the higher-DA chains, they form more h-bonds with water, which makes the solvated form energetically more stable. We suggest the difference in the h-bond formation with water as a major contributor to the solubility increase with decreasing DA. The amine–water h-bond formation is in direct competition with the chain-chain (NH–O6) h-bond formation which contributes to self-assembly. Thus, the h-bond formation between the primary amine and water also contributes to the decrease in self assembly with the decreased DA.

We next examined the effect of DA on the flexibility of chitosan chain by calculating the distribution of the end-to-end distances of the chains in solution. For the 20%-DA chains, the distribution is slightly wider and the speak intensity is lower than the 50%-DA chains (Figure 3c), which indicates that the lower-DA chains are more flexible compared to the higher-DA chains in solution. This data is consistent with the persistent length data from the experiments [21,23] but disagrees with the coarse-grained simulations which suggested that increased DA leads to increased chain flexibility [14]. This data supports the notion that the bulkier acetamido groups (as compared to the amino group) limit the chitosan chain rotation such that the chain stiffness increases with DA [21].

2.3. Pattern of acetylation influences the chitosan’s self assembly process.

We examined and compared the self-assembly process of chitosan at 50% DA with the block and alternating PAs. Interestingly, regardless of the PA, at the end of the simulation runs when the number of interchain h-bonds and solvent exposure plateaus, the chitosan chains formed a single nanofibril with at most 2 chains left in solution (Figure 4a and b). Regardless of the PA, the percentage of solvent exposure and the degree of order of the self-assembled nanofibrils are also similar (Figure 4a and b). The latter was measured by the $P_2$ order parameter (Equation 1) as in our previous work [22,24]. However, according to the number of interchain h-bonds and the total relative solvent exposure, the self-assembly process converged much faster (within ∼2 $\mu$s) for the block PA chains in all three runs as compared to the all three runs of the alternating PA chains (more than 3 $\mu$s, Supplemental Figure S2). The recent experiment of Moerschbacher and coworkers demonstrated that the gelation kinetics of the biotech chitosan that contains block PA is distinct from the conventional random-PA chitosan [8]. Specifically, they found that in contrast to the random-PA chitosan, the block-PA chitosan did not show a single critical gelation time [8]. Although the difference in the kinetics of self-assembly simulations cannot be directly compared to the difference observed in the experiment [8], they may be related. Moerschbacher and coworkers also found that the block-PA chitosan nanoparticles showed a higher dispersity in size than random-PA nanoparticles, which was attributed to the formation of small aggregates [8]. The more rapid self-assembly of block-PA chitosan revealed by the simulations is consistent with this observation.

To quantify the self-assembly process we used a low dimensional free energy surface (FES) spanned by the radius of gyration ($R_g$) of all the chitosan chains and the number of interchain h-bonds per monomer unit (Figure 4c and d), following our previous work [22]. When the chitosan chains are in solution (disassembled state), the radius of gyration of all chains is high and the number of interchain h-bonds is very low, while in the fully assembled state, the radius of gyration is low (the chains are associated) and the number of h-bonds is maximized. There are two extreme cases for the shape of the FES. A diagonally-shaped FES suggests that the chain-chain hydrophobic association and h-bond formation are concomitant, whereas an L-shaped FES indicates that hydrophobic collapse proceeds the h-bond formation [22]. Our recent simulations of chitin’s self-assembly process [22] demonstrated that chitin’s fibril formation can follow either diagonally- or L-shaped mechanism. Here, the chitosan chains with both block and alternating PAs appear to more closely follow the diagonal mechanism, in which the hydrophobic association along with water expulsion occur concomitantly as the interchain h-bond formation (Figure 4c and d; Supplemental Figure S3). Interestingly, the free energy difference between the maximum (disassembled state) and global minimum (assembled state) is about 3 or 4 kcal/mol, similar to that found for the self-assembly of chitin [22]. It is interesting to note that, regardless of the PA, the increased of the number of interchain h-bonds is correlated with the decrease of solvent exposure, consistent with the self-assembly simulations of chitin [22]. Nonetheless, the self-assembled chitosan nanofibrils with either block or alternating PA are hydrated, with sheet-sheet stacking stabilized by water-mediated interactions (Figure 4e), similar to the self-assembled chitin nanofibrils seen in our previous simulations [22].

2.4. The self-assembled nanofibrils are comprised of nearly exclusively antiparallel chains regardless of the acetylation pattern

Experimental studies of chitosan (presumably with a very low degree of acetylation) demonstrated several different crystalline polymorphs, with the hydrated forms (often referred to as tendon) being more abundant than the anhydrous forms (referred to as annealed), which are created by heating or exposure to carboxylic acid [25]. The X-ray fiber diffraction measurements of the hydrated chitosan crystals [26,27] suggested that the fiber sheets are formed by antiparallel chitosan chains. Fiber sheets comprised of antiparallel chains were also suggested by the X-ray [28] and electron diffraction [29] studies of the anhydrous chitosan crystals prepared by high temperature. This arrangement was further confirmed by a subsequent study where the anhydrous chitosan crystals were prepared by using acetic acid [30] as well as a more recent high-resolution X ray fiber diffraction measurement [25].

To examine the effect of PA on chitosan’s polymorphism, we analyzed the relative orientations of chitosan chains in the self-assembled nanofibrils with the block and alternating PAs. Note, the nanofibrils are hydrated. The chain orientation was characterized by calculating the angles between the adjacent chains. Angles of near 0° and 180° indicate parallel and antiparallel chains, respectively. The probability distributions of the angles show that the antiparallel orientation is nearly exclusive for the nanofibrils with either block or alternating PA (Figure 5a-c), which is in agreement with the aforementioned X-ray diffraction results for the hydrated chitosan crystals [25,28,29,30].

Importantly, the simulations of partially acetylated chitosan is in contrast to the recent simulations of the fully acetylated chitosan (chitin) [22], which resulted in nanofibrils comprised of both parallel and antiparallel chains, and the antiparallel orientation is more populated at high temperatures. In nature, chitin exists in three polymorphs: the most common $\alpha$-chitin contains antiparallel chains [31], while $\beta$-chitin contains parallel chains [32] and the rare $\gamma$-chitin contains mixed antiparallel and parallel chains [33]. Thus, both our present and previous simulations [22] are consistent with the experimental observations and together they suggest that deactylation promotes the antiparallel arrangement of chains.

2.5. Intermolecular hydrogen bonding pattern is similar for the nanofibrils with the block and alternating PAs.

To examine the intermolecular h-bonds that drive the self-assembly of chitosan and the potential effect of PA, we calculated the occupancies of specific h-bonds formed between the antiparallel chains for the nanofibrils with the block and alternating PAs (Figure 6a). Interestingly, the nanofibrils with both PAs show the same intermolecular h-bonds, between O6H and O7, between NH and O6, and between O3H and O7 (Figure 6a). This is consistent with our previous simulations of chitin [22], which found the same three h-bonds in the self-assembled chitin nanofibrils comprised of antiparallel chains, whereas a different NH–O7 h-bond was formed between the parallel chains. It is worthwhile noticing that the O6H–O7 and NH–O6 h-bonds are dominant in the chitosan nanofibrils with the block PA (Figure 6a), which is consistent with the chitin nanofibrils of antiparallel alignment [22]. In contrast, the O6H–O7 h-bond dominates the chitosan nanofibrils with alternating PA (Figure 6a). This analysis suggests that chitosan nanofibrils with the block PA is more similar to the chitin nanofibrils with antiparallel alignment (i.e., $\alpha$-chitin). The X-ray or electron diffraction measurements of both hydrated [26,27] and anhydrous chitosan crystals [26,27] suggested that the chains on the fibril sheets are antiparallel and they are joined together by the NH–O6 h-bonds. Considering that the NH–O6 h-bonds have a higher relative occupancy in the simulated block-PA vs. alternating-PA nanofibrils, we suggest that the block-PA chitosan is more similar to the fully deacetylated chitosan.

To further understand how PA influences the h-bonding pattern of the self-assembled chitosan nanofibrils, we examined the formation of the two dominant intermolecular h-bonds for the blockNH–O6 and O6H–O7 h-bonds between two acetylated (GlcNAc–GlcNAc), two deacetylated (GlcN–GlcN), and mixed acetylated and deacetylated (GlcNAc–GlcN) glucosamine units (Figure 6b and c). In the nanofibrils with the block PA, the majority (70%) of the NH–O6 h-bonds are between GlcNAc and GlcNAc, followed by 25% between GlcNAc and GlcN and 5% between GlcN and GlcN units (Figure 6b left and Figure 6c). The finding that the h-bonds prefer to form between the acetylated units was initially surprising, as it is unlikely for in-register antiparallel chitosan chains with the block PA. However, it is consistent with our observation that the chitosan chains with block PA tend to be off register in the nanofibrils (Figure 5c). In contrast and as expected, for the nanofibrils with the alternating PA, the majority (∼60%) of the NH–O6 h-bonds are formed between GlcNAc and GlcN units. Nonetheless, 30% of the NH–O6 h-bonds are formed between GlcNAc and GlcNAc and 10% between GlcN and GlcN units, which can be partly attributed to the slightly off-register alignment between the two antiparallel chains in the nanofibrils with alternating PA (Figure 5c). Similar to the intermolecular NH–O6 h-bonds, the intermolecular O6H–O7 h-bonds also prefer to form between GlcNAc and GlcNAc units in the nanofibrils with block PA and between GlcNAc and GlcN units in the nanofibrils with alternating PA (Figure 6b right and Figure 6c). Taken together, the analysis of the h-bonding suggests that glucosamine acetylation plays a central role in stabilizing the nanofibril formation. We note that with longer chains (i.e., higher DP), the occupancies of the intermolecular h-bonds, including the NH–O6 h-bonds between two GlcN groups, are expected to be much higher due to the more stable nanofibril formation (Romany and Shen, unpublished data); however, it would not affect our conclusion that acetylation strengthens the NH–O6 h-bonds, which is related to the increased solvation of amino vs. acetamide group (Figure 3b).

3. Concluding Discussion

The MD simulations employing the newly optimized force field parameters for $\beta$-glucosamine were conducted de novo to investigate the self-assembly of 10-mer chitosan chains with different DAs and PAs. Consistent with the experimental solubility data [21], the MD simulations showed that the solvent accessibility of chitosan chains decreases with increasing DA. At 50% DA and with either block or alternating PA, the chitosan chains associated to form hydrated nanofibrils comprised of nearly exclusively antiparallel chains, in agreement with the fiber diffraction measurements of the hydrated chitosan crystals [26,27]. Since our previous simulations and experiments [31] demonstrated that chitin nanofibrils can display both antiparallel and parallel chain arrangements, the new simulation results suggest that acetylation promotes the parallel chain arrangements.

Further analysis of the self-assembled nanofibrils of 50% DA chitosan found that the nanofibrils of block and alternating PAs share three common intermolecular h-bonds, O6H–O7, NH–O6, and O3H–O6, with the first two being most stable. Interestingly, for the block PA, the intermolecular h-bonds mostly occur between two acetylated units, which was initially surprising given the antiparallel chain arrangements; however, this can be rationalized by the off-register alignment of two antiparallel chains with the block PA. In contrast, for the alternating pattern, the intermolecular h-bonds most frequently occur between one deacetylated and one acetylated unit, which is expected from the acetylation pattern. Another related finding from the MD simulations is that the NH–O6 h-bond between two deacetylated units is significantly weaker than between two acetylated units or between one acetylated and one deacetylated unit. This difference suggests that acetylation promotes self assembly, which is consistent with the decreased solubility of chitosan with increased DA [21]. The latter can be attributed to the decreased solvation of acetamide relative to (neutral) amino group.

A main caveat of the current study is use of very short chitosan chains (10 glucosamine units), which weakens the chain-chain association and the intermolecular h-bonds. For example, the NH–O6 h-bonds are unstable here but very stable between the 20-mer chitosan chains (Romany and Shen, unpublished data). Another caveat is the limited simulation time, which may result in the underestimated order for the self-assembled nanofibrils. The extent of off registry for the block-PA nanofibrils may be exaggerated due to the limited sampling time as well as the short chain length. Nonetheless, the comparison between the block-PA and alternating-PA nanofibrils is robust, given similar simulation time and multiple trajectories. Taken together, our simulations provide insights on the dynamics and mechanisms of self-assembly of chitosan that are difficult to observe by experiments. The analysis of the nanofibril structures and h-bonding patterns deepens the understanding of how PA influences chitosan’s physicochemical properties. We demonstrated that regardless of the PA, chitosan’s remarkable self-assembly is still observed. With the growing demand for biodegradable and sustainable materials, understanding the sequence-structure relationship is a key to the rational design of novel chitosan-based materials.

4. Materials and Methods

4.1. Force field parameterization

The starting parameters for GlcN were based on the sugar parameters available for GlcNAc ($\beta$(1-4)acetyl-glucosamine) and the modified CGenFF methyl group parameters [19]. The parameters for two dihedral angles related to the amine group, C1–C2–N–H11 and H2–C2–N–H12, were optimized (Figure 2). Staring from the diffraction structure of Okuyama and Ogawa [26], the model compound (Figure 2) first underwent geometry optimization using the quantum mechanical (QM) calculations at the MP2 level of theory and 6-31G(d) basis set. The QM potential energy surface (PES) for the H2-C2-N-H12 and C1-C2-N-H11 were then calculated at the MP2 level of theory and 6-31G(d) basis set. The parameters were manually tuned until the QM PES and the molecular mechanics (MM) PES were in good agreement. All QM calculations were performed using Gaussian-03 [34]. Parameterization was facilitated using the FFparam software [35].

4.2. Molecular dynamics simulation

All-atom molecular dynamics (MD) simulations were performed using the AMBER20 program [36]. The chitosan chains were represented by the CHARMM36 carbohydrate force field [15,16] and the modified CHARMM General force field parameters [19] (see Section Optimization of the molecular mechanics force field). Different acetylation patterns and degrees of acetylation were produced by modifying a 10-mer chitin chain using Avogadro molecular editor [37]. For each distinct system PACKMOL [38] was used to randomly distribute 24 chains in a 85 Å x 85 Å x 85 Å cubic box solvated with CHARMM style TIP3P water model[39]. The simulation systems were comprised of approximately 92493 atoms and the weight percentage of chitosan was 8.6% for 20% DA and 92493 atoms with a weight percentage of 9.2% for 50% DA.

A total of 9 simulations were conducted of the system. Three independent simulations were performed for each system starting from different random velocity seeds. The systems first underwent energy minimization with a harmonic restraint potential (force constant of 10 kcal/mol/A) placed on the heavy atom positions. The system was then heated over 1 ns to 300K in the NVT ensemble. The system was then equilibrated for a total of 100 ns in which harmonic potential force constant was gradually reduced from 5, 2.5, 1, 0.1 to 0 kcal/mol/A under constant NPT conditions. The temperature was maintained using a Langevin thermostat[40] and the pressure was maintained at 1 bar using the Monte Carlo barostat.[41] Each production run lasted between 3 $\mu$s and 4.5 $\mu$s under the constant NPT condition (Table 1). The van der Waals interactions were smoothly switched to zero from 10 to 12 A. The particle mesh Ewald method was used to calculate long-range electrostatic energies with a sixth-order interpolation and 1 Å grid spacing. Bonds involving hydrogens were constrained using the SHAKE algorithm [42] which allows the use of a 2-fs time step for integration. Subsequent analysis of simulations were done using cpptraj [43] and VMD [44].

The P₂ order parameter was defined as [24]

$$P_2=\frac{1}{N}\sum _{i=1}^N{({\overrightarrow{v}}_i.\overrightarrow{d})}^2,$$

(1)

where ${\overrightarrow{v}}_i$ is the normalized molecular vector connecting the center of mass of the second and ninth sugar rings in an individual chain. $\overrightarrow{d}$ is the director, which is a normalized average vector of the unsigned orientations of all individual ${\overrightarrow{v}}_i$. and N is the number of molecular vectors (number of chitosan chains).