Prerpints.org logo
Preprint
Article

Effect of Acetylation on the Nanofibril Formation of Chitosan

Altmetrics

Downloads

153

Views

95

Comments

0

A peer-reviewed article of this preprint also exists.

supplementary.pdf (645.19KB )

This version is not peer-reviewed

Submitted:

15 December 2023

Posted:

18 December 2023

Read the latest preprint version here

Alerts
Abstract
Chitosan is a copolymer of β-glucosamine and N-acetylglucosamine. Recent experiments showed that the degree and pattern of acetylation along the chitosan chain modulate the biological and physicochemical properties; however, the molecular mechanism is unknown. Here we apply the de novo all-atom molecular dynamics (MD) simulations to study chitosan’s self-assembly process at different degrees of acetylation and two different acetylation patterns. The simulations revealed that the 10-mer chitosan chains with 50% acetylation in either block or alternating patterns associate to form ordered nanofibrils comprised of mainly antiparallel chains in agreement with the fiber diffraction data of deacetylated chitosan. Surprisingly, regardless of the acetylation pattern, the same intermolecular hydrogen bonds mediate the fibril sheet formation while water-mediated interactions stabilize the sheet-sheet stacking. Moreover, the acetylated units are involved in forming the strong intermolecular hydrogen bonds (NH–O6 and O6H–O7), which offers an explanation for the experimental observation that increased acetylation lowers chitosan’s solubility. Taken together, the present study provides the atomic-level understanding the role of acetylation plays in modulating chitosan’s physiochemical properties, contributing to the rational design of chitosan-based materials with the ability to tune by its degree and pattern of acetylation. Additionally, we disseminate the improved molecular mechanics parameters that can be applied in MD studies to further understand chitosan-based materials.
Keywords: 
Subject: Chemistry and Materials Science  -   Biomaterials

1. Introduction

Chitosan is a linear copolymer of D-glucosamine (GlcN) and N-D-acetylglucosamine (GlcNAc) linked by β -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 13C-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 β -linked carbohydrates, including β -acetylglucosamine (GlcNAc), the parent molecule of β -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 μ s for DA of 20% and ∼3-4.5 μ 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 μ 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 μ 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 μ 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 μ 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 α -chitin contains antiparallel chains [31], while β -chitin contains parallel chains [32] and the rare γ -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., α -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 β -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 ( β (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 μ s and 4.5 μ 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 P2 order parameter was defined as [24]
P 2 = 1 N ∑ i = 1 N ( v → i . d → ) 2 ,
where v → i is the normalized molecular vector connecting the center of mass of the second and ninth sugar rings in an individual chain. d → is the director, which is a normalized average vector of the unsigned orientations of all individual v → i . and N is the number of molecular vectors (number of chitosan chains).

Supplementary Materials

Supplemental figures contain additional analysis of the data.

Data Availability Statement

All simulation input files, analysis scripts, and optimized force field parameters are freely available for download at https://github.com/janashenlab/chitosan.

Acknowledgments

Financial support from the National Science Foundation (CBET1932963) is acknowledged.

References

  1. Rinaudo, M. Chitin and Chitosan: Properties and Applcations. Progress in Polymer Science. Prog. Polym. Sci. 2006, 31, 603–632. [Google Scholar] [CrossRef]
  2. Kumar, M.N.V.R.; Muzzarelli, R.A.A.; Muzzarelli, C.; Sashiwa, H.; Domb, A.J. Chitosan Chemistry and Pharmaceutical Perspectives. Chem. Rev. 2004, 104, 6017–6084. [Google Scholar] [CrossRef]
  3. Jimenez-Gomez, C.P.; Cecilia, J.A. Chitosan: A Natural Biopolymer with a Wide and Varied Range of Applications. molecules 2020, 25, 1–43. [Google Scholar] [CrossRef] [PubMed]
  4. Weinhold, M.X.; Sauvageau, J.C.; Kumirska, J.; Thoming, J. Studies on Acetylation Patterns of Different Chitosan Preparations. Carbohydr. Polym. 2009, 78, 678–684. [Google Scholar] [CrossRef]
  5. Wattjes, J.; Sreekumar, S.; Richter, C.; Cord-Landwehr, S.; Singh, R.; El Geuddari, N.; Moerschbacher, B.M. Patterns Matter Part 1: Chitosan Polymers with Non-Random Patterns of Acetylation. React. Funct. Polym. 2020, 151, 104583. [Google Scholar] [CrossRef]
  6. Wattjes, J.; Niehues, A.; Cord-Landwehr, S.; Hoßbach, J.; David, L.; Delair, T.; Moerschbacher, B.M. Enzymatic Production and Enzymatic-Mass Spectrometric Fingerprinting Analysis of Chitosan Polymers with Different Nonrandom Patterns of Acetylation. J. Am. Chem. Soc. 2019, 141, 3137–3145. [Google Scholar] [CrossRef] [PubMed]
  7. Kaczmarek, M.B.; Struszczyk-Swita, K.; Li, X.; Szczęsna-Antczak, M.; Daroch, M. Enzymatic Modifications of Chitin, Chitosan, and Chitooligosaccharides. Front. Bioeng. Biotechnol. 2019, 7. [Google Scholar] [CrossRef] [PubMed]
  8. Sreekumar, S.; Wattjes, J.; Niehues, A.; Mengoni, T.; Mendes, A.C.; Morris, E.R.; Goycoolea, F.M.; Moerschbacher, B.M. Biotechnologically Produced Chitosans with Nonrandom Acetylation Patterns Differ from Conventional Chitosans in Properties and Activities. Nat. Commun. 2022, 13, 7125. [Google Scholar] [CrossRef] [PubMed]
  9. Basa, S.; Nampally, M.; Honorato, T.; Das, S.N.; Podile, A.R.; El Gueddari, N.E.; Moerschbacher, B.M. The Pattern of Acetylation Defines the Priming Activity of Chitosan Tetramers. J. Am. Chem. Soc. 2020, 142, 1975–1986. [Google Scholar] [CrossRef]
  10. Foster, L.J.R.; Ho, S.; Hook, J.; Basuki, M.; Marçal, H. Chitosan as a Biomaterial: Influence of Degree of Deacetylation on Its Physiochemical, Material and Biological Properties. PLOS ONE 2015, 10, e0135153. [Google Scholar] [CrossRef]
  11. Cao, W.; Jing, D.; Li, J.; Gong, Y.; Zhao, N.; Zhang, Z. Effects of the Degree of Deacetylation on the Physicochemical Properties and Schwann Cell Affinity of Chitosan Films. J. Biomater. Appl. 2005, 20, 157–177. [Google Scholar] [CrossRef] [PubMed]
  12. Chatelet, C.; Damour, O.; Domard, A. Influence of the Degree of Acetylation on Some Biological Properties of Chitosan "lms. Biomaterials 2001, 22, 261–268. [Google Scholar] [CrossRef] [PubMed]
  13. Franca, E.F.; Freitas, L.C.G.; Lins, R.D. Chitosan Molecular Structure as a Function of N-acetylation. Biopolymers 2011, 95, 448–460. [Google Scholar] [CrossRef]
  14. Tsereteli, L.; Grafmüller, A. An Accurate Coarse-Grained Model for Chitosan Polysaccharides in Aqueous Solution. PLoS ONE 2017, 12, e0180938. [Google Scholar] [CrossRef]
  15. Guvench, O.; Hatcher, E.; Venable, R.M.; Pastor, R.W.; MacKerell, A.D. CHARMM Additive All-Atom Force Field for Glycosidic Linkages between Hexopyranoses. J. Chem. Theory Comput. 2009, 5, 2353–2370. [Google Scholar] [CrossRef] [PubMed]
  16. Guvench, O.; Mallajosyula, S.S.; Raman, E.P.; Hatcher, E.; Vanommeslaeghe, K.; Foster, T.J.; Jamison, F.W.; MacKerell, A.D. CHARMM Additive All-Atom Force Field for Carbohydrate Derivatives and Its Utility in Polysaccharide and Carbohydrate–Protein Modeling. J. Chem. Theory Comput. 2011, 7, 3162–3180. [Google Scholar] [CrossRef] [PubMed]
  17. Morrow, B.H.; Payne, G.F.; Shen, J. pH-Responsive Self-Assembly of Polysaccharide through a Rugged Energy Landscape. J. Am. Chem. Soc. 2015, 137, 13024–13030. [Google Scholar] [CrossRef] [PubMed]
  18. Tsai, C.C.; Morrow, B.H.; Chen, W.; Payne, G.F.; Shen, J. Toward Understanding the Environmental Control of Hydrogel Film Properties: How Salt Modulates the Flexibility of Chitosan Chains. Macromolecules 2017, 50, 5946–5952. [Google Scholar] [CrossRef]
  19. Vanommeslaeghe, K.; Raman, E.P.; MacKerell, A.D. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. Journal of Chemical Information and Modeling 2012, 52, 3155–3168. [Google Scholar] [CrossRef]
  20. Guvench, O.; Greene, S.N.; Kamath, G.; Brady, J.W.; Venable, R.M.; Pastor, R.W.; Mackerell, A.D. Additive Empirical Force Field for Hexopyranose Monosaccharides. J. Comput. Chem. 2008, 29, 2543–2564. [Google Scholar] [CrossRef]
  21. Schatz, C.; Viton, C.; Delair, T.; Pichot, C.; Domard, A. Typical Physicochemical Behaviors of Chitosan in Aqueous Solution. Biomacromolecules 2003, 4, 641–648. [Google Scholar] [CrossRef] [PubMed]
  22. Romany, A.; Payne, G.F.; Shen, J. Mechanism of the Temperature-Dependent Self-Assembly and Polymorphism of Chitin. Chem. Mater. 2023, 35, 6472–6481. [Google Scholar] [CrossRef]
  23. Lamarque, G.; Lucas, J.M.; Viton, C.; Domard, A. Physicochemical Behavior of Homogeneous Series of Acetylated Chitosans in Aqueous Solution: Role of Various Structural Parameters. Biomacromolecules 2005, 6, 131–142. [Google Scholar] [CrossRef] [PubMed]
  24. Wallace, J.A.; Shen, J.K. Probing the Strand Orientation and Registry Alignment in the Propagation of Amyloid Fibrils. Biochemistry 2010, 49, 5290–5298. [Google Scholar] [CrossRef] [PubMed]
  25. Naito, P.K.; Ogawa, Y.; Sawada, D.; Nishiyama, Y.; Iwata, T.; Wada, M. X-Ray Crystal Structure of Anhydrous Chitosan at Atomic Resolution. Biopolymers 2016, 105, 361–368. [Google Scholar] [CrossRef] [PubMed]
  26. Okuyama, K.; Noguchi, K.; Miyazawa, T.; Yui, T.; Ogawa, K. Molecular and Crystal Structure of Hydrated Chitosan. Macromolecules 1997, 30, 5849–5855. [Google Scholar] [CrossRef]
  27. Lertworasirikul, A.; Yokoyama, S.; Noguchi, K.; Ogawa, K.; Okuyama, K. Molecular and Crystal Structures of Chitosan/HI Type I Salt Determined by X-ray Fiber Diffraction. Carbohydr. Res. 2004, 339, 825–833. [Google Scholar] [CrossRef] [PubMed]
  28. Yui, T.; Imada, K.; Okuyama, K.; Obata, Y.; Suzuki, K.; Ogawa, K. Molecular and Crystal Structure of the Anhydrous Form of Chitosan. Macromolecules 1994, 27, 7601–7605. [Google Scholar] [CrossRef]
  29. Mazeau, K.; Winter, W.T.; Chanzy, H. Molecular and Crystal Structure of a High-Temperature Polymorph of Chitosan from Electron Diffraction Data. Macromolecules 1994, 27, 7606–7612. [Google Scholar] [CrossRef]
  30. Okuyama, K.; Noguchi, K.; Hanafusa, Y.; Osawa, K.; Ogawa, K. Structural Study of Anhydrous Tendon Chitosan Obtained via Chitosan/Acetic Acid Complex. Int. J. Biol. Macromol. 1999, 26, 285–293. [Google Scholar] [CrossRef]
  31. Jang, M.K.; Kong, B.G.; Jeong, Y.I.; Lee, C.H.; Nah, J.W. Physicochemical Characterization of $\alpha$-Chitin, $\beta$-Chitin, and $\gamma$-Chitin Separated from Natural Resources. J. Polym. Sci. A Polym. Chem. 2004, 42, 3423–3432. [Google Scholar] [CrossRef]
  32. Fan, Y.; Saito, T.; Isogai, A. Preparation of Chitin Nanofibers from Squid Pen β-Chitin by Simple Mechanical Treatment under Acid Conditions. Biomacromolecules 2008, 9, 1919–1923. [Google Scholar] [CrossRef] [PubMed]
  33. Kaya, M.; Mujtaba, M.; Ehrlich, H.; Salaberria, A.M.; Baran, T.; Amemiya, C.T.; Galli, R.; Akyuz, L.; Sargin, I.; Labidi, J. On Chemistry of γ-Chitin. Carbohydr. Polym. 2017, 176, 177–186. [Google Scholar] [CrossRef] [PubMed]
  34. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Montgomery, Jr., J.A.; Vreven, T.; Kudin, K.N.; Burant, J.C.; et al. Gaussian 03, Revision C.02. Gaussian, Inc., Wallingford, CT, 2004.
  35. Kumar, A.; Yoluk, O.; MacKerell, A.D., Jr. FFParam: Standalone package for CHARMM additive and Drude polarizable force field parametrization of small molecules. Journal of Computational Chemistry 2020, 41, 958–970, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.26138. [Google Scholar] [CrossRef] [PubMed]
  36. Case, D.A.; Ben-Shalom, I.Y.; Brozell, S.R.; Cerutti, D.S.; Cheatham, III, T.; Cruzeiro, V.W.D.; Darden, T.A.; Duke, R.E.; Ghoreishi, D.; Gilson, M.K.; et al. AMBER 2020, 2020.
  37. Hanwell, M.D.; Curtis, D.E.; Lonie, D.C.; Vandermeersch, T.; Zurek, E.; Hutchison, G.R. Avogadro: An Advanced Semantic Chemical Editor, Visualization, and Analysis Platform. J. Cheminform. 2012, 4, 17. [Google Scholar] [CrossRef] [PubMed]
  38. Martinez, L.; Andrade, R.; Birgin, E.; Martinez, J. Packmol: A package for building initial configurations for molecular dynamics simulations. JOURNAL OF COMPUTATIONAL CHEMISTRY 2009, 30, 2157–2164. [Google Scholar]
  39. Brooks, B.; Brooks, C.; MacKerell, A.; Nilsson, L.; Petrella, R.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The Biomolecular Simulation Program. Journal of computational chemistry 2009, 30, 1545–1614. [Google Scholar] [CrossRef]
  40. Grønbech-Jensen, N.; Farago, O. A simple and effective Verlet-type algorithm for simulating Langevin dynamics. Molecular Physics 2013, 111, 983–991. [Google Scholar] [CrossRef]
  41. Chow, K.H.; Ferguson, D.M. Isothermal-isobaric molecular dynamics simulations with Monte Carlo volume sampling. Computer Physics Communicationo 1995, 91, 283–289. [Google Scholar] [CrossRef]
  42. Miyamoto, S.; Kollman, P.A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. Journal of Computational Chemistry 1992, 13, 952–962. [Google Scholar] [CrossRef]
  43. Roe, D.R.; Cheatham, T.E. PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. Journal of Chemical Theory and Computation 2013, 9, 3084–3095. [Google Scholar] [CrossRef] [PubMed]
  44. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. Journal of Molecular Graphics 1996, 14, 33–38. [Google Scholar] [CrossRef] [PubMed]
Preprints 93506 g001
Figure 1. Schematics of a simplified “chitosan matrix” [5]. (a) The degree of polymerization (DP), degree of acetylation (DA) and pattern of acetylation (PA) describe the sequence of chitosan [5,8]. Blue and red circles represent the glucosamine (GlcN) and acetylglucosamine (GlcNAc) monomers, respectively. The effects of DA and PA are studied in this work. (b) and (c) Two GlcN or GlcNAc monomers connected by the glycosidic bond. Atom names are given.
Figure 1. Schematics of a simplified “chitosan matrix” [5]. (a) The degree of polymerization (DP), degree of acetylation (DA) and pattern of acetylation (PA) describe the sequence of chitosan [5,8]. Blue and red circles represent the glucosamine (GlcN) and acetylglucosamine (GlcNAc) monomers, respectively. The effects of DA and PA are studied in this work. (b) and (c) Two GlcN or GlcNAc monomers connected by the glycosidic bond. Atom names are given.
Preprints 93506 g001
Preprints 93506 g002
Figure 2. Optimization of the dihedral parameters for β -glucosamine in chitosan. (a,b) Definitions of C3-C2-N-HN1(HN2) and H2-C2-N-HN1(HN2) torsion angles in β -glucosamine. (c,d) Comparison of the potential energy scans for the two torsion angles calculated by quantum mechanics (QM) method (red) and CHARMM36 force field (before optimization in black and after optimization in blue). For the QM calculations, the MP2/6-31g method was used following the protocol for CHARMM36 development [15,20].
Figure 2. Optimization of the dihedral parameters for β -glucosamine in chitosan. (a,b) Definitions of C3-C2-N-HN1(HN2) and H2-C2-N-HN1(HN2) torsion angles in β -glucosamine. (c,d) Comparison of the potential energy scans for the two torsion angles calculated by quantum mechanics (QM) method (red) and CHARMM36 force field (before optimization in black and after optimization in blue). For the QM calculations, the MP2/6-31g method was used following the protocol for CHARMM36 development [15,20].
Preprints 93506 g002
Preprints 93506 g003
Figure 3. Chitosan solubility and self-assembly is modulated by the degree of acetylation. (a) The representative snapshots of 50% DA chitosan after 3 μ s and 20% DA chitosan after 2 μ s simulations. Probability distribution of the % solvent exposure for chitin (green), chitosan at 50% DA (purple) and 20% DA (yellow) with the block PA. Final 100 ns were used for the calculations. Data for chitin ( D A of 100%) was taken from our previous work [22]. (b) Bar plot comparing the occupancies of the h-bonds between water and the acetylated nitrogen (N*) from acetylglucosamine or amine nitrogen (N) from glucosamine at 50% DA. Note, the nitrogen is acting as a h-bond acceptor. Final 100 ns of was used for calculation and results are an aggregate of three independent runs at 50% DA. All chains were used. (c) Probability distribution of the end-to-end distance of the 50% (purple) and 20% (brown) DA chitosan chains with the block PA in solution. Only the chains in solution were considered in the calculation. The end-to-end distance is defined as the distance between the center of mass of the second sugar unit and center of mass of the ninth sugar unit in a chitosan chain. The green dashed line represents the contour length (length in the built chitosan model).
Figure 3. Chitosan solubility and self-assembly is modulated by the degree of acetylation. (a) The representative snapshots of 50% DA chitosan after 3 μ s and 20% DA chitosan after 2 μ s simulations. Probability distribution of the % solvent exposure for chitin (green), chitosan at 50% DA (purple) and 20% DA (yellow) with the block PA. Final 100 ns were used for the calculations. Data for chitin ( D A of 100%) was taken from our previous work [22]. (b) Bar plot comparing the occupancies of the h-bonds between water and the acetylated nitrogen (N*) from acetylglucosamine or amine nitrogen (N) from glucosamine at 50% DA. Note, the nitrogen is acting as a h-bond acceptor. Final 100 ns of was used for calculation and results are an aggregate of three independent runs at 50% DA. All chains were used. (c) Probability distribution of the end-to-end distance of the 50% (purple) and 20% (brown) DA chitosan chains with the block PA in solution. Only the chains in solution were considered in the calculation. The end-to-end distance is defined as the distance between the center of mass of the second sugar unit and center of mass of the ninth sugar unit in a chitosan chain. The green dashed line represents the contour length (length in the built chitosan model).
Preprints 93506 g003
Preprints 93506 g004
Figure 4. Dynamics of the self-assembly process of chitosan with the block and alternating PAs. (a,b) Representative snapshots of chitosan with 50% DA and block (a) or alternating (b) PA along the local minima of the free energy surface (self-assembly pathway) shown in (c) and (d). Simulation run 1 of the block PA and run 3 of the alternating PA were used. The associated P 2 order parameters and average percent solvent exposure (relative to the fully solvated state) are given. (c,d) Free energy as a function of the radius of gyration and interchain h-bonds per monomer unit for chitosan with 50% DA and block (c) or alternating (d) PA. Only chains that formed part of the nanofibril were used for calculation. The runs with the highest number of interchain hydrogen bonding for both the block (run 1) and alternating (run 3) PAs are shown here, while the rest are shown in Supplemental Figure S4. The P 2 order parameter was calculated using Eq. 1, following our previous work [22,24]. (e) Snapshot of a nanofibril with water shown as red spheres. A zoomed-in view shows water-mediated interactions between two fibril sheets.
Figure 4. Dynamics of the self-assembly process of chitosan with the block and alternating PAs. (a,b) Representative snapshots of chitosan with 50% DA and block (a) or alternating (b) PA along the local minima of the free energy surface (self-assembly pathway) shown in (c) and (d). Simulation run 1 of the block PA and run 3 of the alternating PA were used. The associated P 2 order parameters and average percent solvent exposure (relative to the fully solvated state) are given. (c,d) Free energy as a function of the radius of gyration and interchain h-bonds per monomer unit for chitosan with 50% DA and block (c) or alternating (d) PA. Only chains that formed part of the nanofibril were used for calculation. The runs with the highest number of interchain hydrogen bonding for both the block (run 1) and alternating (run 3) PAs are shown here, while the rest are shown in Supplemental Figure S4. The P 2 order parameter was calculated using Eq. 1, following our previous work [22,24]. (e) Snapshot of a nanofibril with water shown as red spheres. A zoomed-in view shows water-mediated interactions between two fibril sheets.
Preprints 93506 g004
Preprints 93506 g005
Figure 5. Chitosan chains are predominantly antiparallel in the self-assembled nanofibrils with either block or alternating PA. (a,b) Probability distribution of the angles between the adjacent chains for the chitosan with 50% DA and block (a) or alternating (b) PA. An angle around 0° indicates a parallel-chain alignment while an angle around 180° indicates an antiparallel-chain alignment. The final 100 ns of simulation time of each of the three trajectories was used for the calculation. (c) Representative snapshots of two neighboring chains in the nanofibrils with the block and alternating PAs. The respective nanofibrils are displayed in Figure 4a and Figure 4b.
Figure 5. Chitosan chains are predominantly antiparallel in the self-assembled nanofibrils with either block or alternating PA. (a,b) Probability distribution of the angles between the adjacent chains for the chitosan with 50% DA and block (a) or alternating (b) PA. An angle around 0° indicates a parallel-chain alignment while an angle around 180° indicates an antiparallel-chain alignment. The final 100 ns of simulation time of each of the three trajectories was used for the calculation. (c) Representative snapshots of two neighboring chains in the nanofibrils with the block and alternating PAs. The respective nanofibrils are displayed in Figure 4a and Figure 4b.
Preprints 93506 g005
Preprints 93506 g006
Figure 6. Intermolecular hydrogen bonding pattern is similar for the nanofibrils with the block and alternating PAs. (a) Percentage occupancies of the intermolecular h-bonds in the self-assembled chitosan nanofibrils with 50% DA and block (maroon) or alternating PA (brown). (b) Percentage of the NH–O6 or O6H–O7 intermolecular h-bonds formed between two acetylated (NAc–NAc) units, between one acetylated and one deacetylated (NAc–N) units, and between two deacetylated units (N–N). The final 100 ns of simulation time in each of the three trajectories for each PA was used for the calculation. Only chains aligned antiparallel to each other were considered. (c) Two dimensional structure showing the intermolecular h-bonds between two acetylated units (top) and between one acetylated and one deacetylated units (bottom).
Figure 6. Intermolecular hydrogen bonding pattern is similar for the nanofibrils with the block and alternating PAs. (a) Percentage occupancies of the intermolecular h-bonds in the self-assembled chitosan nanofibrils with 50% DA and block (maroon) or alternating PA (brown). (b) Percentage of the NH–O6 or O6H–O7 intermolecular h-bonds formed between two acetylated (NAc–NAc) units, between one acetylated and one deacetylated (NAc–N) units, and between two deacetylated units (N–N). The final 100 ns of simulation time in each of the three trajectories for each PA was used for the calculation. Only chains aligned antiparallel to each other were considered. (c) Two dimensional structure showing the intermolecular h-bonds between two acetylated units (top) and between one acetylated and one deacetylated units (bottom).
Preprints 93506 g006
Table 1. Summary of de novo self-assembly MD simulations of chitosan chains in water
Table 1. Summary of de novo self-assembly MD simulations of chitosan chains in water
DA PA DP Chains Simulation time
20% Block 10 24 3 x 2 μ s
50% Block 10 24 3 x 3 μ s
50% Alternating 10 24 2 x 4 μ s, 4.5 μ s
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated