Effect of Structure and Hydration Level on Water Diffusion in Chitosan Membranes

DOI: 10.1002/mats.202000064 traditional semiconductor devices rely exclusively on electronic currents rather than molecular or ionic movement, they are unsuitable for communicating with living organisms. To address this challenge, bioprotonic devices interface with living organisms by incorporating components that facilitate both molecular/ ionic as well as electronic transport.[8] The bioprotonic devices reported to date, have focused on delivering one type of ion, H+, to biological systems.[9,10] Expanding these devices range of molecular/ionic species is an anticipated future advancement. A standard bioprotonic device consists of a conductive polymer membrane deposited onto an inert substrate (typically SiO2) that lies between two palladium hydride (PdHx) protodes. The device is turned on when a voltage is applied to its gate electrode, and the flow of protonic current is initiated.[11] Protons are generated at the source protode through hydrogen oxidation (PdH → Pd + H+ + e−)[8,12] and flow across the conductive polymer membrane toward the drain protode. Because the flow of protonic current in these devices depends on the transport properties of the polymer membrane, designing membranes with suitable molecular/ionic transport properties and biocompatibility is essential for the advancement of bioprotonics. One candidate polymer for bioprotonic devices is chitosan, a de-acetylated derivative of chitin.[5,6] Its chemical structure is shown in Figure 1. Chitosan is an appealing candidate because its biologically derived parent material, chitin, is widely available and can be extracted from the shells of crustaceans.[13] The design refinement of chitosan-based membranes for bioprotonic applications necessitates an understanding of the factors that influence the polymer’s transport properties. While there have been several experimental and theoretical studies of ionic diffusion in chitosan-based membranes,[14–17] the diffusion of molecular species in these membranes is less well-characterized. Understanding molecular diffusion is essential for the design of not only devices that deliver small molecules but also those that deliver ions, for example, in cases where ionic transport depends on a network of solvent molecules. In this article, we present results from ab intio molecular dynamics (AIMD) simulations showing how membrane structure and hydration level influence water diffusion in chitosan membranes. We consider a set of six model chitosan membranes at varying levels of hydration. Based on our AIMD simulations, we calculate the structure of each membrane and the diffusion coefficient of water within each hydrated membrane. The factors that influence water diffusion are explored, and Bioprotonic devices show promise for medicinal/therapeutic treatments because of their ability to deliver a controlled flow of small molecules and ions to living organisms. These devices rely on biocompatible, conductive polymer membranes to facilitate molecular/ionic transport. Herein, ab initio molecular dynamics simulations are used to probe the effects structure and hydration level on water diffusion in chitosan-based polymer membranes. The diffusion coefficient of water is shown to be highest in membranes with low overall densities and weak non-covalent interactions. Insight from these results may aid in optimizing the molecular/ionic transport properties of polymer membranes for bioprotonic applications.


Introduction
Many physiological processes depend on the movement of small molecules and ions. Two key examples are cells within the human body pumping sodium and potassium ions to maintain their volume, [1] and aquaporin proteins facilitating water transport between kidney cells to maintain urinary function. [2] Adverse health effects are known to occur in cases where the body's intrinsic molecular/ionic transport mechanisms fail. [3,4] Because signaling in physiological processes depends on the movement of small molecules and ions, devices that drive the movement of molecular and/or ionic species could in principle communicate with living organisms and mimic physiological functions. These devices fall under the emerging field of bioprotonics, which shows promise for advancing medicinal/ therapeutic treatments. [5][6][7] Establishing device-organism communication is the first step toward producing an effective bioprotonic system. As

Structural Model Generation and Potential Energy Curve Mapping
A structural model of the OH-terminated chitosan dimer shown in Figure 2a was constructed using the GaussView 5.0 GUI software program. [18] The molecular geometry of this starting structure was optimized at the PBE-D3/6-31G* level of theory [19,20] using Gaussian. [21] Two optimized chitosan dimers were then arranged in a stacked configuration and taken to be a model representation of the unhydrated chitosan membrane shown in Figure 2b. To determine the structure of this model system, a potential energy curve (PEC) was mapped as a function of the separation between the dimers. For the stacked system, the separation between each atom in one dimer and its equivalent atom in the other dimer, was systematically increased from 3.00 to 8.00 Å, in 0.25 Å steps. At each step, the system was optimized at the PBE-D3/6-31G* [19,20] level of theory with a total of three constraints: the distances between the oxygen atoms that link and terminate the pyranose rings and their equivalent atoms in the neighboring dimer, were held constant. Twenty-one structures were used to map the PEC.

Sampling the Isothermal-Isobaric Ensemble
AIMD simulations sampling the isothermal-isobaric (NPT) ensemble were performed using the Quickstep [22,23] module within CP2K [24] to probe the effects of water content on chitosan membrane structure. Six model systems were considered, each containing one chitosan dimer and a varying number of water molecules. The weight percent (wt%) of water in these models varies between 0 and 64. Periodic boundary conditions were enforced to simulate continuous polymer strands, rather than OH-terminated dimers. Table 1 summarizes the compositions of the six model systems. Because the level of membrane hydration is a design parameter that will depend on the device's envisioned application, the model systems were chosen to represent a range from low (18 water wt%) to intermediate (64 water wt%) hydration.
The initial length of the lattice constant in the polymer's stacking direction was chosen based on the preferred value of separation determined by PEC mapping, as described in Section 2.1. The initial values of the other lattice parameters were chosen to enforce periodicity in the polymer strands and to accommodate water molecules within the simulation cells. Water molecules were randomly placed in the vacant space within each simulation cell. They were not arranged in any chemically meaningful way, for example, to predispose the membranes to a chosen hydrogen-bonding configuration. For each model system, a 10-ps trajectory with a time-step of 0.5 fs, a target temperature of 300 K, and a target pressure of 1 atm, was calculated. The CSVR thermostat [25] with a time constant of 100 fs was used for these trajectories. To make Born-Oppenheimer molecular dynamics simulations of these large systems (56-140 atoms per unit cell) more tractable, the electronic structure was described using the Gaussian and plane waves method. [26] The revPBE [27] exchange-correlation functional together with Goedecker-Teter-Hutter (GTH)  Figure 1. Chemical structure of chitosan.

Figure 2.
Structural representations of a) the optimized chitosan dimer and b) the unoptimized dimer stack (not to scale). Carbon atoms are grey, nitrogens are blue, oxygens are red, and hydrogens are white. Green stars indicate oxygen atoms involved in the constraints during constrained optimizations of the dimer stack. The green arrow indicates the axis of separation used in potential energy curve mapping. GaussView [18] was used to visualize the chemical structures.
pseudopotentials, [28,29] the TZV2P basis set, and Grimme's D3 dispersion correction [20] were used for all calculations. Energy cutoffs of 600 and 60 Ry were used for the plane-waves and real-space grid for Gaussian mapping, respectively. The structure of the chitosan membrane at each level of hydration was determined from the results of these simulations.

Sampling the Canonical Ensemble
AIMD simulations sampling the canonical (NVT) ensemble were carried out to probe the effects of membrane structure and hydration level on water diffusion. Forty 20-ps trajectories with a time-step of 0.5 fs and a target temperature of 300 K were calculated for each hydrated membrane, systems 2-6 in Table 1. The CSVR thermostat [25] with a time constant of 1000 fs was used for these trajectories. The electronic structure was described at the revPBE-D3/GTH/TZV2P [20,[27][28][29] level of theory with plane-wave and real-space grid cutoffs of 300 and 30 Ry, respectively. The initial geometries for these trajectories were sampled from separate 5-ps trajectories and the lattice parameters were taken from the results of simulations sampling the NPT ensemble, as described in Section 2.2.1. The first picosecond of each trajectory was discarded for equilibration, and the subsequent 19 ps were used for analysis. In total, 4 ns of simulation time sampling the NVT ensemble were calculated.
To determine the diffusion coefficient of water within each hydrated membrane, the mean-squared displacement (MSD) of the water molecules was calculated from the positions of their oxygen atoms at each time-step. The expressions to calculate the MSD in each direction x, y, z from simulations of a periodic system are, , , where m is the trajectory index, i is the index oxygen atom index, COM refers to the center of mass of the simulation cell, t is time, N is the total number of water molecules, and M is the total number of trajectories. The water diffusion coefficient D within each hydrated membrane was calculated from 〈x 2 (t)〉 + 〈y 2 (t)〉 by fitting, using the expression, [30] ( ) ( ) 4 Further details about the expressions used to obtain the diffusion coefficients are described in the Supporting Information.

Mapped Potential Energy Curve
The computed potential energy points and a best-fit curve for separation of the stacked chitosan dimers are shown in Figure 3. This curve is based on constrained optimizations at the PBE-D3/6-31G* level of theory. The minimum corresponds to a dimer separation of ≈5.0 Å. Visual inspection of the lowest energy structure identified reveals that the dimers slide with respect to each other to facilitate intermolecular hydrogen-bonding. For reference, the energy associated with a moderate strength hydrogen bond is ≈0.17-0.43 eV. [31] Three inter-dimer hydrogen bonds: NH···O, NH···N, and OH···O are present in the lowest energy structure identified. This result is consistent with the X-ray diffraction study by Yui and colleagues, [32] that describes the role of hydrogenbonding in stabilizing the structure of anhydrous, crystalline chitosan. The predicted minimum energy structure for the dimer stack was used to inform choices for initial molecular geometries in a set of AIMD simulations sampling the NPT ensemble.
Macromol. Theory Simul. 2021, 30,2000064 Figure 3. PBE potential energy points and best-fit fourth degree polynomial for the separation of two OH-terminated chitosan dimers. The separation between dimers is defined by the distance between the oxygen atoms that link the pyranose rings. The energy scale is relative to the lowest energy structure identified.

Structure of Hydrated Chitosan Membranes
The structure of the chitosan membrane was determined at six different levels of hydration, 0-64 water wt%, from a set of trajectories sampling the NPT ensemble. Convergence with respect to the average unit cell volume was achieved for all systems considered. These results are reported in the Supporting Information. The density of each hydrated membrane was calculated from the masses of the polymer and water molecules, and the average cell volume ρ = m/V. The density of unhydrated chitosan is predicted to be 1.38 g cm −3 , consistent with the experimental study of Sakurai et al. [33] that reports chitosan film densities within 1.36-1.47 g cm −3 . Table 2 shows that the membrane density decreases with increasing water content, and approaches the bulk density of water.
Visual inspection of each hydrated system reveals the formation of water channels perpendicular to the planes of the pyranose rings within the chitosan strands. Figure 4 shows the membrane's molecular geometry at a hydration level of 18 water wt% (4 waters per dimer). The formation of water channels is consistent across all levels of membrane hydration considered. The water channel width as defined by the shortest nitrogen-nitrogen distance between adjacent polymer strands ranges between ≈8-25 Å.
Similar to the results of PEC mapping shown in Figure 3, these simulations predict that the chitosan strands maintain a close separation, ≈4-5 Å, to facilitate hydrogen-bonding. The waters do not insert themselves between the strands at any level of membrane hydration considered because doing so would require breaking at least one hydrogen bond. Recall that the energy associated with a moderate strength hydrogen bond is ≈0.17-0.43 eV, [31] while k B T at 300 K is 0.026 eV.
It is important to note that our calculations are carried out with perfectly crystalline polymers; these model systems represent ideal membranes. Although there are experimental reports [35,36] of highly ordered polymeric systems with crystalline and semi-crystalline structures, it is unlikely for a real polymer membrane to have an ideal structure like those represented by our models. Rather than representing examples of membranes that could be prepared experimentally, our systems serve as models to explore the factors that influence water diffusion. We expect the water channels in a real chitosan membrane to be more tortuous.

Water Diffusion Coefficients
The diffusion coefficient of water within each hydrated chitosan membrane was computed from the results of trajectories sampling the NVT ensemble. For each system, the MSD of the water molecules was calculated from the positions of their oxygen atoms using Equation (1). We define an orthogonal coordinate system such that z is parallel to the planes of the polymer's pyranose rings and perpendicular to the lengthwise direction of the polymer strands, y is perpendicular to the planes of the pyranose rings and the lengthwise direction of the strands, and x is parallel to the planes of the pyranose rings and the lengthwise direction of the strands. The water molecules are free to move in the xy-plane but are inhibited by chitosan in the z-direction. Figure 5 shows the MSD of water in the xy-plane and in the z-direction, as a function of time. From the results shown in Figure 5b, it is clear that water mobility is limited in the z-direction because the chitosan strands serve as walls that inhibit the movement of diffusing water molecules. We observe normal diffusion in the xy-plane, that is, 〈x 2 (t)〉 + 〈y 2 (t)〉 ∝ t.
The diffusion coefficient of water within each hydrated membrane was calculated from the results of Figure 5a according to Equation (2), by fitting the region of each curve between 4 and 19 ps to a line. Bootstrapping analysis [37] with 10 000 samples per hydration level was performed to quantify the errors in the diffusion coefficients. The calculated water diffusion coefficients and associated bootstrap standard errors are summarized in Table 2. Based on the results of t-tests, [38] the differences in Macromol. Theory Simul. 2021, 30,2000064 Figure 4. Snapshot of water channels formed within a 3 × 3 × 3 supercell of the chitosan membrane at 18 water wt% (4 waters per dimer). Carbons are grey, nitrogens are blue, oxygens are red, and hydrogens are white. The chemical structure was visualized using Jmol. [34] the computed diffusion coefficients are not statistically significant at a 95% confidence level for any pair of membrane hydration levels considered. A comparison of diffusion coefficients at different hydration levels is therefore not meaningful for this set of systems. This is unsurprising given the magnitudes bootstrap standard errors associated with the computed diffusion coefficients. The results of t-tests are reported in the Supporting Information. An additional set of MD trajectories sampling the NVT ensemble were calculated at the revPBE-D3/GTH/ TZV2P [20,[27][28][29] level of theory for a cubic simulation cell containing 64 water molecules. The water diffusion coefficient computed from these trajectories is 0.164 ± 0.00831 Å 2 ps −1 . For comparison, the experimentally reported diffusion coefficient in bulk water at 298 K is 0.2299 Å 2 ps. −1 [39] Nearly all computed membrane water diffusion coefficients reported in Table 2 are lower than the theoretical value for bulk water, and all are lower than the experimental value.

Sensitivity to Water Channel Orientation
In a real polymer membrane, the system's structure and consequently its transport properties are extremely sensitive to the method of preparation. [40] In cases where nanoscale control of the membrane's structure is challenging to achieve experimentally, the effect of water channel orientation with respect to the polymer strands on water diffusion may be relevant to the design of membranes for bioprotonic applications. To explore the sensitivity of water diffusion to water channel orientation, we considered a new set of model chitosan membranes with different geometries than those previously described. Their chemical compositions are the same as those listed in Table 1.

Structure and Water Channel Orientation in Chitosan
AIMD simulations sampling the NPT ensemble were carried out for all systems according to the procedure described in Section 2.2.1. A potential energy curve was mapped as a function of the separation between two in-plane chitosan dimers. Similar to the results shown in Figure 3 the dimers are predicted to maintain a separation (as defined by the distance between pyranose-linking oxygens) of approximately 8 Å to facilitate hydrogen-bonding. The lowest energy dimer separation was used to inform choices for the initial lattice parameters of the periodic polymer systems. They were chosen such that the chitosan strands have a close in-plane separation and are more widely spaced in the stacking direction to accommodate water molecules within the simulation cells. AIMD simulations of membranes with this new starting configuration predict the formation of water channels parallel to the planes of chitosan's pyranose rings as shown Figure 6. The water channel width as defined by the shortest nitrogen-nitrogen distance between adjacent polymer strands ranges between ≈8-18 Å. In addition to a different water channel orientation, the membrane densities reported in Table 3 are lower than those for systems with water channels perpendicular to the planes of the pyranose rings.

Diffusion and Water Channel Orientation in Chitosan
AIMD simulations sampling the NVT ensemble were carried out according to the procedures described in Section 2.2.2. The diffusion coefficient of water within each hydrated membrane was calculated using Equations (1) and (2). We define an orthogonal coordinate system such that z is perpendicular to the planes of the polymer's pyranose rings and the lengthwise direction of the polymer strands, y is parallel to the planes of the pyranose rings and perpendicular to the lengthwise direction of the strands, and x is parallel to the lengthwise direction of the strands. The waters move freely in the xy-plane but are inhibited by chitosan in the z-direction. Figure 7 shows the mean-squared displacements of the water molecules as a function of time. Compared to the membranes with water channels perpendicular to the polymer strands, these membranes exhibit larger MSDs. We attribute this to the lower membrane densities in the parallel orientation. 〈z 2 (t)〉 is again small because the chitosan strands inhibit water diffusion in this direction. We observe normal diffusion in the xyplane, that is, 〈x 2 (t)〉 + 〈y 2 (t)〉 ∝ t with some slight deviations at 18 and 57 water wt%.
The water diffusion coefficient in each system was calculated according to Equation (2), by fitting the region of each curve between 4 and 19 ps to a line. Bootstrapping analysis [37] with 10 000 samples per hydration level was performed to quantify the errors in the diffusion coefficients. The calculated diffusion coefficients and associated bootstrap standard errors are summarized in Table 3. Based on the results of t-tests, [38] the differences in the computed water diffusion coefficients are statistically significant at a 95% confidence level for the membrane at hydration levels of: 1) 31 and 57 water wt%, 2) 31 and 64 water wt%, 3) 47 and 57 water wt%, and 4) 47 and 64 water wt%. The results of t-tests are reported in the Supporting Information.
The water diffusion coefficients are overall greater for this set of membranes because the lower membrane densities result in larger water MSDs. They are also somewhat inversely correlated with the overall densities. It is clear from the results in Table 3 however, that density is not the only factor contributing to the predicted relationship between the membrane's hydration level and its water diffusion coefficient.

Non-Covalent Interactions
Given the chemical nature of these systems, we hypothesize that attractive interactions between non-bonded groups that can form hydrogen bonds are one factor that can contribute to a decrease in mobility with increasing water content. To test this hypothesis, the average water-water and amine-water distances within each membrane were calculated from the computed NVT trajectories. These distances are based on the positions of the oxygen atoms within the waters and the nitrogen atoms within chitosan. Figure 8a shows the average water-water distances in the hydrated chitosan membranes. The first nearest distances fall within the expected range for moderate strength hydrogenbonding (2.5-3.2 Å [41] ) at hydration levels of 47-64 water wt% and weak hydrogen bonding (3.2-4.0 Å [41] ) at 18-31 water wt%. The second nearest distances fall within the range for moderate hydrogen-bonding at hydration levels of 47 and 64 water wt%, weak hydrogen-bonding at 57 water wt%, and electrostatic interactions at 18-31 water wt%. The third nearest distances fall within the range for weak hydrogen-bonding at hydration levels of 47-64 water wt%, and electrostatic interactions at Macromol. Theory Simul. 2021, 30, 2000064   Figure 6. Snapshot of water channels formed within a 3 × 3 × 3 supercell of the chitosan membrane at 18 water wt% (4 waters per dimer). Water channels form parallel to the planes of the pyranose rings. Carbons are grey, nitrogens are blue, oxygens are red, and hydrogens are white. The chemical structure was visualized using Jmol. [34]  18-31 water wt%. At higher levels of hydration, the water molecules are overall closer and the distances corresponding to first, second, third nearest neighbors for each individual system are similar. The close proximity of water molecules at higher levels of hydration suggests that the non-covalent interactions between water molecules are overall stronger.
Consider the membrane at a hydration level of 47 water wt%. The low water diffusion coefficient in this system as reported in Table 3 can be explained in part, by the small distances between its water molecules as shown in Figure 8a. Smaller water-water distances are with consistent stronger hydrogenbonding/electrostatic interactions which in turn can result in an decrease in water mobility and lower water diffusion coefficient. Palin et al. [42] describe a similar finding based on experimental results for a set of chemically similar systems, hydrated cellulose acetate membranes. They report that an increase in water clustering with increasing water concentration leads to a decrease in the measured water diffusion coefficient in these systems. [42] The similarities in the first, second, and third nearest neighbor distances in our chitosan membrane at 47 water wt% could indicate similar water aggregation. Figure 8b shows the average amine-water distances in the hydrated chitosan membranes. These distances fall within the expected ranges for moderate hydrogen-bonding to electrostatic interactions. [41] As in the case of the water-water distances shown in Figure 8a, the amine-water distances are also small at a hydration level of 47 water wt%. This further supports the conclusion that small separations between certain non-bonded units within the membrane can lead to lower water mobility as a consequence of stronger non-covalent interactions. We expect that other hydrated membranes based on chemically similar polymers, that is, those that contain functional groups that can participate in hydrogen bonding, would exhibit similar behavior.   A second point of consideration is the membrane water diffusion coefficients at 47 and 57 water wt%. Although the waterwater and amine-water distances in these two systems differ by <2%, the diffusion coefficient at 57 water wt% is over two times greater than at 47 water wt%. This increase in water diffusion coefficient can be explained in part by the decrease in membrane density from 1.11 g cm −3 at 47 water wt% to 0.95 g cm −3 at 57 water wt%. One change in membrane structure that accompanies this decrease in membrane density is a widening of the water channels by ∼10 Å. With more free space for diffusing water molecules within the channels, the water diffusion coefficient increases at lower membrane density. Moreover, there may be fewer interactions between chemical groups other than those considered in Figure 8 in cases where adjacent chitosan strands are more widely spaced. Overall, a lower membrane density is one quantity that can reflect a change in membrane structure and potentially indicate a higher water diffusion coefficient.
From these results, it is clear that the membrane's density and its non-covalent interactions are two experimentally detectable features that influence water diffusion. These findings suggest that a membrane with a low overall density and weakly interacting non-bonded groups would be most effective in facilitating diffusion. When designing a membrane for a bioprotonic device, these two features can be used to assess the membrane's suitability for facilitating the diffusion of molecules. It is necessary however, to consider both features together when making predictions about the transport properties of a candidate membrane.
To illustrate, consider the chitosan membrane at hydration levels of 18 and 57 water wt% (4 and 24 waters per dimer, respectively). The predicted water diffusion coefficients for these two systems are similar, and their difference is not statistically significant. At a hydration level of 18 water wt%, the membrane's density is 1.08 g cm −3 and the water-water distances are large. Conversely, at 57 water wt%, the membrane's density is 0.95 g cm −3 , and the water-water distances are small. This finding shows that a high membrane density can be compensated for by weak non-covalent interactions and vice versa. It is therefore necessary to consider both polymer membrane features together, rather than either feature individually.
In addition, water diffusion in polymer membranes may also influence ionic diffusion in cases where the mechanisms of ionic transport depend on a network of water molecules. Proton movement for example, involves both structural (H + hopping) and vehicular (center of mass H 3 O + movement) diffusion. [43] A membrane with a low water diffusion coefficient may also exhibit low H 3 O + mobility, whereas a membrane with a widely spaced network of water molecules may be ineffective at facilitating H + hopping.

Conclusions
We determined the effects of membrane structure and hydration level on water diffusion in chitosan membranes using AIMD simulations. Our results show that the un-hydrated membrane is stabilized by hydrogen-bonding between chitosan strands. Upon hydration, the polymer maintains its hydrogen bonds and water channels form. Rather than a direct relationship between membrane hydration level and water diffusion, we find that two features of the chitosan membrane: 1) density and 2) non-covalent interactions, influence water diffusion. Specifically, the diffusion coefficient of water is higher in membranes with low densities, and widely spaced water-water and amine-water pairs.
Our results provide a basis to make qualitative inferences about the transport properties of candidate polymer membranes from these experimentally detectable features. It is important to bear in mind however, that these two membrane features must be considered together. In terms of influencing diffusion, a high density can be compensated for by weak noncovalent interactions and vice versa. We believe that by demonstrating how these two features of the chitosan membrane influence water diffusion, this work provides insight that could help guide the design of polymer membranes with transport properties that are suitable for bioprotonic applications.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.