A multiscale coarse-grained model to predict the molecular architecture and drug transport properties of modified chitosan hydrogels

Hydrogels constructed with functionalized polysaccharides are of interest in a multitude of applications, chiefly the design of therapeutic and regenerative formulations. Tailoring the chemical modification of polysaccharide-based hydrogels to achieve specific drug release properties involves the optimization of many tunable parameters, including (i) the type, degree (χ), and pattern of the functional groups, (ii) the water-polymer ratio, and (iii) the drug payload. To guide the design of modified polysaccharide hydrogels for drug release, we have developed a computational toolbox that predicts the structure and physicochemical properties of acylated chitosan chains, and their impact on the transport of drug molecules. Herein, we present a multiscale coarse-grained model to investigate the structure of networks of chitosan chains modified with acetyl, butanoyl, or heptanoyl moieties, as well as the diffusion of drugs doxorubicin (Dox) and gemcitabine (Gem) through the resulting networks. The model predicts the formation of different network structures, in particular the hydrophobically-driven transition from a uniform to a cluster/channel morphology and the formation of fibers of chitin chains. The model also describes the impact of structural and physicochemical properties on drug transport, which was confirmed experimentally by measuring Dox and Gem diffusion through an ensemble of modified chitosan hydrogels.


Introduction
Polysaccharides have received considerable interest as base materials in the development of hydrogels for biomedical applications (e.g., drug delivery and tissue engineering), [1][2][3] owing to their natural abundance, biocompatibility, and chemical versatility. 4 Among them, chitosan plays a prominent role as inexpensive, FDA-approved, biocompatible, and bioresorbable biopolymer. 5,6 Tailoring pharmaceutical hydrogels requires careful tuning of the physicochemical properties of the polymer chains, which determine the microstructure of the gel, and ultimately the kinetics of drug transport and release.
The network morphology of native chitosan hydrogels is determined by the ratio of acetylglucosamine (GlcNAc) and glucosamine (GlcN) monomers, expressed as degree of acetylation (w Ac ), 7,8 the water/chitosan ratio, and pH. 9 Tailored chemical modification of the free amine groups on GlcN monomers provides additional control over the properties of the chitosan hydrogels. 10,11 Chemical modification provides several design parameters, such as type of functional group, degree of modification (w), conjugation strategy, polymer/water ratio, etc. 1,[12][13][14] Understanding and predicting the impact of these design factors on the morphology and physicochemical properties of the chain network -and therefore its permeability to bioactive compounds -is key to develop drug delivery formulations with optimal therapeutic activity.
Computational modeling represents a powerful tool to gain a molecular-level understanding of the correlations between molecular interactions, the structure of the polysaccharide network, and drug diffusion therein. Molecular dynamics (MD) simulations with atomistic resolution provide a detailed picture of the interactions. The length-and time-scales of the dynamics in polysaccharide assemblies, however, far exceed the power of atomistic models. Coarse-grained (CG) models overcome these limitations by reducing the degrees of freedom of the system and accelerate molecular dynamics by creating a smoother energy landscape. 15,16 As a result, CG models can simulate larger systems and for longer time-scales. This gain in efficiency, on the other hand, comes at the cost of losing some precision in the molecular detail. To maintain predictive power in CG models, different strategies have been developed to find effective interaction potentials between CG sites, which preserve the molecular information of the material and provide a realistic representation of its large-scale behavior. CG strategies applied to polysaccharide systems include (i) parametrization to reproduce bulk thermodynamic data, 17 (ii) sampling the polymer conformations based on the conformational space available to the glycosidic angles j and c, [18][19][20] and (iii) interaction potentials that reproduce features of the atomistic system. [21][22][23][24][25][26] To date, a limited number of studies have focused on the MD simulation of chitosan-based materials. The conformations of single chitosan chains have been explored both at atomistic 19,[27][28][29] and CG resolution. 20 Other studies have probed the properties of nanoscale aggregates at the all-atom scale. [30][31][32] Finally, the aggregation of chitosan with different w Ac and different modification patterns has been simulated with a MARTINI-like model. 7,33 In the present study, we implement a novel and promising approach for elucidating the structure-property relationship in chemically modified chitosan hydrogels. Our method employs (1) multi-scale coarse-graining (MS-CG), also known as ''force matching'', 34 to derive non-bonded interactions for solute-solute, solute-solvent, and solvent-solvent interactions that reproduce the forces sampled in the all-atom system; and (2) Boltzmann Inversion 35 to derive the bonded interaction potentials. This combination enables polysaccharide models that are transferable to simulate polymer chains with different length and concentrations, 25 and reproduces the aggregation behavior and osmotic pressure of the atomistic systems. 26 In this study, we applied this method to model hydrated networks of chitosan chains modified with acetyl, butanoyl, and heptanoyl moieties. The CG model was first validated against atomistic MD simulations and subsequently applied to describe the effect of type, degree (w), and pattern of modification (i.e., alternating or blocky), as well as water/chitosan ratio on the molecular architecture of the network. Our results demonstrate that hydrophobic modifications significantly alter the conformation and spatial arrangement of the chains, causing a transition from a homogeneous distribution to a cluster-channel morphology featuring large chain aggregates and pores. The impact of network structure and properties on drug transport was investigated using two model drug molecules, gemcitabine (Gem) and doxorubicin (Dox), which differ by polarity, hydrophobicity, and size. The calculated diffusion coefficient of Dox was found to vary significantly with the type, degree (w), and pattern of modification, whereas Gem diffusion was approximately constant in all systems. Both predictions were confirmed experimentally by drug release tests conducted with an ensemble of modified chitosan gels. Collectively, these results demonstrate the validity of the proposed CG model and its value in pharmaceutical chemistry to accelerate the design, and therefore the transition to clinics, of chitosan hydrogels for drug delivery.

Atomistic simulations
Systems containing 10 chitosan chains with degree of polymerization (DP, defined as the sum of the number of N-acetylglucosamine (GlcNAc) and glucosamine (GlcN) monomers) of 16 and surrounded by 2000 water molecules (i.e., 200 water molecules per chain) were simulated at atomistic resolution in a cubic simulation box with side lengths between 4.44 and 4.62 nm. The following systems were modeled: (i) full chitosan (w = 0%); (ii) acetylated chitosan with w Ac = 16%, 24%, 32%, and 50%; (iii) butanoyl-chitosan with w But = 16%, 24%, and 32%; and (iv) heptanoyl-chitosan with w Hep = 8%, 16%, 20%, and 24%. The modification groups were distributed along the chitosan chains either uniformly spaced or in blocks of four. Single chains from systems (i)-(iv) were also simulated in a cubic 6.97 nm water box as controls. Solvated drug systems containing 50 drug molecules (Dox or Gem) and 10 water molecules per drug were modeled in a cubic box with side length of either 3.85 nm (Dox) or 3.56 nm (Gem). Finally, systems (i)-(iv) combined with either 10 Dox or 10 Gem molecules were simulated in a 4.55 nm cubic water box. All simulations used periodic boundary conditions. The initial structures of the chitosan molecules were constructed with tleap, 36 their topologies were converted to GROMACS format using the glycam2gmx script, 37,38 and solvated in GROMACS. 39 The coordinates for the drug molecules were obtained from the GROMOS ATB repository. 40 All simulations were performed in GROMACS 5.1.2. 41 The GLYCAM06 TIP5P OSMO,r14 force field 42,43 was used for chitosan, and the TIP5P model for the water molecules. 44 The parameters for Dox and Gem were taken from the general AMBER (GAFF) force field. 45 The partial charges for the modified GlcN monomers, Dox, and Gem were calculated following the GLYCAM06 protocol 42 using the R.E.D. scripts. 46 A cut-off of 1.4 nm was used for Lennard-Jones and electrostatic interactions. Long-range electrostatic interactions were evaluated using particle mesh Ewald. 47 Covalent bonds involving hydrogen atoms were constrained with the LINCS algorithm, 48 while water molecules were kept rigid using SETTLE. 49 Following energy minimization, a 50 ns NPT equilibration run was conducted at the temperature of 300 K and pressure of 1 bar using the Nosé-Hoover thermostat 50,51 and the Parrinello-Rahman barostat. 52 Subsequently, a 400 ns NVT equilibration run was performed using the average box size extracted from the NPT trajectory and the Nosé-Hoover thermostat, 50,51 followed by a 100 ns production MD run. A time step of 2 fs was used, and energy and pressure dispersion corrections were applied where appropriate. To extract forces for the CG procedure, separate reruns of the MD trajectories containing only solute-solute, solute-solvent, or solvent-solvent interaction were calculated 25 using the reaction-field method 53 for long-range electrostatics.

Coarse-grained model
Following previous CG models, 25 GlcN monomers were mapped to 3 CG sites (A, B and C), as shown in Fig. 1a-c. Acetyl, butanoyl or heptanoyl modifications were mapped into one (M), two (MA and MB), or three (MA, MB, and MC) CG sites, respectively, using center-of-mass mapping. Water molecules were mapped to one CG site located at their center of geometry. The potentials for bond, angle, dihedral interactions, as well as specific non-bonded 1-3 and 1-4 interactions (i.e., between third and fourth bonded neighbors) were obtained from Boltzmann inversion 35 using the VOTCA package. 23,87 The list of all bonded interactions is reported in Sections S1 and S2 (ESI †). The non-bonded interaction potentials were obtained from the MS-CG method to match the forces sampled at the atomic scale 34 using the re-run trajectories with separate solute-solute, solute-solvent, and solvent-solvent interactions (note: these interaction potentials are not by default Lennard-Jones-like and can, in principle, feature different profiles, while resembling the Lennard-Jones shape, as shown in Fig. 2). The bonded and intra-molecular interactions obtained from the Boltzmann inversion were excluded from the MS-CG procedure. Dox was mapped to 11 CG sites using 9 different CG bead types (DA through DI, Fig. 1d) and arranged in a planar geometry. Gem was mapped to 4 CG sites (GA through GD, Fig. 1e). The CG interactions for the drug molecules were obtained by the same procedure as described above for chitosan.

Coarse-grained simulations
The following CG simulations of systems (i)-(iv) were performed using GROMACS 4.6.4: (a) 10 chitosan chains with DP = 16, as simulated at all-atom resolution; (b) 50 chains with DP = 50 and evenly spaced modifications, solvated with 80 000 water beads (32 water molecules per monomer); (c) 50 chains with DP = 50 and modifications grouped in blocks of four, solvated with 80 000 water beads; (d) 20 chains with DP = 50 and evenly spaced modifications, solvated with 100 000 water beads (100 water molecules per monomer); (e) 20 chains with DP = 50 and modifications grouped in blocks of four, solvated with 100 000 water beads. The initial structures of group (a) were obtained from the previous atomistic representations; groups (b)-(d) were initially constructed with atomistic resolution, and the corresponding initial structures and box-sizes for the CG simulations were obtained from 10 ns atomistic NPT simulations. All CG simulations were performed in the NVT ensemble using the leap-frog integrator with time steps of 1 fs and the Nosé-Hoover thermostat. 50,51 CG simulations of group (a) were run for 10 ns to obtain equilibrated radial distribution functions (RDFs) and angle distributions. Simulations of 100 ns were run for groups (b-d), and the data from the last 10 ns was used for analysis.
The network structures were characterized by calculating the distribution of pore size, following a published method; 54 this procedure finds the largest sphere that can be constructed to contain randomly selected points in the network, and was implemented in this work with a constrained nonlinear optimization of the center of the sphere using the SOLVOPT routine. 55 For every distribution, three snapshots of the network structures were taken at 5 ns intervals and analyzed. Finally, the diffusion coefficients of Dox and Gem in modified chitosan networks were obtained from the slope of their mean square displacement (MSD) as a function of time, and the corresponding van Hove functions were calculated using the built-in GROMACS routine. The numbers of contacts within 0.6 nm between the modification beads (M in acetyl-chitosan; MA and MB in butanoyl-chitosan; and MA, MB, and MC in heptanoylchitosan) in the network, as well as between the drug and the modification beads or between the drug and the backbone beads (A, B, and C) were calculated using the GROMACS mindist routine, counting only one contact per bead.

Experimental diffusion coefficient of Dox and Gem in modified chitosan hydrogels
Hydrogels comprising acetyl-chitosan (w Ac = 16% or 50%), butanoyl-chitosan (w Ac = 16% or 32%), or heptanoyl-chitosan (w Hep = 8% or 24%) were constructed and loaded with either Dox or Gem at the same concentrations utilized in the simulations, following a procedure that we previously developed. 56 Every drug-loaded gel was put into contact with an aqueous supernatant, which was sampled at set time points (1,4,8,24,48, and 72 hours) and analyzed by UV spectrophotometry to determine the amount of drug released from the gels into solution. The values of diffusion coefficient were finally calculated from the temporal data of drug release by applying a derived form of the Fick's second law outlined by Fu and Kao. 57

Model validation
To evaluate the ability of our CG model to reproduce the local molecular structure of chitosan chains in solution and their tendency to aggregate, we initially compared the radial distribution functions (RDFs) of chitosan chains with DP = 16 obtained from atomistic and CG simulations. Collectively, the RDFs between all pairs of CG sites, provide information on short-range molecular structure, represented by the position and magnitude of the peaks at short distances, as well as the overall aggregation trends of the chitosan chains, visible in the long range behavior of the curve. Because the structural data was not used in the fitting procedure for non-bonded interactions, comparing the RDFs represents a real test for the quality of the model. To characterize the conformation of single chitosan chains, angle distributions and end-to-end distances were calculated.
First, the RDFs obtained from CG simulations of chitosan chains with low degree of modification (w Ac and w But = 16%, and w Hep = 8%) aligned well with the corresponding RDFs obtained from atomistic simulations. There is an excellent agreement between CG and all-atom data for the RDFs between water beads and the various CG sites of chitosan chains, as shown in Fig. 3. These monomer-water interactions are the most sensitive to long range perturbations and sampling issues, as shown in a previous study. 25 The RDFs for all other pairs of CG sites obtained for modified chitosan chains with low w from atomistic and CG simulations are reported in Fig. S1-S3 (ESI †). The small differences observed in the short-range structure of some of the RDFs are unlikely to have an effect on the larger scale structure of the networks. The simulated acetyl-chitosan system at w Ac = 100% (Section 3.2.2) demonstrates in fact that local details can also be captured remarkably well. The largest differences in overall aggregation trends (i.e., long range tails in the RDFs) are observed in butanoyl-chitosan chains, for which all CG sites are slightly overhydrated as compared to their all-atom counterparts. This means that the polymer-polymer aggregation observed in the simulated butanoyl-chitosan networks may be somewhat underrepresented. We note, however, that the differences are small and other factors such as water content and experimental control over the modification pattern are likely to have a much larger effect on the system, when comparing in silico vs. experimental data.
The RDFs for the hydrophobic beads derived from CG and atomistic simulations at higher w (Fig. 4, black and red lines) deviate from their atomistic counterparts more than those obtained at lower w. For example, the RDF of M-M (acetyl-toacetyl) beads obtained from the CG simulation of acetylated chitosan with w Ac = 32% exhibits a major peak at short distances (o6 Å, Fig. 4a), indicating an over-representation of the attractive interactions between the acetyl groups (M-M).
Similarly, the CG RDFs of MA-MA distances (Fig. 4b) differ substantially from their atomistic counterparts at short distances (o6 Å) and show close contacts (o2.5 Å) that are physically unrealistic, and the close-range MA-MA interactions of heptanoyl-chitosan with w Hep = 16% are over-represented by the CG model (Fig. 4d). Notably, the atomistic RDFs contain several irregular peaks at all distances, which suggests the formation of stable clusters as a result of the strong interactions between the modification groups, which is likely to prevent sufficient sampling of atomistic RDFs and forces at all distances; these sampling issues for the MA-MA interactions were confirmed by comparing the RDFs obtained in separate simulation runs (Fig. S4, ESI †). The differences between atomistic and CG results are more moderate for the RDFs of MB-MB distances in butanoyl-chitosan with w But = 32% ( Fig. 4b and c), and the RDFs of MB-MB, and MC-MC distances in heptanoylchitosan with w Hep = 16% ( Fig. 4d-f).
It is anticipated that the formation of kinetically trapped clusters and the ensuing limitations in adequately sampling the conformational space of the system may lead to nonrepresentative interaction potentials because the all-atom forces do not represent the ensemble average. To avoid these issues, we tested the performance of CG interaction potentials that were obtained using chitosan chains with lower degree of modification (w Ac and w But = 16%, and w Hep = 8%) and transferred them to systems with high degrees of modification (w Ac and w But = 32%, and w Hep = 16%). A similar approach was demonstrated in a previous study, 25 where the potentials obtained with a CG procedure were found to be transferable to systems at higher concentrations. The RDFs obtained with the transferred potentials (blue lines in Fig. 4) demonstrate a much better agreement with atomistic RDFs relative to the initial CG simulations (red lines) in modeling the M-M interactions for acetyl-chitosan (Fig. 4a), and remove the physically unrealistic close-contacts from the MA-MA interactions for w But = 32% (Fig. 4b). Notably, the RDFs of the MA-MA interactions for both butanoyl-chitosan with w But = 32% (Fig. 4b) and heptanoyl-chitosan with w Hep = 16% (Fig. 4d) using the transferred potentials no longer exhibit the irregular peaks found in the all-atom systems, and closely resemble the RDF of the M-M interactions in acetyl-chitosan (note: because heptanoyl modifications are more hydrophobic and the effects of chain clustering are already observed w Hep = 8%, the interaction potentials were obtained for w Hep = 8% and transferred to w Hep of 16% and 24%). Finally, the RDFs of MB-MB and MC-MC interactions remained unchanged by the transferred potentials and in overall agreement with their atomistic counterparts (Fig. 4c, e, and f). Similarly, the RDFs between all other CG sites remain unchanged when the transferred interaction potentials were used ( Fig. S5-S7, ESI †). Collectively, these results demonstrate that the transferred interaction potentials (i.e., from lower to higher values of w) improve the agreement between atomistic and CG simulations in modeling the interactions between hydrophobic sites, while preserving the agreement between CG and atomistic approaches for all other CG sites. Accordingly, in all the simulations presented in this study, we employed the CG interaction potentials obtained at w Ac = 16%, w But = 16%, and w Hep = 8% to model systems at higher w. Notably, this approach confers the additional advantage of rendering the model more versatile for constructing polysaccharides with different modification patterns.
We further investigated the effect of chemical modification of chitosan chains on their conformation by evaluating angle distributions and end-to-end distances. The conformational space of single chitosan chains is governed predominantly by the glycosidic bonds between GlcN monomers, whose spatial arrangement is described by the dihedral angles j and c. In the CG model presented here, the conformations of the j and c angles are reflected by the two angles The probability distributions of these angles sampled in the atomistic and CG simulations are reported in Fig. 5a-c, together with snapshots of the molecular conformations corresponding to the three main minima in the j-c free energy landscape ( Fig. 5d-f). Both atomistic and CG models returned angle distributions showing a bimodal distribution ( Fig. 5a and b).
angle, the atomistic model returned a distribution with three maxima (À1401, 301, and 1401), corresponding to the three conformations shown in Fig. 5d-f, only two of which were reproduced by the CG model. The comparison between atomistic and GG distributions shows that the CG model provides a satisfactory representation of the dominant conformation for all three angles, although the states corresponding to the second energy minima (i.e., the maxima of probability distribution) are underrepresented in the distributions. This is due to the application of the Boltzmann inversion method without further iteration, which does not account for the effects of neighboring bonds and can therefore over-represents the stiffness of the interaction potentials. On the other hand, the second minimum provides a relatively minor contribution to the overall polymer conformation, as reflected for example by the end-to-end distances of the molecules.
The comparison of the distribution of end-to-end distances sampled by atomistic and CG models is shown in Fig. 6. The polymer conformational space of acetyl-chitosan is captured very well by the CG model (Fig. 6a). For butanoyl-and heptanoyl-chitosan, the impact of the more rigid glycosidic links is slightly larger. This derives from the conformations adopted by the glycosidic links immediately adjacent to the substituted monomers, which show significantly shifted angle distributions (Fig. S8, ESI †). As a result, more compact conformations are more frequently sampled with butanoyl-and heptanoyl-chitosan chains ( Fig. 6e and f). These compact states are underrepresented in the CG model, which employs the same bonded potentials for all glycosidic bonds, as shown by the distributions of end-to-end distances ( Fig. 6b and c). Nonetheless, there is a considerable overlap between the atomistic and CG distributions, indicating that the conformational space of the chains is for the most part well represented in the  CG model. It is however possible that the end-to-end distances observed in butanoyl-and heptanoyl-chitosan hydrogels are slightly overestimated, and that the transition to more compact polymer conformations, as observed in butanoylchitosan ( Fig. 6b and e), may be initiated already at slightly lower w But . Next, we proceeded to test the effect of degree of polymerization (DP) and water/chitosan ratio (W/C) on the chitosanchitosan and chitosan-water interactions. Initially, the transferability of the interaction potentials from concentrated systems (low W/C, namely 12 water molecules per GlcN monomer) to diluted systems (high W/C, namely 32 water molecules per GlcN monomer) was tested by comparing atomistic and CG RDFs for chitosan chains with DP = 16. As observed in similar systems, 25 the RDFs obtained for the diluted system with the transferred interaction potentials closely resemble those from the atomistic simulations (all RDF plots are reported in Fig. S9 and S10, ESI †). Based on these results, a larger system comprising acetyl-chitosan chains (w Ac = 16%) with DP = 50 and 32 water molecules per GlcN monomer was constructed and simulated. The resulting RDFs, reported in Fig. 7, indicate that the chitosan-water and chitosan-chitosan interactions are minimally affected by the value of DP. Only the peaks corresponding to the short-range interactions between bonded neighbors, such as those observed in the A-A and A-M RDFs ( Fig. 7c and d), increased at higher DP, due to the presence of more bonded neighbors.

Structure of the chitosan hydrogels
The CG model was applied to study the effects of modification type and w, water/chitosan (W/C) ratio, and distribution (pattern) of the modification groups along the chitosan chains upon the morphology of the hydrated network of modified chitosan chains (hydrogel). The resulting structures were quantitatively characterized by measuring the distribution of pore sizes, the end-to-end distances of the chains, and the average number of contacts between the modification groups.
3.2.1. Simulations of modified chitosan hydrogels with low water/chitosan ratio. An initial set of simulations was performed of concentrated networks comprising 50 modified chitosan chains (DP = 50) and 80 000 water molecules (32 waters/monomer). The acetyl, butanoyl, or heptanoyl modifications were distributed along the polymer chains either in even pattern (Fig. 8) or blocky (Fig. S11, ESI †) pattern. The resulting networks, snapshots for which are presented in Fig. 8a, b, d, e, g and h, feature chitosan chains uniformly distributed throughout the simulation box. The visual similarity between the network structures was confirmed by the distribution of pore sizes in these networks. The pore size distributions (Fig. 8c, f, and i) calculated for acetyl-, butanoyl-, and heptanoylchitosan networks were found to be nearly independent of w, and feature a relatively narrow distribution with most pore diameters in the range of 1.0 to 1.5 nm. Despite the overall similarity, however, the simulation snapshots suggest that the modifications with higher hydrophobicity (i.e., butanoyl and heptanoyl) tend to aggregate into local hydrophobic clusters.
To obtain a more detailed picture of the local molecular interactions within the networks, we measured the number of contacts with distances below 0.6 nm between the modification groups (Table 1). Specifically, the contacts between (i) all M beads in acetyl-chitosan networks, (ii) between MB and all MA or MB beads in butanoyl-chitosan networks, and (iii) between MA and any MA, MB or MC beads in heptanoyl-chitosan networks were counted. As anticipated, the number of contacts was found to increase with the hydrophobicity of the modification group and w. In acetyl-chitosan systems with evenly-spaced modification pattern, the average number of contacts between M beads grows from a minimum of only 0.07 AE 0.02 per modified monomer when w Ac = 16% to 0.22 AE 0.02 when w Ac = 50%. Similarly, the contacts formed by the butanoyl MB beads increases from 1.08 AE 0.02 when w But = 16% to 1.25 AE 0.01 when w But = 32%. Finally, the heptanoyl MC beads form on average 1.53 AE 0.006 contacts per monomer when w Hep = 8%, and 1.63 AE 0.015 contacts when w Hep = 24%.  Together with the contact numbers, the values of end-to-end distance provide insight into the role of chain modification and intermolecular interactions on the conformation of chitosan. The variation in average end-to-end distance with different modification groups, w, and pattern (i.e., evenly spaced or blocky) is presented in Table 2. A consistent trend is observed, where acetylation increases the end-to-end distances, and both butanoation and heptanoation drive a sharp decrease. The amide bond introduced by acetylation, combined with its mild hydrophobicity, promotes the formation of chitosan-water and chitosan-chitosan hydrogen bonds. This in turn drives the alignment of chitosan chains, whether single or in a network, which translates in higher end-to-end distances. Butanoyl and heptanoyl groups, instead, impart a higher hydrophobicity, which promotes intramolecular monomer-monomer interactions and chain de-solvation. These results demonstrate that the modification type and w can be determining parameters for the structure of the chains, and larger changes to the network structure are likely hindered predominantly by the dense packing in the networks.
3.2.2. Simulations of modified chitosan hydrogels with high water/chitosan (W/C) ratio. The hydration of chitosan chains (i.e., W/C ratio) is a major determinant of the network morphology, as it controls chain flexibility and intermolecular interactions. The water content of real hydrogels can reach values of B99% w/w and is typically larger than what is feasible in silico. To investigate the effect of W/C ratio on the network structure, we repeated the simulation of the diluted chitosan systems (20 chains with DP = 50) listed in Table 2, yet with a much higher solvation. Specifically, we utilized 100 000 CG water beads (i.e., 100 water molecules per GlcN monomer), corresponding to a water weight fraction of 90% w/w. Representative snapshots collected from the MD simulations of chitosan chains with even modification and the corresponding distribution of pore sizes are reported in Fig. 9, while the number of contacts between CG beads are listed in Table 3 (note: as the W/C ratio increases, the number of M beads and M-M contacts decrease; therefore, the number of contacts normalized by the number of modifications is also reported to enable a better comparison with the analogous results for more concentrated systems reported in Table 1). The corresponding results for chitosan chains with blocky modifications are in Fig. S12 (ESI †).
These results indicate that the overall structure of the networks formed by acetylated chitosan chains remains independent of w (Fig. 9a-c), as was observed in systems at the higher chitosan concentration (Fig. 8). The acetyl-chitosan chains are evenly distributed throughout the simulation box, and the pore-size distributions are independent of w and feature pore sizes between 1.5-2.5 nm.
The butanoyl-chitosan systems, on the other hand, show a drastic restructuring of the chain network as w But increases from 16% to 32%. This triggers the formation of dense ''micelle-like'' structures, comprising clusters of butanoyl moieties surrounded by chitosan backbones, connected by individual polymer strands (Fig. 9e). The transition of the network morphology from homogeneous to cluster/channel is accompanied by a distinct increase in the measured pore size, from 1-3 nm to 4.0-6.0 nm (Fig. 9f). While partly caused by the larger size of the simulation box, the observed variations of pore size are an important indicator of the morphological change of the network. The value of pore size in the larger systems are determined both by the water content and the length of the polymer bridges.
Reflective of this morphological change, a notable decrease in the end-to-end distances from 12.7 AE 0.8 nm (w But = 16%) to 5.6 AE 0.4 nm (w But = 32%) was observed in this network (Table 2), which results from the chitosan backbones wrapping the butanoyl clusters. A closer inspection of the cluster-channel structure (Fig. 10) indicates that the clusters comprise several hydrophobic domains (or centers) with diameters of B1-1.5 nm, corresponding to approximately twice the length of the hydrophobic moieties. These domains coalesce into larger clusters featuring almost completely hydrophilic surfaces formed by layers of chitosan chain backbone. The persistence length of chitosan chains, which experimental measurements estimate to be in the range of 5-20 nm 58-60 and our simulations indicate it to be B5 nm (based on the end-to-end distance, Table 2), is comparable to the separation between adjacent modification groups on the chains. Forming loops at this scale would lead to a considerable cost in bending energy. Thus, separate hydrophobic moieties of one polymer chain do not typically cluster in one hydrophobic domain, but are rather embedded in separate hydrophobic domains that can either be part of the same larger cluster (Fig. 10, orange) or as bridges between two larger clusters (Fig. 10, red). In both cases, the modified chitosan chains remain approximately linear on the length scale of several monomers (i.e., their persistence length). Both micelles with intermolecular hydrophobic domains, resembling the clusters observed here 61 and the bridged micelles 62 have been proposed as probable models for the structure of hydrophobically modified chitosan. Interestingly, no morphological transition took place in heptanoyl-chitosan systems, despite the higher hydrophobicity of the heptanoyl groups. While small clusters of M beads can be observed in the simulation snapshots, the pore size distribution remains approximately constant with w Hep . Consistent with the absence of hydrophobic clusters, no significant decrease in end-to-end distances was observed ( Table 2).
This different behavior of butanoyl-and heptanoyl-chitosan networks stems from the association between the hydrophobic modification groups, which is much stronger with heptanoyl moieties. Accordingly, the initial heptanoyl clusters are too Fig. 9 Simulation snapshots of networks comprising 20 chitosan chains with DP = 50 and 100 water molecules/monomer, and modified with acetyl groups at (a) w Ac = 16% and (b) w Ac = 50%; butanoyl groups at (d) w But = 16% and (e) w But = 32%; heptanoyl groups at (g) w Ac = 8% and (h) w Ac = 24%. All modifications were distributed evenly on the chitosan chains. The chitosan backbone (A, B, and C beads) is shown in red, the modifications (M, MA, MB, and MC beads) are in yellow, and water molecules are in blue, and the corresponding pore-size distributions in (c) acetyl chitosan, (f) butanoyl chitosan and (i) heptanoyl chitosan. stable to allow a rearrangement of the network into a cluster/ channel morphology, which requires dissociation and rearrangement of some of these initial contacts, i.e., the uniform networks are kinetically trapped. While it is possible that network structures appear to be kinetically trapped due to insufficient simulation time, the experimental values of drug diffusion throughout heptanoyl-chitosan systems (vide infra) indicate that similar effects are present at the experimental time scale. Furthermore, the numbers of contacts per modification bead in diluted networks formed by butanoyl-chitosan chains with w But = 32% and heptanoyl-chitosan chains with w Hep = 24% (Table 3) are nearly identical to the corresponding values found in concentrated networks (Table 1). This showcases the role of hydrophobic modifications in forming contacts among modified chitosan chains.

Simulations of modified chitosan chains with different modification patterns.
Together with the degree of modification (w), the modification pattern (i.e., the distribution of the modification groups along the chitosan chains) is another important design parameter. Polymers with the same global composition but varying modification patterns, in fact, differ significantly in terms of conformation and aggregation of the chains, and ultimately morphology of the hydrogel network. 63,64 Unlike w, which can be easily varied and measured, the modification pattern is harder to control and characterize experimentally. [65][66][67] Different modification patterns can be achieved by performing the chemical modification on the polymer chains dissolved in solvents of different ''quality'' and correlated to the Kerr constant of the resulting modified polymer in solution. [68][69][70][71][72][73] In silico modeling provides full control over the distribution of modification groups, and represents an ideal modality for exploring the effect of modification pattern on the molecular architecture of hydrogels.
A random distribution of acyl groups on the GlcN monomers is generally assumed for modified chitosan and has been supported by NMR analysis of chitosan samples with different degrees of acetylation (w Ac ). 65 However, the acylation pattern of chitosan chains can vary; for example, several studies have shown that deacetylation of chitin under heterogeneous conditions can result in chitosan with blocky modification patterns. [74][75][76] The use of hydrophobic groups for chitosan modification can also promote the formation of blocky patterns: hydrophobic moieties initially installed on the chain act as ''seeding points'', as they attract incoming modification reagents by hydrophobic interaction and promote their local conjugation to the chains. The local accumulation of modification groups around hydrophobic focal points translates into the formation of blocky patterns. The aliphatic chains introduced by butanoic and heptanoic anhydrides can lend themselves to this mechanism.
In this study, we investigated the effect of modification patterns using two model distributions, evenly spaced and regular blocky. Blocky patterns in real polymers likely comprise a distribution of block sizes; similarly, truly random patterns also comprise short blocks of different length. Constructing realistic models of irregular blocky and truly random polymers in silico is not only challenging from the point of view of design, but can also return equivocal results. Therefore, to elucidate the role of modification pattern on chain structure and network morphology, we resolved to compare two idealized modification patterns: (i) ''evenly spaced'' with single modified monomers distributes at equal distances along the polymer and (ii) ''blocky'', where four consecutive monomers bear modifications (Fig. 11a). In the comparison between these different patterns, the largest influence on the network morphology is, once again, observed in butanoyl-chitosan. The networks of butanoyl-chitosan chains, featuring identical w But = 16% and DP = 50 and blocky vs. evenly spaced modification, shown in Fig. 11b, illustrates the effect of pattern on chain aggregation and network morphology. A systematic investigation of the combined effect of modification type, w, and pattern on the number of contacts and end-to-end distances of modified chitosan chains (Tables 1-3) and the pore size distribution (Fig. 11d-i) of the resulting structures indicates that the acetylchitosan networks do not possess any structural dependence upon modification pattern. Acetyl-chitosan features a two-fold helical structure, where neighboring modifications lie on alternating sides of the chain. 77 Therefore, an increase in the number of contacts can occur only if the modifications aggregate into local clusters. The constant number of contacts in both concentrated and dilute acetyl-chitosan systems (Tables 1 and 3) indicates that such clustering at low water/chitosan ratio does not depend on the modification pattern. Similarly, chainchain interactions and network morphology in concentrated butanoyl-chitosan networks (low W/C) do not present any correlation to modification pattern (Table 1). A different scenario appears in diluted (high W/C) butanoyl-chitosan systems with lower w But (16%), where chain aggregation is promoted by the blocky patterning (Fig. 11e); at w But = 32%, on the other hand, the cluster/channel architecture of the Fig. 10 Simulation snapshot of the cluster-channel morphology in butanoyl-chitosan system with w But = 32%. The butanoyl groups are highlighted in yellow, whereas the polymer backbone is in grey. The conformations of three individual polymer chains are highlighted: one that forms part of a single ''cluster'' (orange) and two that form bridges between separate clusters (red). network and its characteristic higher pore sizes are formed irrespective of the modification pattern (Fig. 11h). Finally, in heptanoyl-chitosan systems, the number of contacts is consistently higher between chains with blocky modification patterns (Tables 1 and 3); the stronger tendency of blocky heptanoylchitosan towards the cluster/channel architecture is also indicated by the small, yet notable, increases in pore size for both low (w Hep = 8%) and high (w Hep = 24%) levels of modification ( Fig. 11f and i). Nonetheless, this effect is much less pronounced in heptanoyl-chitosan systems than in butanoyl-chitosan systems. The strong hydrophobic interactions between heptanoyl groups, in fact, make the chain network more rigid and prevents its rearrangement into the cluster/channel morphology.
3.2.4. Simulations of full chitosan (v Ac = 0%) and chitin (v Ac = 100%) networks. To evaluate the predictive power of our model, we simulated two limit cases, namely w Ac = 0% (known as ''full chitosan'') and w Ac = 100% (chitin). While chitin is commercially available, fully de-acetylated chitosan is difficult to produce at large scale, and commercial forms of chitosan typically feature w Ac = 5-25%. Experimental studies report that full chitosan forms hydrogels with homogeneous networks. 78 Similarly, our CG simulations indicate that chitosan chains with w Ac = 0% produce a uniform network spanning the entire simulation box, as shown by the snapshot in Fig. 12a, the pore size distribution in Fig. 12b, and a value of end-to-end distance (14.7 AE 0.65 nm) that resembles those found in acetyl-chitosan networks at low w Ac (Table 2).
Consistent with experimental observations, our CG simulation was also able to predict the formation of thick chitin fibrils consisting of bundles of aligned single chitin chains (Fig. 13a);  coherent with chain alignment, the values of end-to-end distance of the aligned polymers increases to 20.4 AE 0.4 nm, as compared to B15 nm in the gels. Chitin forms crystalline fibrils wherein the alignment of the chains is either anti-parallel (i.e., a-chitin 79,80 ) or parallel (i.e., b-chitin 81 ). Both chitin allomorphs are found in natural chitin, although b-chitin has never been recrystallized in vitro. In the CG simulations, in fact, half of the chitin chains in every fibril form the anti-parallel alignment observed of a-chitin, while the other chains are aligned in parallel as in b-chitin (Fig. 13b), determined predominantly by the assembly kinetics. The simulated fibrils, however, are more strongly twisted compared to the reported crystal structures.
Both effects are likely due to the onset and evolution of the aggregation process in silico, where the outcome of fiber alignment is strongly dependent on the first contacts formed. Nonetheless, the molecular structure and intermolecular contacts in the segments with anti-parallel alignment produced by the CG simulations show remarkable similarity with those detected experimentally in a-chitin 79,80 (Fig. 13c). At the same time, reproducing fully the crystal fibers and the internal orientation of the chains in such large systems is beyond the capability of unbiased simulations, which are kinetically constrained. Collectively, these results demonstrate the ability of the proposed CG model to capture the molecular aggregation and simulate the architectural features of chitosan-based hydrogels across the entire range of w Ac . Remarkably, the model maintains reasonable predictive accuracy for values of w Ac that differ significantly from those upon which the CG interaction potentials were developed.

Correlating hydrogel properties to the diffusion of probe molecules
Chitosan hydrogels have been extensively utilized as materials for tissue engineering, wound healing, and drug delivery. [1][2][3]13,14,82 Their optimization as drug delivery carriers for clinical applications, however, is extremely laborious due to the width of the experimental space of design controlling drug loading and release. 56 Computational models capable of predicting the transport of drugs through chitosan-based hydrogels hold great promise to accelerate pre-clinical development and ultimately improve translation potential. Key to this end is the ability of the model to connect physicochemical properties and architecture of the chain network and drug properties. The characteristic size and complexity of these systems prompts the adoption of CG models, whose smoother energy landscape typically enables faster systems dynamics 15,16 and therefore a more efficient investigation of the dependence of drug delivery upon a wide space of design.
In this work, we adopted the synergistic chemotherapeutic pair doxorubicin (Dox) and gemcitabine (Gem) as model molecules to study diffusion through modified chitosan hydrogels. Gem is a small (263.2 g mol À1 ), hydrophilic and electrically neutral species, whereas Dox is larger (543.5 g mol À1 ) and amphiphilic. CG models of Dox and Gem were initially prepared through the mapping shown in Fig. 1, while the drug-drug, drug-solvent, and drug-chitosan interactions were obtained from initial all-atom simulations following the same steps as for the chitosan interactions. Initially, we calculated the diffusion coefficients of the two drugs in water using atomistic and CG simulations, and compared them with experimental data (Table 4) to characterize the dynamics in our CG systems. The value of water self-diffusion coefficients obtained from atomistic simulations (D AS W,W ) agrees with the diffusion coefficients reported in the literature for the TIP5P water model 83 and is similar to diffusion coefficients determined experimentally (D Exp W,W ). 84 As anticipated, the water diffusion coefficient provided by CG simulations (D CG W,W ) is larger, resulting in a scaling factor (t W = D CG W,W /D AS W,W ) of 6.4. Similarly, scaling factors of 9.6 and 10.0 were found for Gem and Dox, respectively, between the values of diffusion coefficients derived from atomistic vs. CG simulations. We note, however, that Dox molecules tend to form small aggregates, both in atomistic and CG simulations, which may affect the observed scaling factor; a single molecule in a box of water, in fact, has a 10-fold Therefore, upon rescaling, the CG model provides a reasonable prediction of the diffusion coefficients of the model drugs.
To model diffusion through modified chitosan hydrogels, 20 drug molecules -either Dox or Gem -were randomly placed inside an equilibrated network comprising 20 chains (DP = 50 and 5000 water molecules per chain) of either (i) acetyl-chitosan (w Ac = 16% or 50%), (ii) butanoyl-chitosan (w But = 16% or 32%), or (iii) heptanoyl-chitosan (w Hep = 8% or 24%) featuring either evenly-spaced or blocky modification patterns. The diffusion of Dox and Gem through the networks was analyzed by tracking the mean-square displacement (MSD) of the drug molecules vs. time, from which the corresponding diffusion coefficients were calculated. In parallel, following the experimental procedure developed in prior work, 56 the above-listed hydrogels were constructed and loaded with either Dox or Gem at the same concentration as utilized in the simulations, and the drug release was measured as a function of time. From the resultant drug release profiles, the values of diffusion coefficients were calculated by applying a derived form of the Fick's second law for slab-like devices. 57 The in silico and experimental results are collated in Table 5.
The CG simulations indicate that the diffusion of Gem through modified chitosan networks D CG Gem,Chit is not affected by the modification type, w, and pattern, and is only slightly lower than the corresponding value in water (D CG Gem,W , Table 5). This was confirmed by the experimental values of D Exp Gem,Chit . With its hydrophilic character and small size (R gyr = 0.33 nm) compared to the average pore diameter of the simulated networks (1.2 nm), Gem can migrate seamlessly through the chitosan hydrogel. The absence of strong interactions with the chitosan networks is corroborated by the low number of Gem-chitosan contacts ( Table 6).
On the other hand, with a larger size (R gyr = 0.52 nm) and a distinctively higher hydrophobicity compared to Gem, Dox migrates more slowly. The diffusion coefficients of Dox obtained from CG models were found to be lower than the reference value in water (D CG Dox,Chit o D CG Dox,W ) and to vary significantly with the type, w, and pattern of modification. Specifically, Dox diffusion in acetyl-chitosan decreases moderately as w Ac increases from 16% to 50%, whereas it increases more than 3-fold in butanoyl-chitosan between w But = 16% and 32%, and 1.5-fold in heptanoyl-chitosan between w Hep = 8% and 24%. The same trends were observed in the experimental values ( Table 5).
The interactions with the chitosan network are reflected by the average number of Dox-polymer contacts (distance o0.6 nm), which account for both the hydrophobic interactions with the modified monomers (GlcNAc, GlcNBut, and GlcNHep) and the hydrogen bonds with the chain backbone (Table 6). Notably, the number of contacts between Dox and the acetyl groups (M beads) nearly triples as w Ac increases from 16% to 50% in chitosan chains with evenly-modified pattern. Concurrently, the number of contacts with the backbone increases; as Dox molecules are attracted towards GlcNAc beads by hydrophobic interactions, they also form secondary contacts with the intercalated GlcN beads via hydrogen bonding. 56 A different trend is observed with blocky acetyl-chitosan, where the number of interactions between Dox and M beads grows with w Ac , whereas the number of contacts with backbone beads decreases; the blocky modification pattern, in fact, by removing the intercalation between the monomers, lowers the probability of Dox forming secondary interactions with GlcN as a consequence of initial hydrophobic binding to GlcNAc monomers.  The interactions between Gem and acetyl-chitosan chains, instead, is more straightforward. The number of contacts with the M beads is low and independent of the modification pattern; the increase with w Ac is related to the formation of hydrogen bonds with the amide moieties connecting M and A beads. 56 The number of contacts with the backbone beads is higher and rather unaffected by either w Ac or pattern. The observed sharp increase in Dox diffusion coefficient with w But can be attributed to the formation of the large pores in the network following the transition from uniform structure to the cluster/channel morphology. With larger pores, the Dox molecules can diffuse freely through the aqueous medium without interacting with the butanoyl-chitosan chains (Fig. 9e), thereby increasing their mobility and ultimately their diffusion coefficient. On the other hand, it is also noted that Dox forms a higher number of contacts with both the backbone and the modification beads compared to that formed with the acetyl-chitosan chains ( Table 6). The simulation snapshots of Dox diffusing through networks with different w But and pattern, in fact, show that several Dox molecules adsorb onto -and become irreversibly embedded into -the chain clusters. This observation was confirmed by analyzing the mean squared displacement (MSD) of individual Dox molecules (Fig. 14).
These fall into two distinct groups: one moving with normal Brownian diffusion (MSD p DÁDt a , a = 0.98) and the other exhibiting a significantly slower dynamics (a B 0.7). All the diffusion curves of Dox through the uniform networks fall between these two behaviors, with a comprised between 0.76-0.82; all the individual atoms in these systems show the same diffusion characteristics (Fig. S14, ESI †), consistent with the uniform structure of the networks. Furthermore, the van Hove functions (Fig. S13, ESI †) indicate some confinement at the length scale of B5 nm in the uniform networks; however, in the systems formed by butanoyl-chitosan chains with w But = 32%, a small fraction of Dox molecules travel farther on a timescale of 8-10 ns than in any of the uniform networks. Similar larger displacements are also seen for heptanoylchitosan chains with w Hep = 24%, indicating that these networks are less uniform than was suggested by the calculated distribution of pore sizes.
Accordingly, the increase in D CG Dox,Chit with w But resulting from the MSD analysis is uniquely due to unimpeded transport of unabsorbed Dox molecules through the pores of butanoylchitosan cluster/channel networks. So that the diffusion and release of Dox from hydrogels with such cluster/channel morphology will depend on the numbers of molecules diffusing through the pores, N p , and being adsorbed to the clusters, N cl . These are determined by the partition coefficient where C cl and C p are the concentration of drug molecules adsorbed on the clusters and contained in the pores, and V cl and V w are the volumes of the chitosan clusters and water-filled pores. The number of molecules stuck to the clusters is thus determined both by the affinity of Dox for the hydrophobic moieties and the volume ratio of clusters in In hydrogels with high water content (V w c V cl ), as is the case for the modified chitosan hydrogels considered in this study, the majority of Dox molecules will be located in the pores and the observed (average) Dox diffusion is governed by those molecules diffusing freely through the water filled pore space. The experimental values of Dox diffusion coefficient and total release from butanoyl-chitosan hydrogels increase with w But , confirming the results of the CG model. Numerical discrepancies between simulated and measured diffusion likely stem from the difference in network hydration between in silico and the in vitro systems. Nonetheless, the CG Table 6 Average number of contacts (distance o0.6 nm) formed by Dox and Gem with either the modification beads (M in acetyl-chitosan; MA and MB in butanoyl-chitosan; and MA, MB, and MC in heptanoyl-chitosan) or the backbone beads (A, B, and C), as obtained from CG simulations of drug transport through networks of chitosan chains with different type, w, and pattern of modification model captures well the mechanisms determining the dependence of the diffusion coefficient of Dox upon w But . Finally, heptanoyl-chitosan networks present a hybrid behavior, intermediate between those of acetyl-and butanoylchitosan. The partial clustering of the heptanoyl moieties, while not producing a defined cluster/channel morphology, results in larger pores as w Hep increases, especially with blocky heptanoation (Fig. 10f and i). As a results, the values of diffusion coefficient of Dox in heptanoyl-chitosan networks fall between those calculated in acetyl-and butanoyl-chitosan systems. Despite the increase in pore diameter, in fact, Dox retention by the heptanoyl-chitosan chains is significant, as shown by the high number of hydrophobic interactions with the modification beads ( Table 6). As a result, the diffusion coefficient of Dox through heptanoyl-chitosan with evenly spaced modification increases with w Hep , as in butanoylchitosan, whereas it decreases when heptanoation is blocky, as in acetyl-chitosan. These values indicate that Dox diffusion is determined by the effective hydrophobicity of the chains and the morphology of the network, which are interconnected to the design parameters of the chitosan hydrogel (i.e., chemical modification and water fraction). The comparison between these predictions and the experimental data reported in Table 5 showcases the predictive accuracy of the proposed CG model, which stems from its ability to accurately determine diffusion coefficients while accounting for the interplay between modification-based molecular architecture and the chitosan-drug interactions.

Conclusions
Pre-clinical development of polysaccharide-based hydrogels for drug delivery applications consists of the systematic evaluation of the wide range of design parameters, such as base chemical composition and physicochemical modification of the substrate materials, the water to polymer ratio, and the pharmaceutical payload. A purely empirical approach to the exploration of such a wide design space is unfeasible in terms of both times and costs. To accelerate pre-clinical development, we have developed and demonstrated a systematic bottom-up CG model for evaluating drug transport through hydrogels constructed with chitosan chains featuring different types, degrees, and patterns of modifications. This CG model does not require any empirical data input to make predictions of the molecular architecture or the drug transport properties through the hydrogels. Numerous types of modification and therapeutic payloads, other than those explored here, can therefore be integrated into this model by implementing the CG procedure outlined herein.
This study demonstrates that the proposed CG model (i) is flexible and can be transferred to large systems comprising many long polysaccharide chains and different water contents; (ii) accurately captures the supramolecular structures (e.g., clusters/channels, fibers) and connects them to basic design parameters; and (iii) predicts very accurately the transport of small drugs through these materials by virtue of a comprehensive interpretation of the fundamental physicochemical properties and molecular interaction mechanisms. In so doing, this model provides a study toolbox with two-fold significance. In terms of fundamental science, our model enables the use of simple measurements of diffusion of small organic molecules (e.g., fluorescent Dox molecules) in lieu of laborious and expensive experimental techniques (e.g., small angle neutron or light scattering) to evaluate the molecular and supramolecular architecture of modified polysaccharide hydrogels. In terms of technological relevance, this model provides pharmaceutical chemists with a reliable rationale for focusing their developmental efforts towards a specific subset of the wide space of experimental design inherent to drug-loaded soft materials, allowing them to rapidly tailor therapeutic materials that meet the need of different therapeutic treatments, covering the entire span from first-line to consolidation chemotherapy.