A MULTI-SCALE TOOLBOX TO PREDICT
STRUCTURE AND FUNCTION OF
POLYSACCHARIDES AGGREGATES
vorgelegt von
Master of Science (M.Sc.)
Ankush Singhal
an der Fakultät II - Mathematik und Naturwissenschaften
der Technischen Universität Berlin
zur Erlangung des akademischen Grades
Doktor der Naturwissenschaften
Dr. rer. nat.
genehmigte Dissertation
Promotionsausschuss:
Vorsitzender: Prof. Dr. Michael Lehmann
Gutachterin: Dr. Andrea Grafmüller
Gutachterin: Prof. Dr. Sabina Klapp
Gutachterin: Prof. Dr. Maria Andrea Mroginski
Tag der wissenschaftlichen Aussprache: 18th February 2020
Berlin 2020
Abstract
Carbohydrates are class of biomolecules- their functions and properties covers a vast field
that still needs to be explored. Many biological polysaccaharides form aggregates and
their structures and properties are very versatile and depends on the aggregate structure
and molecular interactions. Natural polysaccahrides can also form hydrogels, porous
network of polymers, that can take up a high percentage of water. Their properties
can be additionally tuned by the introductions of chemical modifications to a fraction of
monomers. To make efficient use of their versatile properties, understanding the relation
between the molecular structure and interactions of polysaccahrides and the properties
of the aggregates formed is essential.
Computational modeling provides an efficient tool for understanding interactions at
the molecular level, thus providing a qualitative direction for future experiments. Hence
modeling offers a cost and time efficient method for their study. In this thesis, the
aggregates and structures formed by different glucose and chitosan oligomers were simu-
lated and the resulting solution or aggregates structures are characterized using all-atom
and coarse-grained molecular dynamics. Usually, polymers have slow dynamics making
all-atom simulation computationally inefficient. Therefore a coarse-grained model was
developed to study the properties of the polysaccahrides at the required length and time
scale.
Chitosan hydrogels with various hydrophobic modification were modeled. The trans-
ferability of short oligomers with respect to different water concentration, degree of poly-
merization and modification was explicitly established. Different morphological network
structures of longer polymer were obtained corresponding to different degree, type, and
pattern of modification. In particular, different morphological transition from a uniform
polymer network to a structure containing dense hydrophobic cluster and large pores
was found for certain conditions.
Finally, one of the principle applications of the chitosan hydrogel as a drug carrier
was explored. The molecules Doxorubicin(DOX) and Gemicitabine(GEM) were chosen
1
as model drugs and their interactions with the different modified chitosan polymers have
been thoroughly studied at all-atom and coarse-grained resolution. The diffusion of DOX
and GEM through the different network morphologies formed by the hydrophobically-
modified chitosan was found to show quite different, network dependent trends. Whereas
GEM migrates through all chitosan hydrogels freely irrespective of type and degree of
modification. Placing the drugs together in the networks affects the diffusion behavior of
both. The results demonstrate the potential of this computational tool in the systematic
development of drug-loaded hydrogels for pharmaceutical applications.
2
Zusammenfassung
Kohlenhydrate sind eine klasse von biomolekülen- ihre funktionen und eigenschaften
umfassen ein weites feld, das noch nicht erforscht ist. Viele biologische polysaccaharide
bilden aggregate und ihre strukturen und eigenschaften sind sehr vielseitig und hän-
gen von der aggregatstruktur und den molekularen wechselwirkungen ab. Natürliche
polysaccharide können auch hydrogele bilden, ein poröses netzwerk von polymeren, die
einen hohen anteil an wasser aufnehmen können. Ihre eigenschaften können durch die
einführung chemischer modifikationen an einem bruchteil der monomere zusätzlich opti-
miert werden. Um ihre vielseitigen eigenschaften effizient nutzen zu können, ist es uner-
lässlich, den zusammenhang zwischen der molekularstruktur und den wechselwirkungen
von polysaccahriden und den eigenschaften der gebildeten aggregate zu verstehen.
Die rechnergestützte modellierung stellt ein effizientes werkzeug zum verständnis
von wechselwirkungen auf molekularer ebene dar und liefert so eine qualitative orien-
tierung für zukünftige experimente. Daher bietet die modellierung eine kosten- und
zeiteffiziente methode für ihre Studie. In dieser arbeit wurden die aggregate und struk-
turen, die aus verschiedenen glukose- und chitosanoligomeren gebildet wurden, simuliert
und die resultierenden lösungs- oder aggregatstrukturen werden durch eine vollatomige
und grobkörnige molekulardynamik charakterisiert. Normalerweise weisen polymere eine
langsame dynamik auf, was die simulation von atomen ineffizient macht. Daher wurde
ein grobkörniges modell entwickelt, um die eigenschaften der polysaccharide in der er-
forderlichen länge und Zeit zu untersuchen maßstab.
Chitosan-hydrogele mit verschiedenen hydrophoben modifikationen wurden model-
liert. Die ubertragbarkeit von kurzen oligomeren in bezug auf unterschiedliche wasserkonzen-
tration, polymerisationsgrad und modifikation wurde explizit festgelegt. Verschiedene
morphologische netzwerkstrukturen aus längerem polymer wurden erhalten, die unter-
schiedlichem grad, typ und muster der modifikation entsprechen. Insbesondere wurde
unter bestimmten bedingungen ein unterschiedlicher morphologischer ubergang von einem
einheitlichen polymernetzwerk zu einer Struktur mit dichtem hydrophoben cluster und
3
großen poren gefunden.
Schließlich wurde eine der hauptanwendungen des chitosan-hydrogels als wirkstoffträger
untersucht. Die moleküle Doxorubicin(DOX) und Gemicitabin(GEM) wurden als mod-
ellmedikamente ausgewählt und ihre wechselwirkungen mit den verschiedenen modi-
fizierten chitosanpolymeren wurden gründlich in atomarer und grobkörniger auflösung
untersucht. Die diffusion von DOX und GEM durch die verschiedenen netzwerkmor-
phologien des hydrophob modifizierten chitosans zeigte ganz unterschiedliche, netzw-
erkabhängige Trends. Während GEM durch alle chitosan-hydrogele frei wandert, un-
abhängig von Art und Grad der modifikation. Das zusammenstellen der medikamente
in den netzwerken beeinflusst das Diffusionsverhalten beider. Die ergebnisse zeigen das
potenzial dieses rechenwerkzeugs für die systematische entwicklung von medikamenten-
beladenen hydrogelen für pharmazeutische anwendungen.
4
Declaration
I declare that this thesis is an original report of my research, has been written by me
and has not been submitted for any previous degree. I confirm that the work submitted
is my own, except where work which has formed part of jointly-authored publications
has been included. My contribution and those of the other authors to this work have
been explicitly indicated below. I confirm that appropriate credit has been given within
this thesis where reference has been made to the work of others.
The work presented in Chapter 2 was previously published in Systematic hydro-
gen bond manipulations to establish polysaccharide structure-property cor-
relations published in Angewandte Chemie and a manuscript Multi-scale modelling
study of self-assembly of cellulose based derivatives under preparation.
The work presented in Chapter 3 to 5 was previously published in Tailoring the
Chemical Modification of Chitosan Hydrogels to Fine Tune the Release of a
Synergistic Combination of Chemotherapeutics in ACS Biomacromolecules and
a manuscript Predicting Chitosan Hydrogel Properties with Multiscale
Coarse-Grained Simulations is about to be submitted in Soft Matter.
Ankush Singhal
Potsdam, 04 December 2019
5
Contents
1 Introduction 9
1.1 Carbohydrates ................................. 9
1.1.1 Polysaccharides: Cellulose and chitosan ............... 10
1.2 Hydrogels ................................... 12
1.3 Molecular simulation ............................. 13
1.3.1 Molecular dynamics .......................... 14
1.4 Coarse-grained simulation .......................... 17
1.4.1 Multi-scale coarse graining (Force Matching) ........... 17
1.4.2 Boltzmann inversion ......................... 19
1.4.3 Coarse-grained models for polysaccharides ............. 20
1.5 Aim and overview of the thesis ....................... 20
2 Tailor-made cellulose derivatives 22
2.1 Introduction .................................. 22
2.2 Methods: MD details ............................ 24
2.3 All-atom simulation results .......................... 25
2.3.1 Glucose and its derivatives ...................... 25
2.3.2 Chitosan ................................ 30
2.4 CG model and mapping ........................... 32
2.4.1 CG simulations for short oligomers ................. 32
2.5 Model validation ............................... 33
2.5.1 Cellulose ................................ 33
2.5.2 Methyl modified cellulose ....................... 34
2.5.3 Amine (-NH2and -NH+
3) groups .................. 36
2.5.4 Single polymer conformation ..................... 37
2.6 Aggregation and aggregate properties .................... 38
6
2.6.1 Cellulose and chitin structure .................... 39
2.6.2 Methylated and fluorinated cellulose hydrogels ........... 40
2.6.3 Chitosan hydrogel ........................... 43
2.7 Conclusion ................................... 45
3 A multiscale model for hydrophobically modified chitosan 46
3.1 Introduction .................................. 46
3.2 Methods .................................... 47
3.2.1 All-atom simulations ......................... 47
3.3 CG Model ................................... 48
3.3.1 CG Simulations ............................ 50
3.4 Model validation ............................... 50
3.5 CG model for drug .............................. 59
3.6 Conclusion ................................... 62
4 Effect of hydrophobic modifications on chitosan hydrogel properties 63
4.1 Introduction .................................. 63
4.2 Methods: CG network model ........................ 64
4.3 Low water content .............................. 65
4.4 High water content .............................. 69
4.5 Influence of modification pattern ...................... 71
4.6 Full range of substitution: 0% and 100% .................. 76
4.7 Conclusion ................................... 78
5 Effect of chitosan hydrogel properties on drug diffusion 79
5.1 Introduction .................................. 79
5.2 Analysis of drug-chitosan interactions at all-atom resolution ....... 80
5.3 Dynamic properties of the CG system .................... 84
5.4 Simulation of single drug migration through modified chitosan hydrogels:
Lower drug concentration .......................... 85
5.4.1 Effect of higher drug concentration on diffusion trends ...... 92
5.5 Simulation of dual drug migration through modified chitosan hydrogels . 94
5.6 Conclusion ................................... 97
6 Conclusion 98
7
Chapter 1
Introduction
1.1 Carbohydrates
Carbohydrates are a class of biologically significant compounds consisting of chemically
bonded carbon, hydrogen, and oxygen, that, together with lipids, proteins and nucleic
acids, belong to the fundamental building blocks of living systems. Their biological
functions reach from structural stability1and energy storage2, 3 to cell communication4, 5
and interactions with bacteria and viruses. Some major examples of carbohydrates
include sugar (glucose), starch6and cellulose6. Based on the number of sugar units,
carbohydrates can be classified into sub-categories, namely mono-, di-, oligo- and poly-
saccharides. Polysaccharides are the most abundant organic matter on the earth. These
have applications in textile manufacture7, food8, paper9and pharmaceutical industry10
.
Despite their versatile biological functions, carbohydrates remain the least studied and
understood class of biomolecules, owing to their complexity and a lack of appropriate
experimental methods to study them systematically.
The basic building blocks of polysaccharides are monosaccharides, as the examples
shown in Figure 1.1 a) and b). Most biologically relevant monosaccharides are pentoses
or hexoses, containing five or six carbon atoms, respectively. They form rings in which
the carbon atoms are typically labeled C1-C6 starting from the ring oxygen and count-
ing in the clockwise direction. Monosaccharides can be classified as two stereoisomers,
namely Levorotatory(L) and Dextrorotatory(D). They basically used to describe the ro-
tation of plane polarized light around the chiral center. Each asymmetric carbon, e.g.
linked to an OH group, is a chiral center. Based on the conformation of the -OH group
9
Figure 1.1: Chemical structures of (a). β-D- glucose, (b). β-D-Acetyl-glucosamine
attached to the anomeric carbon(C1) are defines the αor βforms of glucose. Here, α
correspond to equitorial position of -OH with respect to ring oxygen while -OH will be in
an axial position in βform. In Figure 1.1,βforms of glucose and N-acetyl-glucosamine
are shown.
Monosaccharides can be connected to form polymers by glycosidic bonds, covalent
bonds that form between the anomeric carbon and any hydroxyl group of the next
monomer. Due to the many possible configurations, the resulting molecular structures
range from linear to highly branched and from oligomers consisting of a few monomers
to polysaccharides with a degree of polymerization (DP) of several thousand. In addi-
tion, the -OH groups can be chemically modified both in nature or artificially. Several
modification such as methyl11 , fluorine11 , acetyl, etc. can be introduced to alter their
chemical properties.
The conformational parameters of the carbohydrates depend on their possibilities
to form intra-molecular hydrogen bonds between adjacent monomers, as shown in Fig-
ure 1.2. In addition, many possibilities of forming inter-molecular hydrogen bonds be-
tween polymer chains shown in Figure 1.2 leads to the formation of diverse aggregate
structures. Their flexibility, as well as the DP and monomer sequences of these molecules
also contribute to determine their solubility.
1.1.1 Polysaccharides: Cellulose and chitosan
Natural polysaccharides form the basis of many biomaterials. The two most abundant
polysaccharides are cellulose and chitin. Cellulose is a polysaccharide consisting of lin-
10
Figure 1.2: Intra- and inter- molecular hydrogen bonds between adjacent monomer and
polymer chains
ear chains of β(1-4) linked D-glucose units as shown in Figure 1.1a). It is an important
structural component of the primary cell wall of green plants, and many forms of algae.
Cellulose can reach a DP up to 1500012 . The many -OH groups present on the cellulose
polymers can make multiple inter-polymer hydrogen bonds between different polymer
chains as shown in Figure 1.2, forming crystalline microfibrils. The insolubility of cel-
lulose results from this dense network of inter-molecular hydrogen bonds, as well the
stiffness of the molecules produced by intra-molecular hydrogen bonds between the ring
oxygen (O5) and the -OH group at the C3 carbon. Altering the chemical constitution of
cellulose can alter the physical properties of the resulting network. Specific derivatives
can be designed that disrupt hydrogen bond networks and alter the physical properties
such as enhancing or decreasing the water solubility or change in ionic character13 .
Chitin is a polymer very similar to cellulose composed of β(1-4) linked N-acetyl
glucosamine as monomeric unit (as shown in Figure 1.1b). It has a similar structure
to cellulose with the hydroxyl group at the C2 carbon on each monomer replaced with
an acetyl amine group as shown in Figure 1.1b). The presence of the acetyl groups
increases the hydrogen bonding between the polymer chains further compared to that
of the hydroxyl groups. This result in the formation strong chitin polymer fibrils.
The N-acetylglucosamine monomers can be deacetylated chemically, to produce chi-
tosan. As a result, chitosan is composed of a random sequence of β-(1,4) linked glu-
11
cosamine(GlcN) and N-acetyl glucosamine(GlcNAc) monomers. GlcN comprises a pri-
mary amine group which can be protonated at pH value lower than its pKa. This results
in an electrostatic repulsion between the monomers, rendering chitosan chains to be sol-
uble. The primary amine group also provides a site for chemical modifications as well.
The presence of the N-acetyl group in GlcNAc, on the otherhand allow hydrophobic and
hydrogen bond interactions, leading to self-association between the monomers. As a
result, the self-assembly of chitosan primarily depends on the pH14, 15 and the degree of
acetylation(DA)16, 17 of the polymer. It had been suggested that with protonation more
than 75%, chitosan behaves as a poly-electrolyte, as chains have minimal association
due to high charge density. While chains with protonation below 50% or a DA of more
than 50% behave as hydrophobic polymer with only isolated positive charges. Neutral
chitosan polymers form different aggregate structures with varying hydrophobic mod-
ifications. These varying chitosan hydrogel has been useful in various application like
water purification18 , oil spill remediation17 , wound healing19 and drug delivery20, 21 .
1.2 Hydrogels
Hydrogels are highly porous networks of polymers that can contain as much as around
99 wt% water. The chemical and physical properties of the hydrogels can be modified
by chemically changing the constituents of the monomers, or the chemical nature and
amount of linker holding the polymers together. Due to their flexible micro structure,
they become a suitable candidate for drug delivery module20, 21, 22, 23 . Further, the load-
ing and release rates of drugs in and out of hydrogels are mostly governed by diffusion,
which depends on the molecular interactions between the polymer chains and the drugs,
as well as the network morphology. These features have been harnessed for tuning the
release kinetics and the scheduling of drugs, especially through chemical modification
of the polymer backbone. Due to its many desirable properties, chitosan has received
a lot of attention for potential applications as pharmaceutical hydrogel. With proper
modifications of chitosan can be used to alter the absorption, diffusion and the release
of small molecules from its hydrogels20, 21 .
Theoretical and computational models provide a method to explore the system pa-
rameters of such hydrogels on a large scale, while being comparatively cost effective.
Some of the methods that have been applied to predict structure and interactions of
polysaccharide assemblies and hydrogels are introduced in the next sections.
12
1.3 Molecular simulation
Molecular simulations provide essential tools to understand macromolecular structure
based on the molecular interactions. Simulations can be regarded as in-silico experiments
performed on the molecules. These computer based experiments can provide details
about the conformations and interactions of the molecules and the physical and chemical
properties of their aggregates.
At the most detailed level, quantum mechanics (QM) provides the most accurate
and fundamental description of matter. However, QM simulations can only be run for
at most few hundred atoms. Often systems need to be treated on larger length and
longer times scales than QM methods can achieve. On these scales, the fluctuations
of the electronic degrees of freedom play only a minor role. At even larger scales, the
same can be said for the motion of individual atoms. Many modeling approaches to
study a system on the required length and time scales have been developed, as shown
in Figure 1.3.
Figure 1.3: Different approaches to study the molecular systems with various resolution
depending on the properties of interest. In a bottom-up model development, coarser
resolution simulations are guided by detailed level studies like quantum mechanics. In
top-down approaches, macroscopic properties are used to guide finer-resolution simula-
tions like classical atomic or coarse grained molecular simulations.
13
The resolution required depends on the properties of interest and on the type of
phenomenon that needs to be analyzed. All-atom molecular dynamics(MD) and coarse
grained (CG)modeling are two such resolution as shown in Figure 1.4 for the solution of
GlcNH2monomers. These models will be discussed in the next sections.
Figure 1.4: Atomistic representation of (a) β-D-glucosamine solution and the corre-
sponding (b) coarse-grained representation.
1.3.1 Molecular dynamics
MD is a technique used to describe the positions and velocities of the molecules in the
system based on the Newton’s equations of motion as shown in Eq 1.2.
mi¨ri=fi(1.1)
fi=−∂U
∂ri
(1.2)
Here, miand riare the mass and position of atom i. To study the evolution of the
system with time, we need to calculate the forces fiacting on atom i. The forces are
usually described by a potential energy U(rN), where rN=(r1, r2, ...rN)represents the
complete set of atomic coordinates.
U(rN)is usually described by a set of interaction functions and corresponding pa-
rameters that is referred to as a force-field. Force-field functions and parameter sets are
derived from both experimental data and high-level quantum mechanical calculations.
Force-fields can be based on different parametrization principles and are specialized for
different applications and give varied results24, 25 .
14
Typically U(rN)is decomposed into bonded terms, relating atoms that are linked by
covalent bonds and non-bonded (also called “non-covalent”) terms, describing the long-
range electrostatic and van der Waals forces. A general form for the total energy in an
additive force field can be written as Utotal =Ubonded +Unon−bonded.
The non-bonded terms are computationally more costly, as they include many more
interactions per atom. The non-bonded interaction are most commonly given by the
Lennard-Jones potential (ULJ ) and the Coulomb potential (Ucoloumb). The components
of the non-bonded contributions are given by the following summations, Unon−bonded =
ULJ +Ucoulomb.
ULJ (rij)=4ϵij [(σij
rij )12
−(σij
rij )6](1.3)
Both, σij and ϵij in Eq 1.3 correspond to the equilibrium separation of two atoms i, j
and the depth of the energy minimum and rij =|ri−rj|is the distance between the
two atoms.
The coulomb interaction between charges or partial charges on the atoms is expressed
as:
Ucoulomb(rij) = qiqj
4πϵorij
(1.4)
where, the qi,qjare the charges and ϵois the permittivity of the free space.
The bonded potential typically comprises bond, angle, dihedral and improper dihedral
as shown in Figure 1.5. The components of the covalent contributions are given by the
summations: Ubonded =Ubond +Uangle +Udihedral +Uimproper.
Typically, a covalent bond between two atoms is modeled as harmonic potential and
expressed as,
Ub(rij) = 1
2kij(rij −req)2(1.5)
where, rij =|ri−rj|is the distance between the two atoms, req and kij are equilibrium
distance and spring constant.
A covalent angle is described by a harmonic angular potential of the form:
Ua(θijk) = 1
2kθ
ijk(θijk −θeq)2(1.6)
where θ=arccos rij rkj
rij rkj is the angle between atoms i, j and k. A simplified form can be
U(θijk) = 1
2kθ
ijk (cos(θijk)−cos(θeq)) (1.7)
15
Figure 1.5: Bonded interaction potentials include (a) bond, (b) angle, (c) dihedral, and
(d) improper dihedral
The dihedral angle ϕis formed by four atoms with indices i, j, k, and l.ϕis an angle
between the normal n and m to the two planes of i, j, k and j, k, l.
ϕ= arccos n m
|n || m |(1.8)
where n =rij ×rkj and m =rjk ×rlk. The dihedral angle potential is represented
as,
Ud(ϕijkl) = kϕ
ijkl (1 + cos(mϕijk −ϕo)) (1.9)
The another type of dihedral angle, i.e. improper dihedral is used to keep the groups
planar and prevent molecules from flipping over to their mirror images. This type of
dihedral is defined by a harmonic potential,
Uimproper =1
2kξ(ξ−ξo)2(1.10)
where ξis an improper dihedral angle and ξoits equilibrium value.
In this thesis, GLYCAM06TIP5P
OSMOr14
26, 27 force-field was used along with the TIP5P28 water
model throughout. As GLYCAM06TIP5P
OSMOr14 had shown good agreement with experimen-
tal free energy of hydration data for small saccharides29 .
16
1.4 Coarse-grained simulation
CG models helps to overcome the limitations of accessible length and time scales of
all-atom MD by grouping together atoms into CG interaction sites. This reduces the
number of degrees of freedom in the system and in addition typically creates a smoother
energy landscape, which leads to faster dynamics. As a consequence, such CG models
can be used to simulate larger systems for longer times. The gain in efficiency comes
at the cost of losing some of the chemical detail in the system. One of the challenges
lies in retaining enough information about the chemical details of the system in the CG
representation.
In creating predictive CG models, different strategies have been developed to find
effective interaction potentials between the CG sites. These potentials have enough
information of the underlying system to predict their large scale behavior reliably.There
are primarily two ways to transfer the information to obtain molecular interaction for
CG model as shown in Figure 1.3 :
1. bottom-up: fundamental physical principles at the more detailed scale are used
to parametrize a model at a CG scale
2. top-down: the behavior at larger scales is used to inform the interactions at more
detailed scales
In this thesis, a bottom-up approach is used to obtain the interaction potentials
for the CG model. The CG models consists of CG sites i.e. group of several atoms,
and the interaction potentials are derived form atomistic molecular simulations. The
two bottom-up coarse-graining methods used for that model, and had been successfully
applied in the development of CG models of polysacharides, the Multi-Scale Coarse
Graining (MS-CG)30 method and Boltzmann inversion,31 are described in more detail
in the next section.
1.4.1 Multi-scale coarse graining (Force Matching)
The idea of the force matching strategy is to reproduce the average force acting on CG
sites that are sampled in the all-atom system for the CG system. The force experienced
by the group of atoms in the all-atom system, averaged over the all-atom conformation
that correspond to the same CG sites is given as.
17
⟨fI⟩A=FI,CG for all CG sites I = 1, ..., NCG (1.11)
Here, ⟨fI⟩Ais the average atomic force and FI,CG is the CG force acting on site I. The
FI,CG can be used to generate UCG according to the relation given below,
FI,CG =−∂UCG(R)
∂Ri
(1.12)
where, UCG is the CG potential energy.
Both forces ⟨fI⟩Aand FI,CG here correspond to those on the CG-sites for the same
CG-sites configuration as shown in Figure 1.632 . The average forces on each CG sites
in the all-atom system are constructed using the mapping function and propagating the
individual atom forces to the CG sites. This procedure is applied to the entire reference
trajectory in the all-atom system to extract the interaction potential between the CG
sites.
Figure 1.6: Demonstrating force matching32 procedure by showing set of atomistic forces,
⟨fI⟩and its corresponding resultant CG force FIfor single water molecule.
The MS-CG method30 is used throughout to derive non-bonded interactions for
solute-solute, solute-solvent and solvent-solvent interactions separately. The separa-
tion ensured that solvent-solute interaction does not perturb the sensitive solute-solute
interaction. This was achieved by separating the atomistic trajectory into three sepa-
rate trajectories. Each individual trajectory has specific interactions i.e. solute-solute,
solute-solvent and solvent-solvent. A MS-CG method30 was applied on all the three
trajectory to obtain specific interaction potential.
18
1.4.2 Boltzmann inversion
Boltzmann inversion(BI)31 is based on matching the all-atom structure to obtain inter-
action potential for a given degree of freedom. It is based on the probability distribution
of a canonical independent degree of freedom which obeys the Boltzmann distribution
i.e.
P(r) = Z−1exp[−βU(r)] (1.13)
From the above equation, the probability distributionP(r)can be calculated from the
partition function i.e. Z =∫exp[−βU(r)]dr with β= 1/kBT(where kBrepresents the
Boltzmann constant and T is the temperature). The corresponding interaction potential
U(r)can be calculated from this probability distribution by
U(r) = −kBT ln(P(rN)) (1.14)
where P(rN)correspond to the probability of a CG configuration rNobtained from
the atomistic trajectory.
Bonded interactions such as bond, angle and dihedral potentials are calculated using
BI from canonical sampling of the system, which is given by
U(r)bond =−kBT ln(P(r)) (1.15)
U(θ)angle =−kBT ln(P(θ)) (1.16)
U(ϕ)dihedral =−kBT ln(P(ϕ)) (1.17)
where r, θand ϕare bond-length, angle and dihedral angle respectively. Pis the
probability distribution function for each degree of freedom.
In principle, these BI based interaction potential provide a good initial estimate for
bonded interactions and could be used directly to perform CG simulations provided
that the assumptions that the degree of freedoms are independent are approximately
fulfilled. If too many other interactions influence their conformations, the corresponding
interactions are typically overestimated and iterative optimizations steps are required.
19
1.4.3 Coarse-grained models for polysaccharides
A few CG models for polysaccahride systems have been developed until now. These
include the parametrization to reproduce bulk thermodynamic data for the popular
MARTINI model33 , sampling polymer conformations based on the conformational space
available to the glycosidic angles ϕand ψ15, 34, 35 , or deriving interaction potentials based
on the MS-CG method36, 27, 37 . In the latter model, a hybrid procedure was employed,
where force matching was used to obtain non-bonded interactions while bonded interac-
tions are obtained by BI. Non-bonded interactions for solute-solute, solute-solvent and
solvent-solvent interactions were derived separately. This method offers a promising ap-
proach to elucidate structure-property relations in saccharides systems, because it can
produce polysaccharide models that can be transferred to other concentrations as well
as to longer polymer chains27 and reproduces aggregation behavior and osmotic pressure
of the atomistic system37 .
To date, only few molecular modeling studies have addressed these for cellulose and
chitosan assemblies. However, the conformations of cellulose27 and single chitosan chains
at atomistic15 as well as in CG resolutions15 , have been simulated and the aggregation
of chitosan with different DA38 and different monomer sequences38 has also been studied
for charged polymers with a MARTINI33 like model.
1.5 Aim and overview of the thesis
The aim of this thesis is to better understand the self-assembly in these polysaccha-
rides systems based on the effects of various modification and physical conditions on
the aggregate structure. MD simulations have been employed to observe the forma-
tion and behavior of various aggregate structures and hydrogels formed from different
monosaccharide units. Simulations with all-atom resolution were employed to provide a
detailed picture of the local interactions and CG models for the different molecules were
developed and extensively validated to study these systems at much larger scale.
First, cellulose, chitosan and similar molecules with various chemical modifications
are simulated to understand how different perturbations of the intra- and inter-molecular
hydrogen bond network affect the aggregation of these molecules. Structural analysis of
the different morphologies generated by these polymer systems was performed. Then,
CG models for different chitosan hydrogels are introduced and tested and finally, the
potential applications of these models to optimize such systems for applications in fields
20
such as targeted drug delivery is demonstrated.
The rest of the thesis is structured as follows:
In chapter 2, the structure function relations for various cellulose derivatives and
chitosan derivatives are analyzed. First, all-atom results for single polymers and dense
solutions are shown. Then the development of CG models for all molecules and their
application to study the aggregate structures were described. Aggregate structure were
characterized by pore-size distribution, and contacts formed. Methylated, flourinated
and chitin analogues with various hydrophobic modification were prepared with different
patterns of substitution. Different aggregate structure were obtained for different types
of substitution as well as pattern of substitution.
Chapter 3 describes the development and validation of an efficient and transferable
CG model for chitosan. The model transferability across different concentrations, poly-
merization and degree of modification is explicitly tested. CG models for Doxorubicin
(DOX) and Gemicitabine (GEM)were developed and polymer-drug interaction potential
were also obtained.
In Chapter 4, chemically modified chitosan hydrogels were modeled for various con-
ditions and system parameters. The effect of water concentration, type, degree and
pattern of hydrophobic modification was investigated and their effects on the structural
characteristics of the hydrogel such as the pore-size distribution, the average number of
contacts and end-to-end distance of the polymers were analyzed.
Chapter 5, describes the migration of two model drugs molecules through the different
hydrated network structures of modified polysaccharide chains that were obtained in
chapter 4. Here DOX and GEM were chosen as model drugs. Initially, an all-atom
analysis was performed to study the dependence on the type of interaction between
both drugs and the polymers. CG simulations of their motion in the different chitosan
networks were performed with both the drugs separately or in combination.
Finally, Chapter 6 provides a conclusions and summary of the results.
21
Chapter 2
Tailor-made cellulose derivatives
2.1 Introduction
The mechanical and structural properties of polysaccharides depend upon the intra-
and inter- molecular hydrogen bonding of the polymers. Alterations in the monomeric
units provide a tool to modify the structural properties by disrupting these hydrogen
bonds. This approach can help us to understand what determines structure formation
and thus guide the development of novel, tailor made carbohydrate based materials.
It also provides insight to further understand the structure property relations of these
polymers. Modifications such as methylation, fluorination and acetylation and charged
glucosamine were used as shown in Figure 2.1. These modification have either hydropho-
bic or hydrophillic nature and were designed to selectively disrupt the hydrogen bond
between the hydroxyl groups and monomeric oxygens of the polysaccharides i.e. cellulose
and chitosan derivatives. Experiments with full control over the length and degree of
substitution were constructed using Automated Glycan Assembly method39, 40 and had
verified that the solubility and the gelation properties vary with these modifications as
compared to pure cellulose11 .
In this chapter, first all-atom MD simulations of cellulose and chitin derivatives are
described. Simulations were performed to understand the effect of different functional
groups on the molecular geometry and on the polymer-polymer association. In particu-
lar, we have also analyzed how the monomer sequence can lead to different conformations
of the oligosaccarides with the same monomer composition. An-other important factor
is the flexibility of the molecules which depends on the glycosidic dihedral angles shown
in Figure 2.2. In that regard, the effects of substitution on the torsion angles (ϕand ψ)
22
Figure 2.1: Example of modified cellulose structures (a) Methylation of alternat-
ing monomers (b) Alternative fluorination of alternating monomers increases the hy-
drophilicity. (c) Chitosan containing amine groups increases the positive charge in the
polymer (d) Chitin, acetylation of all monomers increases the hydrophobicity.
were calculated. The change in ψpopulation are directly related to the presence of a
hydrogen bond, and are therefore most affected by the modifications.
Then, to be able to study the self-assembly of these polymers on large time- and
length scales, bottom-up CG models were developed for all molecules. A CG force field
was developed using a hybrid approach based on MS-CG FM30 and BI31 to calculate
23
Figure 2.2: Definition of the dihedral angels ψ(C1,O4,C4,H4) and ϕusing the atoms
(H1,C1,O4,C4)
the non-bonded and bonded interactions, respectively. The CG model was validated by
comparing the radial distribution function (RDFs), end-to-end distances and radius of
gyration to data from the respective all-atom system. The CG force-field was trans-
ferred to longer polymers (DP=12) and lower polymer concentration and used to follow
aggregation in these systems. Different structural morphologies of the aggregates for the
different modification types, but also for different modification patterns were obtained.
2.2 Methods: MD details
The following atomic system were simulated with Gromacs 5.1.241 :
1. Single β-D 1-4 linked glucose with DP=6 in water.
2. Single β-D 1-6 linked glucose with DP=6 in water.
3. Single β-D 1-4 linked glucose with DP=6 in water with methylation and fluorinated
with different pattern i.e. alternate(A) and blocky(B) pattern.
4. Single N-acetyl-glucosamine. with DP=6 in water.
5. Single glucosamine polymer (charged and uncharged monomer) with DP=6 in
water.
24
6. 25 chains of system 1-5 with 2100 water molecules with DP=6.
Initial structures for the different cellulose and chitosan based molecules were con-
structed with tleap42 . The topologies were converted to gromacs format using the gly-
cam2gmx script43, 44 and subsequently solvated in GROMACS45 . The GLYCAM06TIP5P
OSMOr14
26, 27
force-field was used together with the TIP5P28 water model. Parameters for existing
modifications namely methyl, and acetylation were taken from the GLYCAM06 force-
field26 while parameter for fluorine were take from GAFF force field46 . Partial charges
for the modified monomers were calculated using the R.E.D scripts47 and following the
GLYCAM06 protocol26 . A cut-off of 1.4 nm was used for Lennard Jones and electro-
static interaction. Long range electrostatics were evaluated using Particle Mesh Ewald48
. Covalent bonds involving hydrogen atoms were constrained with the LINCS 49 algo-
rithm while, water molecules were kept rigid using SETTLE50 .
Energy minimization was performed following a standard protocol and a 50ns NPT
equilibration at 300 K and 1 bar, using the Nóse-Hoover thermostat51, 52 and Parrinello-
Rahman barostat53, 54 . Subsequently, a 400ns NVT equilibration run using the average
box size extracted from the NPT trajectory and the Nosé-Hoover thermostat51, 52 was
performed, followed by a 100 ns production MD run. A time-step of 2 fs was used and
energy and pressure dispersion corrections where appropriate have been applied.
To extract forces for the coarse-graining procedure, separate reruns of the MD tra-
jectories containing only solute-solute, solute-solvent or solvent-solvent interaction were
conducted27 . Long range electrostatics were calculated in the using the reaction-field
method55 using the same 1.4 nm cutoff as in the original simulations.
2.3 All-atom simulation results
2.3.1 Glucose and its derivatives
Two natural glucose oligomers having β1-4 and β1-6 linkage with DP =6 were simulated
to study the effect of the glycosidic linkage. As expected, the simulation snapshots shown
in Figure 2.3 represent linear and coiled structure for 1-4 and 1-6 linkage respectively.
The end-to-end distance also changed from 2.58±0.3 nm for 1-4 linkage to 1.28±0.2 nm
for (1-6) linkage portraying more coiled configuration with 1-6 linkage.
Next, derivatives of cellulose, i.e. β 1-4 linked glucose with different well defined
substitution patterns were modeled to study the effect of these moleculer structure and
25
Figure 2.3: Simulation snapshots in (a). β-D 1-4 linked Glucose and (b) for β-D 1-6
linked Glucose.
the properties of the aggregates formed. Modifications were introduced by replacing
the hydroxy group at C3 as shown in Figure 2.4 with methyl and fluorine. These
selectively disrupt the intra-molecular hydrogen bonds reducing the rigidity of cellulose
oilgosaccharides.
Figure 2.4: Chemical structure of (a). β-D- glucose methyl modified at C3 atom, (b).
β-D-glucose fluorine modified at C3 atom.
Methylation tends to disrupt the hydrogen bond between the O(5) and OH(3) desta-
bilizing the linear conformation of the cellulose. In addition, it is bulkier and represents
an additional steric hindrance. Two different patterns namely alternate methylation(AB)3
and di-block methylation (A3B3) were simulated to study the effect of the modification
26
pattern on the results.
The alternate(AB)3methyl substitution patterns shows very similar configuration as
pure cellulose. As in alternated case, there is slight decrease in the in the ψpopulation
as shown in Figure 2.5a), which result due to decreased tendency to form hydrogen
bond between methyl and monmeric oxygen. However, the same degree of methylation
with a block distribution leads to very different configurations with two diferent maxima
for different linkage as shown in Figure 2.5b). The distribution again show negative ψ
values for the links involving methylated monomers, however the sharp increase in the
negative ψvalues of the non-modified block indicates an increase of OH3-O5 hydrogen
bond formation. The corresponding simulation snapshots show a bent structure for the
Figure 2.5: Analysis of ψdistribution for (a) alternating methyl modified cellulose (b)
blockwise methyl modified cellulose, (c)alternating fluorine modified cellulose, and (d)
blockwise fluorine modified cellulose. The residues are numbered from the nonreducing
end to the reducing end.
block methyl pattern as compared to alternate pattern and pure cellulose as shown in
Figure 2.7. The end-to-end distance measured for the block methyl modification, reduces
from 2.7±0.2 nm for alternate methyl to 2.4 ±0.9 nm. The decrease in the average end-
27
to-end distance confirm overall structural changes in block methyl pattern. The time
trace of the end-to-end distance also show in Figure 2.6, highlights the more flexible
nature of the blockwise methyl modified cellulose as compared to the alternating methyl
modified cellulose.
Figure 2.6: Analysis of end-to-end distance (a) alternating methyl modified cellulose
(b) blockwise methyl modified cellulose, (c)alternating fluorine modified cellulose, and
(d) blockwise fluorine modified cellulose. The end-to-end distance was monitored over
100ns. The residues are numbered from the nonreducing end to the reducing end.
Similar to methylation, fluorination also prevents the hydrogen bond formation as the
OH group forming the hydrogen bond is no longer there, and affects the electron density
of the monomer. The replacement of -OH by electron withdrawing fluorine affects the
population of the ψangles as shown in Figure 2.5. However, the effect is small, the
fluorine modified cellulose molecules still have an overall linear conformation with both
patterns of modification i.e. alternate or block as shown in Figure 2.7. However, the
large distribution in average end-to-end distance shown in Figure 2.6c), present a very
28
flexible system for the alternated fluorine modified cellulose. The average end-to-end
distance between alternate 2.70±0.2 nm to block pattern 2.75±0.2 nm does not show
significant differences from each other. However, in alternate fluorination end-to-end
distance vs time plot show more flexible structure as compared to blockwise fluorination
pattern.
Figure 2.7: Simulation snapshots of system 1 and 3 with DP = 6 and 2000 water
molecules. The simulation snapshots show the carbon atom in gray, oxygen in red,
hydrogen in white, and fluorine in pink. The hydrophobic modification are encircled in
yellow and hydrophilic modification in blue for (a) glucose (b) alternated methyl modified
glucose (c) block methyl modified glucose (d) alternated fluorine modified glucose, and
(e) block fluorine modified glucose.
Thus overall the two chemically different substitutions did not shown drastic struc-
tural variations as compared to pure cellulose. The strongest change was found in the
configuration with block pattern methylation were seen.
29
2.3.2 Chitosan
Chitosan polymers are made up from three monomeric building blocks GlcNH2, GlcNH+
3
and GlcNAc as shown in Figure 2.10. Three polymers with single monomeric unit i.e. β
1-4 linked GlcNAc, GluNH2, and GluNH+
3were simulated and their snapshots are shown
in Figure 2.9.
Figure 2.8: Chemical structure of (a).β-D- glucosamine(NH2), (b) β-D-
glucosamine(NH+
3), and (c). β-D-Acetyl-glucosamine
Figure 2.9: Simulation snapshots of system 4 and 5 with DP = 6 and 2000 water
molecules. The simulation snapshots show the carbon atom in gray, oxygen in red,
hydrogen in white, and nitrogen in blue for (a) N-acetyl glucosamine (b) neutral glu-
cosamine (c) charged glucosamine.
As has been described before15 we see that both, the charge and the acetylation of the
monomers can have significant influence on the flexibility of the polymer.The maps in all
30
Figure 2.10: Analysis of end-to-end distance and conformational maps of ϕand
ψof chitosan (a) N-acetyl glucosamine (b) neutral glucosamine(NH2) (c) charged
glucosamine(NH+
3) obtained by MD simulations. The end-to-end distance was moni-
tored over 100ns.The residues are numbered from the nonreducing end to the reducing
end. The dihedral angles ϕand ψare shown on x- and y- axes, respectively.
cases shown in Figure 2.10 have a main minimum present at the same angle. However, the
charged monomer(NH+
3) had been shown a slight reduction of conformational flexibility
in comparison to neutral monomer(NH2)15 . The charge of the monomers also have
31
significant influence on the flexibility of the link. As, there is a existence of a second
minimum in case of full uncharged glucosamine as compared to charged one as shown in
Figure 2.10. The average polymer conformation is linear in all the cases in Figure 2.9 and
the end-to-end distance does not show significant difference as varying from 2.84±0.12
nm for Glc-NH2to 2.83±0.11 nm in Glc-NH+
3and 2.82±0.14 nm for GlcNAc.
2.4 CG model and mapping
Polymers behave differently in a crowded environment and different molecular properties
modulates the aggregation and structure of these polymer system. Along with that,
polymer systems usually have slow dynamics. An efficient approach was required to
study the self assembly of these system at longer length and time scale. A CG simulation
method was proposed as an appropriate tool to study their network structure formation.
In the CG model, each monosaccharide was mapped onto three coarse-grained in-
teraction sites. The modifications were represented by modified sites for the methyl,
amine and fluorine containing sites while the larger acetyl groups were modeled as an
additional site, as shown in Figure 2.11. Water molecules were represented as a single
CG site. The interaction potentials between these sites were generated using Boltz-
mann inversion31 for the bonded interactions and the MS-CG method30 for non-bonded
interactions, following the procedure described by Sauter et al27 .
2.4.1 CG simulations for short oligomers
CG simulations were performed with Gromacs 4.6.441 . Systems 1-5 containing single
polymer and 25 polymer chains with DP=6 and 2,500 water molecules of glucose, methyl
and fluorine modified glucose and charged as well as uncharged glucosamine were sim-
ulated. Two different patterns, alternate and blocks of three were simulated for methyl
and fluorine modified cellulose. The initial CG structures for all the short oligomers
were obtained from their corresponding atomistic representation. All CG simulation
were conducted in the NVT ensemble using the Leap-Frog integrator56 with the Nóse-
Hoover thermostat51, 52 . A 1.4 nm cut-off was used for CG simulations throughout.
A 10ns production run was done and results were used for validation of the model by
comparing the RDFs, end-to end distances and radius of gyration obtained from the CG
and atomistic system.
32
Figure 2.11: All atom and coarse-grained representation of (a). β-D- glucose, (b). methyl
modified glucose (c). fluorine modified glucose (d) N-acetyl glucosamine
2.5 Model validation
2.5.1 Cellulose
To validate the performances of CG force field, RDFs between all interaction sites of
the CG model were compared and showed excellent agreement with their equivalent
obtained from all-atom simulations. The RDFs of the atomistic and CG trajectories of
β-D glucose with DP=6 and 2100 water molecules are shown in Figure 2.12 for the CG
sites A and B with themselves and with water.
The RDFs show good overall agreement of the short range structural features, and
the A-WAT and B-WAT RDFs, which were previously found to be the most sensitive to
perturbations27 , show that aggregation behavior is captured correctly in the CG-FFs.
Overall, local structure and features are well represented in all CG RDFs. A CG model
for β1-6 linked glucose was also developed and showed similar agreement of RDFs as β
1-4 link glucose.
33
Figure 2.12: Comparsion of AA and CG RDFs for interaction involved with CG site
type A and B with itself and water. (a) A-A (b) A-WAT, (c) B-B, and (d). B-WAT
2.5.2 Methyl modified cellulose
Next, CG-FFs were developed for methyl modified cellulose to study the network struc-
ture formed by these molecules, as compared to pure cellulose. As for all-atom system,
two patterns, block and alternate were analyzed as shown in Figure 2.13. The CG in-
teractions were developed for both patterns, alternate and block, separately. The model
performance was analyzed by comparing RDFs between M(modified) beads and between
M beads and water. Due to the hydrophobic nature of the methyl modification, it was
expected to plays an important role in polymer-polymer association.
Comparing the RDFs shown in Figure 2.13, differences in the short range structure
between the two modification patterns become apparent. The short range peak at 5Å and
10Å contain large contributions from the first and second bonded neighbors, respectively.
In the alternate pattern, all modifications with in the same molecules are 2nd neighbor,
so the second neighbor peak at 10Å is very pronounced. The peak at 5Å in this case
comes entirely from the non-bonded assembly. It is noticable, that in the CG model this
34
Figure 2.13: Comparison of AA and CG RDFs and their corresponding polymer snapshot
for interaction involved for (a) M-M (b) M-WAT (c) snapshot, for alternate methyl
modification and (d) M-M (e) M-WAT and (f) snapshot for block methyl modification.
peak is reduced compared to the atomistic case, whereas the peak at 10Å is very well
captured.
Figure 2.13d) has a less pronounced second peak compared to Figure 2.13b) because
in the block pattern as there is only one second neighbor pair. Instead, the 5Å peak
has drastically increased as a result of the bonded neighbors and very good agreement
between the CG and all atom representation is found. The water RDFs show similar
behavior irrespective of modification pattern and show overall good resemblance with
the atomistic counterpart.
Similarly, the fluorine modified cellulose with both pattern of modification was de-
veloped. Overall, it showed good agreement between the RDFs for all-atom and coarse-
grained simulation respectively.
35
2.5.3 Amine (-NH2and -NH+
3) groups
The systems with amine groups (-NH2and -NH+
3) were evaluted as shown in Figure 2.14.
Similar to the methylated system, the RDFs clearly demonstrate that the developed force
field was able to capture the overall aggregation behavior and the local structure of the
atomistic system well.
Figure 2.14: Comparison of AA and CG RDFs of neutral and charged glucosamine and
their corresponding polymer snapshot for interactions (a) A-A (b) A-WAT (c) snapshot
for uncharged system (d) A-A (e) A-WAT and (f) snapshot for charged system.
The RDFs shown in Figure 2.14 (a & b) correspond to uncharged glucosamine, while
Figure 2.14 (d & e) correspond to charged glucosamine. Comparing the RDFs, the
minimum bead distance between the charged A sites has increased in Figure 2.14d) as
compared to that of the neutral molecules in Figure 2.14a). This shift is properly cap-
tured by CG FFs. In addition, the charged monomer leads to a more ordered structure
of water beads around the charged beads. Figure 2.14e) clearly shows water shell for-
36
mation around the charged bead. Overall, the location of water shells in the CG have
good agreement with the atomistic structure, but are less pronounced and are smoothed
out after the third shell.
2.5.4 Single polymer conformation
The end-to-end distances for single polymers were calculated and compared with their
corresponding atomistic end-to-end distances. They provide useful information about
the overall conformations of the molecules and about the effect of type and pattern of
modification on the individual molecules. The end-to-end distances are summarizes in
Table 2.1. In most cases, atomistic and CG distances show good resemblance with each
other. However, in the case of β1-6 linked glucose, a significant difference between the
end-to-end distance for atomistic and CG simulation appear. The 1-6 linked glucose tend
to form a coiled structure in atomistic simulation, which seems to linearize the structure
in the CG simulation. A possible explanation is that, the BI bonded interaction are too
over-determined, making 1-6 linkage more stiffer. An solution can be to use iterative BI
to make molecule more flexible and get a better glycosidic link sampling. Despite this
observation, also in the CG model the 1-6 linked glucose has the shortest end-to-end
distance and very large fluctuations, indicating that the high flexibility of the 1-6 link is
at least partially captured by the CG model.
Atomistic Coarse-Grained(CG)
Modification End-to-End
distance(nm)
Radius of Gy-
ration(nm)
End-to-End
distance (nm)
Radius of Gy-
ration (nm)
Glu(1-4) 2.58 ±0.32 0.90 ±0.00 2.86 ±0.12 0.93 ±0.02
Glu(1-6) 1.28 ±0.26 0.60 ±0.04 1.96 ±0.62 0.78 ±0.14
Glu-OMe(A) 2.73 ±0.22 0.93 ±0.04 2.79 ±0.19 0.92 ±0.03
Glu-OMe(B) 2.47 ±0.88 0.93 ±0.06 2.70 ±0.22 0.89 ±0.03
Glu-F(A) 2.70 ±0.24 0.91 ±0.04 2.73 ±0.21 0.90 ±0.03
Glu-F(B) 2.75 ±0.19 0.93 ±0.02 2.75 ±0.19 0.91 ±0.03
Glu-NAc 2.82 ±0.14 0.97 ±0.01 2.65 ±0.14 0.95 ±0.01
Glu-NH22.80 ±0.12 0.94 ±0.01 2.81 ±0.15 0.92 ±0.02
Glu-NH+
32.83 ±0.11 0.94 ±0.01 2.87 ±0.08 0.93 ±0.01
Table 2.1: Comparison between atomistic and CG end-to-end distance and radius of
gyration for cellulose and chitosan and their derivatives. Here A and B correspond to
alternating and blockwise modification respectively
37
2.6 Aggregation and aggregate properties
After validation, the CG models were applied to model the self assembly of polymer
systems at greater length and longer time scales. Network structures of these polymers
with different type and pattern of modifications were generated in the CG simulations
and subsequently analyzed in terms of their pore sizes, the end-to-end distances of the
polymers and the average number of contacts formed between the polymers as well as
between the modifications. The distribution of pore sizes, was calculated by the method
described in Ref57 . This procedure finds the largest sphere that can be constructed to
contain randomly selected points in the network. This is achieved with a constrained
nonlinear optimization of the center of the sphere using the SOLVOPT routine58 .
CG systems containing 100 polymer chains with DP = 12, and 40 000 water beads
(400 water molecules per chains) with different modification, were generated by fol-
lowing the same steps as described in the previous section for the hexamers. Initial
CG structures were obtained from a 10 ns NPT equilibration simulation run using the
Nosé-Hoover thermostat51, 52 and Parrinello-Rahman barostat53, 54 with all-atom resolu-
tion. Then, the atomistic structures were mapped to their CG representations using the
VOTCA package59, 60 . All the CG simulation were simulated for 100ns in NVT at the
optimal volume obtained in the NPT run with the Nosé-Hoover thermostat51, 52 .
The following system were simulated with Gromacs 4.6.441 .
1. 100 cellulose chains with DP =12 and 40,000 water beads.
2. 100 chitin chains with DP =12 and 40,000 water beads.
3. 100 cellulose chains with DP =12 having methyl modification and 40,000 water
beads with alternate modifications pattern.
4. 100 cellulose chains with DP =12 having methyl modification and 40,000 water
beads with modifications in blocks of three.
5. 100 cellulose chains with DP =12 having fluorine modification and 40,000 water
beads with alternate modifications pattern.
6. 100 cellulose chains with DP =12 having fluorine modification and 40,000 water
beads with modifications in block of three.
7. 100 chitosan chains with neutral glucosamine with DP =12 and 40,000 water beads.
38
8. 100 chitosan chains with charged glucosamine with DP =12 and 40,000 water
beads.
2.6.1 Cellulose and chitin structure
First we characterized the aggregation behavior of the natural polysaccarides, cellulose
and chitin. Both polymers self-assemble into fibrils from different initial conditions.
Simulation snapshots of the cellulose and chitin cluster are shown in Figure 2.15. All
polymers have straightened out and the average end-to-end distance for cellulose is 6.03
±0.11 nm while for chitin it is 5.74 ±0.07 nm. On further analysis, approx. 50% of
the polymer are aligned anti-parallel to each other while the rest are aligned parallelly.
The self-assembled structure of both, cellulose and chitin show some twisting of the
fibrils. The crystal structure of α-chitin and cellulose have been resolved61, 62 and can
be compared to the self-assembled fibrils in our simulation model. Alignment is not
as optimal in crystal structure, but structure with both alignments had been found
experimentally for different conditions and thus neither orientation is not completely
unfavorable. Regarding twisting of the fibrils, it could be brought about by the assembly
kinetics e.g. to accomodate initial contacts formed with different molecules. On the other
hand, twist is also observed in many atomic simulations of cellulose fibrils29, 63 and was
found to depend sensitively on the interplay of different force field contributions64 .
Figure 2.15: Simulation snapshots of (a) cellulose and (b) chitin networks with 100
polymers of DP=12 and 32 water molecules/monomer.The simulation snapshots show
the polymer backbone (A,B,C beads) in red, modifications (M beads) in yellow and
water molecules as blue dots
39
A similar model for chitin was used for longer polymer with DP=50, as described in
chapter 4 and showed to form fibrils as well. The fibril structures formed by the cellulose
and chitin oligomers show that the CG model closely captures the aggregation behavior
that is expected for these molecules.
2.6.2 Methylated and fluorinated cellulose hydrogels
Simulation snapshots of the systems with 100 polymers, DP =12 and 40,000 water beads
of cellulose with methyl and fluorine modification are presented in Figure 2.16. These
polymer have either an alternate or a block (3 modifications) pattern of methyl and
fluorine modifications. Strikingly, the simulation snapshots and pore-size distributions
differ dramatically for the different modification patterns demonstrating that pattern can
have a dramatic effect on the aggregation behavior of the molecules. In case of alternated
methyl modification, polymer aggregate to aligned fibrils was observed similar to that
found for unmodified cellulose while with block modification, the polymer was soluble
and distributed evenly through the box.
This is reflected in the pore size distribution in Figure 2.16c), which shows a large
shift in the pores sizes, with large pores corresponding to the separate solvent phase
for the alternate pattern. The pore size decreases from 4.0 nm in alternate pattern to
1.5 nm in block pattern. The end-to end distance had increased from 5.1 ±0.0 nm in
blockwise to 5.7±0.1 nm in alternate patterns, showing the effect of polymer aggregation
which leads to alignment and thus a more linear configuration for alternate pattern. The
average number of contacts between the M beads was also calculated. As obvious from
the snapshots the average number of contacts decreased drastically from 638±86 for the
alternate pattern to 270±23 in the blockwise pattern. The average number of contact
increased with alternate pattern due to the aggregation of the polymer. The solubility
in block pattern can be due to the steric hindrance of the adjacent methylated monomer,
causing non-linearity in the polymer. The non-linearity was found in the snapshots of
both atomistic as well as coarse-grained snapshot of single polymer shown in Figure 2.7
c) and Figure 2.17c) respectively. Overall, pattern of modification plays a bigger role in
polymer aggregation.
Surprisingly for fluorine modifications, the opposite effect was observed. The polymer
were soluble for the alternate modification pattern but formed highly aligned aggregates
with the block pattern. The pore size distributions shown in the Figure 2.16 illustrate
an increase in the pore size from alternate to blockwise fluorinated pattern. Similarly,
40
Figure 2.16: Simulation snapshots and pore-size distributions of cellulose networks with
100 polymers of DP=12 and 32 water molecules/monomer. The simulation snapshots
show the cellulose backbone (A,B,C beads) in red, methyl modifications (M bead) in
yellow, fluorine (F bead) modifications in blue, and water molecules as blue dots for
(a) alternately methylated cellulose (b) blockwise methylated cellulose (d) alternately
fluorinated cellulose, and (d) blockwise fluorinated cellulose.
the end to end distance changes from 4.8±0.01 nm in the alternate to 5.6±0.1 nm in the
blockwise pattern. The polymer has become more linear with block modification, due to
the aggregation of the polymers as shown in Figure 2.16e). Again, the average number
of contacts between the F beads, had increased from the alternate to block patterns i.e.
167±23 to 1000±93 portraying the strong polymer aggregation in block pattern. It is
noticable, that the increase of contacts in the aggregates is much higher than the one
observed for methyl modification, by a factor of 5.9 as compared to 2.4.
The data for end-to-end distances of single polymer with different modifications and
modification pattern are shown in Table 2.2 compared to the dense solution.
41
End-to-End Radius of End-to-End
Modification distance(nm) Gyration(nm) distance(nm)
(Single chain) (Single chain) (100 chains)
Glu(1-4) 5.56 ±0.50 1.77 ±0.08 6.0 ±0.1
Glu-OMe(A) 5.13 ±0.72 1.70 ±0.11 5.7 ±0.1
Glu-OMe(B) 5.08 ±0.74 1.67 ±0.12 5.1 ±0.0
Glu-F(A) 4.63 ±0.82 1.58 ±0.14 4.8 ±0.0
Glu-F(B) 5.24 ±0.64 1.70 ±0.10 5.6 ±0.1
Glu-NAc 5.48 ±0.38 1.81 ±0.59 5.7 ±0.1
Glu-NH25.46 ±0.51 1.74 ±0.86 5.4 ±0.5
Glu-NH+
35.78 ±0.25 1.80 ±0.04 5.7 ±0.2
Table 2.2: Polymer end-to-end distance for single polymers in solution. The errors
represent one standard deviation. Here A and B correspond to alternating and blockwise
modification respectively
The single cellulose and chitin polymer resemble closely with the network end-to-end
distance. For the methyl-modified polymers, no great differences were found between
the different modification patterns. Both have an end-to-end distances that are slightly
lower than pure cellulose, and similar standard deviation which mark them as more
flexible than pure cellulose. Although comparison of the hexamers with all-atom re-
sults has shown, that it is possible that the CG model may not correctly capture the
greater flexibility of the blockwise-methylated cellulose, the CG model predicts strong
differences for the aggregation of these two molecules, so there must be some difference
captured by the model. A possible explanation may lie in the different shape of the
two patterns, which can be seen in the snapshots in Figure 2.17c). While the molecule
with an alternate pattern curves smoothly, the blockwise pattern leads to sharp kinks
in the molecular structure, which makes it less suitable for stacking. Differences in
flexibility were found for the different patterns of fluorination. The fluorine modified
polymer with the alternate pattern showed a reduced end-to-end distance and larger
as compared to the blockwise modification pattern. Thus the blockwise fluorine mod-
ification causes the polymer to be stiffer which facilities to aggregation compared to
the more flexible alternate fluorination. In addition, it is possible that the blockwise
pattern leads to enhanced interaction between the modification, as suggested by the
strong increase in average contacts between the F beads in the aggregates. Oligosac-
charides with blockwise and alternate patterns of modifications have been produced
42
experimentally using the Automated Glycan Assembly method39, 40 and have confirmed
the opposing trends observed for the fluorine and methyl modification. While molecules
with blockwise methylation and alternating fluorination appeared completely soluble11
, the XRD spectra of molecules with alternate methylation and blockwise fluorination
showed characteristics resembling those observed for pure cellulose.
Figure 2.17: Simulation snapshots of single cellulose polymers with DP=10 and 40000
water molecules/monomer. The snapshots show the cellulose backbone (A,B,C beads) in
red, methyl and acetyl modifications (M bead) in yellow, fluorine modification F in blue
and water molecules as blue dots for a) pure cellulose (b) 50% methylated cellulose with
alternating pattern (c) 50% methylated cellulose with block pattern (d) 50% fluorinated
cellulose with alternating pattern (e) 50% fluorinated cellulose with block pattern (f)
chitin.
2.6.3 Chitosan hydrogel
Chitosan self-assembly is governed by both charge and DA. The effect of DA will be
discussed thoroughly in chapter 3 and 4. Chapter 3 shows the model development of the
chitosan. Chapter 4 shows the self-assembly for longer polysaccharides with DP =50
and also the effect of degree and pattern of acetylation as well as other hydrophobic
modification namely butyl and heptyl on hydrogel structure. Here, the effect of charges
in the full deacetylated chitosan is investigated.
43
Simulation snapshots of glucosamine polymer with neutral and charged monomers of
100 chains with DP =12 and 40,000 water beads are shown in Figure 2.18 a) & c). On
first sight both solution appear similar, with the polymer distributed evenly through the
simulation box. The difference in the network structure can be observed from the pore
size distributions, also shown in Figure 2.18 c). The charged glucosamine network have
smaller pore size diameter of 1.2 nm, while uncharged glucosamine have 1.5 nm pore-
size diameter. This is caused by the electrostatics repulsion between the charged sites,
causing the chains to distribute as evenly as possible through the box to maximize the
distance between the charged sites, whereas the slight association between the neutral
chains frees up space for larger pores in the solvent phase. This is also reflected in the
minimum distance between polymer, changing from 0.29 nm to 0.24 nm in charged to
uncharged networks. The charged network is soluble, as the box had fix dimensions
this result in the maximum spacing. Thus the model characterize the swelling and de-
swelling of the polymer chains well, as far as the fixed box size allows. To fully capture
the swelling of the charged chains the system would have to be coupled to a water
reservoir, that allow the charged chains to dispense.
Figure 2.18: Simulation snapshots and pore-size distributions of chitosan networks with
100 polymers of DP=12 and 32 water molecules/monomer. The simulation snapshots
show the chitosan backbone (A,B,C beads) in red, charged glucosamine A in blue, and
water molecules as blue dots for (a) Glucosamine(NH2) (b) Glucosamine(NH+
3).
44
2.7 Conclusion
Cellulose can be chemically modified to form various derivatives which poses the desired
properties for specific application. In this chapter, we have shown that a combination
of all-atom simulations and a systematic coarse grained model based on the all-atom
interactions, can be used to efficiently and precisely predict the effect of the modifications
on the aggregates. Chemical alterations such as substitution of one hydroxyl group with
methyl groups or fluorine atoms were modeled to study the self-assembly of these polymer
systems by different functional groups. The CG model revealed insights into properties
of the network such as the pore size, end-to-end distance and the minimum distances
between the polymers.
45
Chapter 3
A multiscale model for
hydrophobically modified chitosan
3.1 Introduction
Chitosan is a polymer of major interest to researchers and clinicians for developing ther-
apeutic hydrogels. It is derived from naturally abundant chitin and is a bio-compatible,
nontoxic polymer that is degradable by human digestive enzymes20 . In addition, the
presence of the primary amine groups on the glucosamine monomers provide a site for
the chemical modification, which has been extensively exploited to tailor the kinetics
of drug release65, 66 . As a result, a variety of chitosan-based formulations have been
developed for oral, ophthalmic, and transdermal applications,20, 67, 68, 69 several of which
have received FDA approval, demonstrating clear feasibility of these materials toward
clinical translation.
In this thesis acetyl, butanoyl, and heptanoyl moieties as shown in Figure 3.1 were
chosen for chitosan modification as they represent similar but increasingly hydrophobic
modifications and therefore allow us to study systematically their effect on the properties
of modified chitosan polymers and network.
To understand how these above mentioned chemical modifications govern the mor-
phology of the hydrogel. It is necessary to model chitosan hydrogels across a large
length and time scale, where a high number of long chitosan polymer chains, and many
water molecules and their dynamic interplay can be simulated, for a sufficient time to
render a physically accurate representation of these systems. This is not possible with
all-atom simulations. We therefore resolved to adopt a CG modeling procedure that was
46
Figure 3.1: Chemical structure of (a). N-acetyl glucosamine, (b) N-butyl glucosamine,
and (c) N-heptyl glucosamine.
informed by the results of the all-atom simulations and that has been shown to capture
the behavior of the polysaccharide-based hydrogel70 .
In this chapter, the development and validation of a bottom-up CG model for the
modified chitosan molecules is described. An atomistic simulations were performed to
obtain the CG potentials for polymer-polymer and polymer-solvent interactions. As
shown in chapter 2, to obtain a predictive CG model, the MS-CG30 procedure was used
to obtain the non-bonded interaction while BI31 was used for bonded interactions. The
model is validated by comparing RDFs and end-to-end-distances and radius of gyration
against all atom data.
3.2 Methods
3.2.1 All-atom simulations
For the development of the CG interactions, all-atom MD simulations of the following
systems were performed using Gromacs 5.1.241 .
1. 10 chitosan chains with degree of polymerization (DP) = 16 and 0% acetylation,
with 200 water molecules per oligosaccharide, i.e. 2000 water molecules in total.
2. 10 acetylated chitosan chains with DP = 16 and 200 water molecules per oligosac-
charide for degrees of acetylation of 16 %, 24 %, 32%, and 50 %.
3. 10 chitosan chains with DP = 16 and 200 water molecules per oligosaccharide with
degrees of butylation 16 %, 24%, 32%, and 40 %.
4. 10 chitosan chains with DP = 16 and 200 water molecules per oligosaccharides
with degrees of heptylation 8 %, 16%, 20%, and 24 %.
47
5. 50 DOX molecules with 10 water molecules per DOX.
6. 50 GEM molecules with 10 water molecules per GEM.
7. Systems 1 to 4 with the addition of 10 DOX or 10 GEM molecules.
The acetyl, butyl, or heptyl modifications in systems (2-4) were uniformly distributed
along the chains.
In addition, single polymers with DP=16 in a water box were simulated to compare
the end-to-end distance and radius of gyration between atomistic and CG simulations.
The computational procedure used here was similar that described in Chapter 2. Param-
eter for the butanoyl and heptnoyl modifications and the drug molecules were obtained
from the gaff force field46 . Partial charges for the modified glucosamine monomer, DOX
and GEM were obtained by the same way mentioned in Chapter 2. Other parameter
for the carbohydrates were taken from the modified GLYCAM06TIP5P
OSMOr14
26, 27 force-field
with the TIP5P28 as water model.
3.3 CG Model
As similar to previous chapter 2, glucosamine monomers were mapped to three CG sites
- A, B and C - as shown in Figure 3.2 using center-of-mass (COM) mapping. Acetyl
or butyl or heptyl modifications were mapped to one, two or three CG M sites, as also
shown in Figure 3.2. DOX was mapped to 11 CG sites, using 9 distinct bead types, to
produce a structure that can capture the planar geometry of the DOX molecule. Four
distinct beads were used to model GEM. Water molecules are mapped to one CG site
located at its center of geometry (COG).
Potentials for bond, angle and dihedral interactions, as well as specific non-bonded
1-3 and 1-4 interactions were obtained from Boltzmann inversion31 , using the VOTCA
package59, 60 . The non-bonded interaction were obtained using the MS-CG Method30 us-
ing the rerun trajectories with separate solute-solute, solute-solvent and solvent-solvent
interactions. All the bonded and intra-molecular interaction were excluded during the
MS-CG procedure.
48
Figure 3.2: All-atom and coarse-grained representations of (a) DOX, (b) GEM, (c)
unmodified glucosamine monomer, (d) acetyl-glucosamine, (e) butanoyl-glucosamine,
and (f) heptanoyl-glucosamine.
49
3.3.1 CG Simulations
CG simulations of systems corresponding to the all-atom systems 1-7 were performed
with Gromacs 4.6.441 .
Initial structures were obtained from the corresponding atomistic representations.
CG simulations for each system were performed in the NVT ensemble using the Leap-
Frog integrator56 with the Nosé-Hoover Thermostat51, 52 and a time step of 1 fs. A
cut-off 1.4 nm was used for all the CG simulations. Simulations were run for 10ns to
obtain equilibrated RDFs, angle distributions, average end-to-end distances and radius
of gyration to compare between all-atom and CG simulation.
3.4 Model validation
To evaluate the CG model’s ability to reproduce the local molecular structure of the
solution, as well as the overall tendency of chitosan molecules to aggregate, we compare
the radial distribution functions (RDFs) obtained from the CG simulations of chitosan
chains with DP = 16 to those obtained in the atomistic simulations of the same system.
The ensemble of RDFs between all pairs of CG sites was calculated and offers information
both on the short range molecular structure, reflected in the position and magnitude of
the peaks at short distances, and the overall aggregation trends of the chitosan chains,
visible in the long range behavior of the curve. were compared to obtain a measure of
the solution structure.. To characterize the conformation of single chitosan chains, angle
distributions and end-to-end distances were also calculated.
First, the RDFs obtained from CG simulations of chitosan chains with low degree of
modification, χAc = 16% (Figure 3.3), χBut = 16% (Figure 3.4) , and χHep = 8%
(Figure 3.5) aligned well with the corresponding RDFs obtained from atomistic sim-
ulations . Of particular notice is the agreement between CG and all-atom data for the
RDFs between water beads and the various CG sites of chitosan chains Figure 3.6, since
carbohydrate-water interactions were previously found to be highly sensitive to long
range perturbations and sampling issues27 . The RDFs for all other pairs of CG sites
obtained for modified chitosan chains with low χfrom CG and atomistic simulations
also showed good resemblance with each other.
For the bonded interactions, the most flexible degrees of freedom of polysaccharide
systems are those of the glycosidic bonds. In the present CG model, the conformations
of the glycosidic dihedral angles, ϕand ψare reflected by the angles Bi−Ai−Ci+1 and
50
Figure 3.3: Comparison of CG and atomistic RDFs of the distances between CG (A, B,
C, M for χAc=16% . Note: A, B, and C beads map the GlcN monomers, whereas M
map the modification(acetyl) group.
51
Figure 3.4: Comparison of CG and atomistic RDFs of the distances between CG (A, B,
C, M for χBut=16% . Note: A, B, and C beads map the GlcN monomers, whereas MA
and MB map the modification groups
52
Figure 3.5: Comparison of CG and atomistic RDFs of the distances between CG (A, B,
C, M for χHep=8% . Note: A, B, and C beads map the GlcN monomers, whereas MA,
MB, and MC map the modification groups
53
Figure 3.6: RDFs for interaction involved with CG site type WAT with other beads for
(a),(b),(c) and (d) 16 acetylation and (e), (f) for 16 % butylation.
54
Ai−Ci+1 −Ai+1 as well as the dihedral angle Bi−Ai−Ci+1 −Bi+1. The probability
distributions sampled for these angles in the all atom as well as simulation snapshots
corresponding to the three minima in the ϕ−ψfree energy landscape are shown in
Figure 3.7. Both angle distributions show two distinct maxima, whereas in the atomistic
system, the dihedral distribution has three maxima at -140◦, 30◦and 140◦, corresponding
to the three conformations shown in Figure 3.7.
Figure 3.7: (a-c) Angle distributions in the atomistic and CG models; (d-f) molecular
conformations corresponding to the three free energy minima of the ϕ−ψdihedral angles.
All-atom models are drawn as grey sticks, CG molecules as red (ABC) and yellow (M)
beads.
55
A comparison of the atomistic and CG distributions shows that the CG model gives
a good representation of the dominant conformation, but that the second energy min-
imum is under-represented in the CG model. This is a consequence of applying the
inverse Boltzmann method without further iteration, which does not account for the
effects of neighboring bonds and may therefore sometimes over-represents the stiffness
of interaction potentials. However, the population of the second minimum is relatively
small in the atomistic model, and will therefore only have a minor effect on the overall
polymer conformation, as illustrated for example by the comparison of the end-to-end
distances of the of the atomistic and CG oligosaccharides shown in Table 3.1.
Atomistic Coarse-Grained (CG)
Modification End-to-End
distance(nm)
Radius of Gy-
ration(nm)
End-to-End
distance (nm)
Radius of Gy-
ration (nm)
16% Acetylation 7.14 ±0.35 2.33 ±0.05 6.80 ±0.87 2.19 ±0.14
16% Butylation 7.11 ±0.44 2.20 ±0.08 5.95 ±1.18 2.02 ±0.19
8% Heptylation 5.83 ±0.68 2.06 ±0.09 5.94 ±1.38 2.00 ±0.21
Table 3.1: Table showing comparison between atomistic and CG end-to-end distance
and radius of gyration for (a). 16% acetylation ,(b). 16% Butylation and (c). 16%
Heptylation
For higher degrees of modification however, it becomes apparent, that the CG RDFs
for the hydrophobic modifications become progressively worse. The RDFs of the M-beads
for 32% acetylated chitosan, shown in Figure 3.8a), exhibit a strongly exaggerated peak
at short distances. Similarly, the RDFs for 32% butyl modification shown in Figure 3.8b)
and c) for MA-MA and MB-MB interactions differ between the atomistic and CG sys-
tems. Whereas the MB-MB CG RDFs in Figure 3.8c) capture the overall atomistic
behavior, although with a reduced magnitude of the short distance peaks, the MA-MA
RDFs differ substantially at short distances, and show un-physically close contacts. The
all-atom RDF on the other hand contains a number of irregular peaks at all distances,
which indicates the formation of clusters and suggests that there may be problems for
accurately sampling the distribution of atomistic forces as a result of the strong inter-
actions between the hydrophobic modifications. Because previously potentials obtained
from a similar coarse-graining procedure were found to be transferable to different con-
centrations,70 we tested the use of M-bead interaction potentials obtained at a lower
degree of acetylation χ= 16%, where they performed well. The results are shown in
56
Figure3.8 a) together with those from the native CG model.
Figure 3.8: RDFs of the distances between modification beads: (a) M-M beads in acetyl-
chitosan with χAc = 32%; (b) MA-MA and (c) MB-MB beads in butanoyl-chitosan with
χBut = 32%; and (d) MA-MA, (e) MB-MB, and (f) MC-MC beads in heptanoyl-chitosan
with χHep = 16%. The RDFs obtained from the atomistic, native CG, and CG with
transferred potential models are in black, red, and blue, respectively.
Comparison shows, that the transferred potentials significantly improve the over-
aggregation of the M-beads for acetylated chitosan, and the CG RDF for M-M inter-
57
actions now closely resembles the atomistic one. Applying the same approach of using
interaction potential from 16% butylated for 32% butylated, the CG RDF for MA-
MA interactions also shows significantly reduced clustering and now resembles that of
the acetylated chitosan, whereas the MB-MB interactions (Figure 3.8c) are unchanged
compared to the CG model obtained explicitly for χ= 32%. All other RDFs remain
unchanged for the transferred modification interactions. Analogous results were found
for the heptyl modified chitosan chains. Because the heptyl modifications are more hy-
drophobic, effects of clustering was already seen for χ=16% modification, so that the CG
interaction potential were obtained for χ=8% heptylation and transferred to χ=16% and
χ=24% heptylation. The comparison of RDFs for the MA-MA interactions in heptanoyl-
chitosan chains again shows a better representation of the short range features in the
atomistic RDF by the transferred interaction potentials as shown in Figure 3.8d). MB-
MB and MC-MC interactions on the other hand, remain unchanged Figure 3.8e) & f).
Similarly, the RDFs between all other CG sites remain unchanged when the transferred
interaction potentials are used.
Collectively, these results showed that transferring the interaction potentials obtained
at lower values of χto systems at higher χimproves the agreement between atomistic
and CG simulations for the hydrophobic sites, without altering the interactions between
other CG sites. Accordingly, the CG interaction potentials obtained at χAc and χBut
of 16% and χHep of 8% were applied in the rest of this study to model the systems with
higher χ. Notably, this approach has the additional advantage of rendering the model
more versatile for constructing polysaccharides with varying degrees of modification and
that can be arranged in different patterns.
Next, the performance of the chitosan model for higher degrees of polymerization and
for systems with a higher water content was tested. As a first step, the transferability
to systems with higher water content, i.e. lower chitosan concentration, was explicitly
tested for chitosan with DP=16. As expected from related systems27 , the RDFs ob-
tained for chitosan with 32 water molecules per chitosan monomer using CG interactions
transferred from a system with 12 water molecules per chitosan monomer, were found to
very closely resemble the native CG model produced at the lower chitosan concentration
in Figure 3.9.
Then, a system of larger polysaccharides with DP=50 and 32 water molecules per
monomer is constructed and simulated with the same CG model with respect to DP
and water concentration. For most interactions, the RDFs of the longer chains, shown
58
Figure 3.9: RDFs for 16% acetylation comparing all-atom results with DP=16 and 32
waters/monomer (black), CG results with DP=16 and 32 waters/monomer (red) and
transferred CG interaction potential from DP=16 and 12 waters/monomer (blue).
in Figure 3.10, strongly resemble those obtained for the 16-mers. Only the short range
peaks of the solute-solute interactions, which correspond to bonded neighbors, increased
in intensity for the higher DP, because there are more bonded neighbors.
3.5 CG model for drug
A CG FFs was developed for system 5 and 6 by following the same procedure as described
in the method section. A feasible mapping schemes were used to map DOX and GEM to
obtain CG representation as compared to all-atom representation. For DOX, 9 distinct
bead types were used to have a planar structure in the chosen mapping which resemble
59
Figure 3.10: RDFs for 16% acetylation comparing all-atom results with DP=16 and 12
waters/monomer (black), CG results with DP=16 and 12 waters/monomer (red) and
CG results for DP=50 and 32 waters/monomer (blue).
to its atomistic structure as shown in Figure 3.2, with a total number of 11 beads.
Different combination of bonds, angles and dihedrals were tested. In the final model,
the bonds DA-DB, DA-DCi, DB-DCi+1, DCi-DCi+1, DCi-DDi, DCi+1-DDi+1, DDi-DE,
DDi+1-DE, DE-DF, DDi+1-DG, DG-DH, DG-DI, and DI-DH are used. The four angles,
DC2-DD2-DG, DD1-DE-DF, DD2-DG-DH and DE-DD2-DG and five dihedral, DC1-
DD1-DD2-DG, DC2-DD2-DG-DH, DD2-DG-DH-DI, DF-DE-DD2-DG, and DE-DD2-
DG-DH are used. For the CG representation of GEM, four distinct bead types were
used as presented in Figure 3.2. Six different bonds namely GA-GB, GB-GC, GC-GD,
GA-GC and two explicit 1-3 intra-molecular interactions GA-GD and GB-GD were used.
Simulation snapshot of DOX in water and GEM in water are shown in Figure 3.11.
Like for the carbohydrates, the CG-FFs was validated by comparing the RDFs for
60
Figure 3.11: Simulation snapshots of DOX and GEM with 10 water molecules/molecules.
The simulation snapshots show DOX in purple, GEM in green and water molecules as
blue dots for a) DOX in water (b) GEM in water
both the drugs. RDFs show overall good agreement with the atomistic RDFs as shown
for two example in Figure 3.12.
Figure 3.12: RDFs for interaction of DOX and GEM with water (a). DA bead from
DOX with water, (b) GA bead from GEM with water
The other RDFs show a similar resemblance with atomistic RDFs. The developed
model was used to provide drug-drug CG interaction potential for futher simulation.
Along with that, system 1-4 were simulated with 10 DOX and 10 GEM molecules sepa-
61
rately to obtain drug-chitosan CG interaction potential.
3.6 Conclusion
In this chapter, a CG model for chitosan with three different hydrophobic modifica-
tions was introduced and analyzed. First, all-atom simulations of the different chitosan
solutions were performed and CG interactions potential were obtained from MS-CG
procedure for non-bonded and BI for bonded interaction to study the polymer network
at larger length and time scale. The CG model can reproduce structural data of the
all-atom systems, such as details of the RDF, angle distribution and end-to-end dis-
tances well. In addition we showed that it can be transferred to systems with high water
content and to molecules with high degree of hydrophobic modification. A CG model
of the two drugs namely DOX and GEM was also proposed, showing good resemblance
with their atomistic counterpart.
62
Chapter 4
Effect of hydrophobic modifications on
chitosan hydrogel properties
4.1 Introduction
This chapter describes the application of the CG model developed in Chapter 3 to ob-
tain equilibrated structures of chitosan hydrogels under different conditions. We had
analyzed the effect of modification type, degree and pattern as well as water content
in the gel on the network structure. The network structure were equilibrated and then
quantitatively characterized in terms of the pore sizes, end to end distance of the poly-
mers and average number of contacts between the modifications. Finally the effect of
pattern of modification on the network structure was also analyzed and the structure
formed compared to those corresponding to the uniformly-spaced modifications. The
CG-FFs developed in previous chapter 3 was used for longer polymer i.e. DP =50. Be-
cause no correlation between the gel fraction and the molecular weight of chitosan exist
and simulating chitosan chains with DP >50 would lead to long computational time,
hence considered unnecessary.
In chitosan hydrogels, usually both modification and their patterns are governed
by the monomers and the reagents utilized for chemical modification, as well as by
the structure of the polymer chain in solution as it evolves through the course of the
chemical modification(χ). Usually, the anhydrides used in the modification process are
amphiphilic molecules, their tendency to form aggregates will increase with increasing
hydrophobic fraction. Whereas χcan be easily measured, evaluating the modification
patterns is much more challenging, and it is often assumed to be random. However, other
63
distributions of the modification groups, an evenly spaced or a blocky pattern can also
be envisioned, and their presence can have profound effect on the macroscopic behaviors
of the polymer. A block-type polymer pattern is more favored during hydrophobic
modifications of the anhydrides (like acetic, butanoic and heptanoic). As, the probability
of creating a modification next to an already modified(hydrophobic monomers) site may
becomes larger, leading to clustered modificationss in aqueous acidic media.
It is likely that both χand the modification pattern, in fact, govern the aggregation
of polymer chains in solution; in particular, chain aggregation is likely to be promoted
by blocky modification patterns as opposed to a random distribution along the polymer
chain71, 72, 73 . While such an effect cannot be easily controlled or measured in exper-
iments, computational models have full control over the distribution of modifications,
and can therefore explore the effect of different patterns on the hydrogel structure. To
evaluate this effect of the modification pattern on the network structure and the molec-
ular interactions, chitosan chains modified with two different patterns, namely, evenly
spaced modifications (i.e., two neighbor modification groups are separated by a number
of unmodified monomers) and blocky (i.e., clusters of four modified monomers separated
by a number of unmodified monomers), were constructed to be simulated and analyzed.
4.2 Methods: CG network model
The larger hydrogel systems were also initially constructed with atomistic resolution,
following the same steps as for the shorter chains ( described in Chapter 3 ). Starting
structures for the CG simulations were obtained from 10 ns atomistic NPT simulations,
to obtain the corresponding box-sizes. Then the systems were mapped to their CG
representation using the VOTCA software59, 60 . The large hydrogel structures were sim-
ulated for 100 ns in NVT. Data from the last 10 ns was used for analysis.
CG simulations of the following systems were performed with Gromacs 4.6.441 :
i. 50 chitosan chains with DP=50 and 80,000 water beads, corresponding to 1600
water molecules per polymer chain (32 water molecules per monomer ), with evenly
spaced modifications.
ii. 50 chitosan chains with DP=50 and 80,000 water beads, with modifications grouped
in blocks of four.
64
iii. 20 chitosan chains with DP=50 and 100,000 water beads or 5000 water molecules
per polymer chain ( 100 water molecules per monomer ), with evenly spaced mod-
ifications.
iv. 20 chitosan chains with DP=50 and 100,000 water beads, with modifications
grouped in blocks of four.
The hydrogel structures were characterized by the distribution of pore sizes, following
the same protocol mentioned in chapter 2. For each distribution, three snapshots of the
network structures at 5 ns intervals were analyzed.
4.3 Low water content
Simulation snapshots of dense chitosan networks with 50 polymers with DP=50 and
80,000 CG water molecules are presented in Figure 4.1. These polymers have evenly
spaced acetyl, butyl and heptyl modifications. For all systems, the chitosan polymers
appear to be overall uniformly distributed throughout the box. Analysis of the pore sizes
in these networks results in a similar picture. The pore size distributions, also shown
in Figure 4.1, reveal no differences between the different modifications, all showing a
relatively narrow distribution with most pore diameters at 1.2 nm.
Though no measurable large scale differences in the network structure are found,
inspection of the simulation snapshots suggest that the more hydrophobic modifications
form local clusters. For a more local picture of the molecular interactions within the
networks we have characterized the number of contacts formed between the hydrophobic
modifications, as summarized in Table 4.1. For the acetyl modifications the all contacts
of M beads have been counted, while butyl and heptyl the number of contacts formed
by MB or MC beads, respectively, with any other modification beads were counted.
As one would expect, the total number of contacts increases with higher χand for
larger modifications, both of which correspond to a higher total number of modification
beads in the network. For acetyl chitosan, the M beads form on average only 28 ±7
contacts for χ= 16% but 273 ±2 for χ= 50%. Similarly the contacts of the butyl MB
beads increase from 433 ±8 for χ= 16% to 1002 ±15 for χ= 32%, and the heptyl
MC beads form 306 ±9 contacts at χ= 8% and 980 ±18 contacts for χ= 24%. Put
in relation to by the total number of modifications in the system, the trends remain
similar: both butyl and heptyl form significantly more contacts per modification than
acetyl, and in all cases the number of contacts increases with χ.
65
Figure 4.1: Simulation snapshots and pore-size distributions of chitosan networks with 50
polymers of DP=50 and 32 water molecules/monomer with (a-c) acetyl, (d-f) butyl and
(g-i) heptyl modifications. The simulation snapshots show the chitosan backbone (A,B,C
beads) in red, modifications (M, MA, MB, MC beads) in yellow and water mlecules as
blue dots for (a) χ= 16% acetylation (b) χ= 50% acetylation, (d) χ= 16% butylation,
(e) χ= 32% butylation, (g) χ= 8% heptylation and (h) χ= 24% heptylation.
66
low water high water
Modification total per modification total per modification
16% Acetylation 28 ±7 0.07 ±0.02 4 ±2 0.025 ±0.012
50% Acetylation 273 ±2 0.22 ±0.002 46 ±10 0.092 ±0.02
16% Butylation 433 ±8 1.08 ±0.02 165 ±3 0.41 ±0.01
32% Butylation 1002 ±15 1.25 ±0.01 561 ±10 1.03 ±0.06
8% Heptylation 306 ±9 1.53 ±0.006 115 ±5 1.44 ±0.06
24% Heptylation 980 ±18 1.63 ±0.015 233 ±7 0.97 ±0.03
Table 4.1: Number of contacts between modification beads formed in the chitosan net-
works
Finally, the end-to-end distance can provide insight as to the effect of the modifica-
tions and of the interactions within the network on the conformations of the polysac-
charides. The data for single polymers with different modifications and different χis
summarized in Table 4.2 and simulation snapshots in Figure 4.2. The addition of acetyl
modifications barely affect the conformational space of the polysaccharide. Butyl and
heptyl modifications on the other hand lead to more compact polymers, with decreasing
end-to-end distances when the degree of substitution is increased, consistent with the
increased hydrophobicity of the molecules.
Evenly-Spaced Pattern Blocky Pattern
Modification single polymer network 50 chains network 20 chains single polymer network 20 chains
16% Acetylation 15.4 ±3.6 14.6 ±0.4 14.4 ±0.6 12.10 ±4.06 15.4 ±0.6
32% Acetylation 15.1 ±0.2 15.4 ±0.7 15.4 ±0.6
50% Acetylation 15.9 ±4.2 16.2 ±0.3 16.1 ±0.9 16.29 ±3.92 15.7 ±0.6
16% Butylation 10.2 ±5.1 12.9 ±0.5 12.7 ±0.8 13.68 ±3.65 13.4 ±0.6
24% Butylation 11.8 ±0.3 10.4 ±0.5 13.0 ±0.4
32% Butylation 8.7 ±3.4 9.0 ±0.3 5.6 ±0.4 10.84 ±3.66 10.2 ±0.3
8% Heptylation 13.2 ±3.7 13.5 ±0.5 13.1 ±0.7 10.42 ±3.78 13.0 ±0.4
24% Heptylation 11.7 ±3.7 11.8 ±0.3 12.4 ±0.9 9.94 ±1.38 12.4 ±1.1
Table 4.2: Polymer end-to-end distance for single polymers in solution and in the net-
work. The errors represent one standard deviation.
The trends observed in the chitosan networks follow the same overall line as for
the single polymer. For acetylated networks, the average end-to-end distance of the
polymers increases slightly from 14.6 ±0.4 nm to 16.2 ±0.3 nm from 16% to 50%
acetylation, which suggests that the contacts in the hydrogel favor alignment and there-
67
Figure 4.2: Simulation snapshots of single chitosan polymers with DP=50 and 100000
water molecules with (a-d) acetyl, (e-h) butyl and (i-l) heptyl modifications. The sim-
ulation snapshots show the chitosan backbone (A,B,C beads) in red, modifications (M,
MA, MB, MC beads) in yellow (a) χ= 16% acetylation with evenly-spaced pattern (b)
χ= 50% acetylation with evenly-spaced pattern, (c) χ= 16% acetylation with blocky
pattern (d) χ= 50% acetylation with blocky pattern, (e) χ= 16% butylation with
evenly-spaced pattern, (f) χ= 32% butylation with evenly-spaced pattern, (g) χ= 16%
butylation with blocky pattern, (h) χ= 32% butylation with blocky pattern, (i) χ= 8%
heptylation with evenly-spaced pattern and (j) χ= 24% heptylation with evenly-spaced
pattern, (k) χ= 8% heptylation with blocky pattern and (l) χ= 24% heptylation with
blocky pattern.
68
fore more linear conformation in the polymers. Both butyl and heptyl modified chitosan
polysaccharides have more compact conformations with higher degrees of modification.
Thus overall, the observed network structure of the dense systems is very similar
across different modification types, degrees and patterns. However, the higher number
of contacts formed and the differences in the polymer conformations reflected in the
altered end-to-end distances suggest that the modifications are likely to affect hydrogel
properties under the right conditions.
4.4 High water content
One factor that can influence the network structure significantly, and that limits polymer
flexibility and the sizes of the pores, is the high chitosan density and low water content
in the systems above. The water content in real hydrogels is often much larger, reaching
weight fractions up to to 0.98-0.99.
To better understand what effect the water content has on the hydrogel properties,
we have simulated hydrogels for all modifications with a higher water content. A 20
chitosan chains with DP = 50 and 100,000 CG water beads, i.e. 100 water molecules
per monomer. This corresponds to a weight fraction to 0.89. Representative simulation
snapshots from the end of the 100 ns trajectories are shown in Figure 4.3. From the
snapshots, the overall structure of the acetylated networks still appears to be indepen-
dent of χ. Similar to the denser systems, chains are evenly distributed throughout the
simulation box. For the butyl-chitosan on the other hand the network-structure changes
drastically between χ= 16% and χ= 32% of modification, forming dense clusters of
hydrophobic modifications, surrounded by the polymer backbone, and connected by in-
dividual polymer strands. The result is a cluster/channel network featuring large pores
between the dense chitosan clusters as shown in Figure 4.3e).
This morphological change is also clearly visible in the pore size distributions shown
in Figure 4.3. While at 16% butylation, the pore size distribution looks very similar to
those observed in the acetylated networks, with the majority of pore diameters between
now 1.5 and 3.0 nm, in the 32 % butylated chitosan network, the pore sizes have markedly
shifted to pore diameters centered around 7.0 nm.
As can be expected from the smaller number of modification beads in the system,
and the larger volume fraction occupied by water, the total number of contacts between
modification beads in the network has decreased. To account for the smaller total number
69
Figure 4.3: Simulation snapshots and pore size distributions of chitosan networks with 20
polymers of DP=50 and 100 water molecules/monomer with (a-c) acetyl, (d-f) butyl and
(g-i) heptyl modifications. The simulation snapshots show the chitosan backbone (A,B,C
beads) in red, modifications (M, MA, MB, MC beads) in yellow and water molecules as
blue dots for (a) χ= 16% acetylation (b) χ= 50% acetylation, (d) χ= 16% butylation,
(e) χ= 32% butylation, (g) χ= 8% heptylation and (h) χ= 24% heptylation.
70
of modifications in the system, the number of contacts has again been normalized by
the number of modifications in the system. The number of contacts and the number
of contacts per modifications for the high-water networks are also listed in Table 4.1.
In the acetylated networks, and for butyl at low χ, the average number of contacts
each modification forms has decreased with the higher water fraction. For χ= 32%
butylation and heptyl on the other hand, the contacts per modification are almost as
high as in the denser networks, indicating their role in forming adhesive contacts.
The end-to-end distances of the polymers in the high water content gels, listed in
Table 4.2, also show the same trends as before. Especially for the butylated chitosan, the
changes with increasing values of χhave become much more pronounced in the hydrated
networks, and the end-to-end distance decreases significantly, changing from 12.7 ±0.8
nm at χ= 16% to 5.6 ±0.4 nm at χ= 32%. These very compact conformations of
the polymers are required to accommodate the aggregation of the butyl modifications
to clusters. In addition, the polymer backbones tend to wrap around these hydrophobic
clusters to form micelle-like structures.
In contrast, a similar morphological transition does not take place in heptanoyl-
chitosan systems, despite the even higher hydrophobicity of modification. Although
small clusters of modification beads are visible in the simulation snapshots (Figure ??,
the resulting pore size distribution does not significantly change. Consistent with the
absence of hydrophobic clusters heptanoyl-chitosan did not show any notable differ-
ence in the end-to-end distances. This different behavior of butanoyl- and heptanoyl-
chitosan networks originates from the association between the hydrophobic modification
groups, which is much stronger for the larger heptanoyl moieties. As a result, the initial
heptanoyl clusters are too stable to allow a rearrangement of the network into a clus-
ter/channel morphology, which requires dissociation and rearrangement of some of these
initial contacts. This suggests that the network structure heptanoyl-chitosan is likely
to be more rigid, and less amenable to adjust to different conditions in the surrounding
aqueous environment, such as variations in ionic strength and pH.
4.5 Influence of modification pattern
Whereas the degree of acetylation or further modification of chitosan is typically well
characterized, the distribution of the modifications on the chain is not easily determined
experimentally and is likely to depend on the preparation process71, 72, 73 . For chitosan,
71
typically a random distribution is assumed, however, it has also been suggested that
deacetylation under heterogeneous conditions favors the formation of a blockwise pat-
tern74, 75, 76 .
A real random pattern will contain a mixture clusters of different sizes, especially
for higher degrees of modification, and single modifications. To gain a clearer under-
standing of the relation between modification pattern and network properties, instead of
simulating random patterns, we opted to compare the effect of individual, evenly spaced
modifications as described in the previous sections to small modification blocks of four
consecutive modifications (see Figure 4.4).
Figure 4.4: Effect of modification pattern: a) scheme of the evenly spaced and blockwise
modification pattern; b) pore size distribution for χ= 16% butylation with the two
patterns; c) and d) simulation snapshots of for χ= 16% butylation with c) evenly
spaced and (d) blockwise modification.
The tendency of the hydrophobic modifications to interact with each other, is al-
ready observed in the RDFs shown in Figure 3.8 as well as in the clusters formed in
72
Figure 4.5: Simulation snapshots and pore-size distributions of chitosan networks with
50 polymers of DP=50 and 32 water molecules/monomer with (a-c) acetyl, (d-f) butyl
and (g-i) heptyl modifications with blockwise modification pattern. The simulation
snapshots show the chitosan backbone (A,B,C beads) in red, modifications (M, MA,
MB, MC beads) in yellow and water molecules as blue dots for (a) χ= 16% acetylation
(b) χ= 50% acetylation, (d) χ= 16% butylation, (e) χ= 32% butylation, (g) χ= 8%
heptylation and (h) χ= 24% heptylation.
73
Figure 4.6: Simulation snapshots and pore-size distributions of chitosan networks with
20 polymers of DP=50 and 100 water molecules/monomer with (a-c) acetyl, (d-f) butyl
and (g-i) heptyl modifications with blockwise modification pattern. The simulation
snapshots show the chitosan backbone (A,B,C beads) in red, modifications (M, MA,
MB, MC beads) in yellow and water mlecules as blue dots for (a) χ= 16% acetylation
(b) χ= 50% acetylation, (d) χ= 16% butylation, (e) χ= 32% butylation, (g) χ= 8%
heptylation and (h) χ= 24% heptylation.
74
the butylated and heptylated chitosan networks. A block pattern of modifications al-
lows the modification beads to interact with several other modifications simultaneously,
and thereby facilitates the formation of hydrophobic clusters while requiring less defor-
mation of the polymer backbone. Thus, an enhancement of interactions between the
modifications is expected.
Low Water High Water
Modification Total per modification Total per modification
16% Acetylation 37 ±9 0.09 ±0.02 5 ±3 0.03 ±0.02
50% Acetylation 319 ±23 0.26 ±0.02 56 ±10 0.11 ±0.02
16% Butylation 456 ±12 1.14 ±0.03 308 ±3 1.92 ±0.02
32% Butylation 1034 ±20 1.23 ±0.03 507 ±16 1.58 ±0.05
8% Heptylation 399 ±11 2.00 ±0.06 151 ±6 1.89 ±0.08
24% Heptylation 1383 ±24 2.35 ±0.04 491 ±15 2.04 ±0.06
Table 4.3: Number of contacts between modification beads formed in the chitosan net-
works with modifications grouped in blocks of four
A comparison of the polymer data for the two patterns is summarized in Table 4.2.
For acetylated chitosan, a comparison to the evenly spaced modifications shows no dis-
cernible effect of the blockwise distribution on the polymer properties, or network struc-
ture.
The average number of contacts between modifications increases compared to a uni-
formly distributed pattern in all the cases as shown in Table 4.3 , because a modification
bead can more easily interact with several others in the immediate vicinity in the block
pattern. For the end-to-end distances, overall the same trends are observed for the block
chitosan, though they become less pronounced in the dense system as shown in Table 4.2
for high-water and low water content respectively.
The most striking difference between the two modification patterns, is observed for
the onset of the network morphological transition to the cluster-pore network, which is
markedly shifted to lower values of χfor a blockwise distribution of butyl groups. For
other cases, it shows similar behaviour as compared to uniformly spaced modifications.
At high (χ= 32%) level of modification, the networks for blockwise and evenly dis-
tributed modifications show a similar structure and and pore size distribution. However,
for the block pattern, the distribution of pore sizes has already shifted to larger pore
diameters, as compared to both the evenly spaced modifications (Figure 4.4b).
75
A change in the hydrogel structure is more easily accessible to experimental quan-
tification than the distribution of modifications, for example by comparing the diffusion
of different probe molecules through the system, as described below, and may therefore
be helpful to distinguish modification patterns. However, for such a purpose, a detailed
analysis of the effects of various cluster sizes as well as truly random distributions would
be required.
4.6 Full range of substitution: 0% and 100%
To complete the structural picture, we have simulated the limiting cases of χ= 0%
and χ= 100% acetylation. The latter corresponds to unmodified chitin, whereas the
former represents fully deacetylated chitosan, which is, in practice, hard to produce and
therefore not often used, so that typically 5% to 15% of monomers retain their acetyl
groups. Experimental studies report that full deacetylated chitosan forms hydrogels with
homogeneous networks77 . Similarly, the fully deacetylated chitosan forms a uniform
network structure spanning the simulation box as shown in the simulation snapshot in
Figure 4.7a) , and the distribution of pore sizes and end-to-end distances that is very
similar to those found in the acetylated chitosan networks.
Figure 4.7: (a) Simulation snapshot and (b) pore-size distribution of deacetylated chi-
tosan networks with 20 polymers of DP=50 and 100 water molecules/monomer. The
simulation snapshots show the chitosan backbone (A,B,C beads) in red.
The fully acetylated polymers on the other hand self-assemble into form thick fib-
rils, as shown in Figure 4.8a). As a result of the alignment, the individual polymers
76
are straightened out and the average end-to-end distance increases to 20.4 ±0.4 nm.
Chitin is known to form crystalline fibrils with anti-parallel alignment of the polymers.
Figure 4.8b) shows a fibril segment color-coded by the polymer direction. In this fib-
ril, as well as the rest of the system, only about 50% of the polymers have formed an
anti-parallel alignment, the rest aligned parallel to each other. In addition, the self-
assembled fibrils are more twisted than the chitin structure. The crystal structure of
α-chitin has been resolved61, 62 and can be compared to the self-assembled fibrils in our
simulation model. Nevertheless, the contacts formed between the polymers with anti-
parallel alignment show a remarkable similarity to those in the α-chitin structure, as
shown for example in Figure 4.8c), which superimposes a segment of self-assembled fib-
ril on the crystal structure. Thus, the CG model, despite being derived at quite different
conditions, and low acetylation is able to reproduce the right crystal packing at least
locally, which speaks for the wide rage of applicability of the model. Reproducing the
overall formation of crystal fibers with the right orientation cannot be expected in an
unbiased simulation of such large molecules, which are kinetically trapped, once they
have aligned in either orientation78 .
Figure 4.8: (a) Simulation snapshot of 100% acetylated chitosan networks with 20
polymers of DP=50 and 100 water molecules/monomer showing the chitosan backbone
(A,B,C beads) in red, modifications M beads in yellow, (b) alignement of polymer in
blue as non-reducing end at the top and orange as non-reducing end at the bottom while
(c) shows the structure of two antiparallely aligned polymers superimpose to α-chitin
structure from Ref61, 62 .
77
4.7 Conclusion
We have studied the morphologies of hydrogels formed by chitosan polymer with different
modifications i.e. acetyl-, butanoyl-, heptanoyl and different water concentration with
CG MD simulations. The structures formed were found to depend significantly on these
material parameters. At low water concentration, the network structure does not show
a measurable difference in hydrogel properties such as pore-size distribution, end-to-end
distance and average number of contacts. The polymers form a uniform network in all
cases. However, in case of high water content, the acetyl modified hydrogels still form a
uniform network while with butanoyl, polymers aggregate around hydrophobic clusters
formed. In heptyl modified hydrogel, polymer show strong association of the heptyl
modification. However, they are so strong that once formed they are unable to rearrange
so that no cluster-channel morphology can form. The pore size distribution shows shifts
in the pore size, confirming the formation of polymer aggregates. The distribution of
modifications across the polymer chains can also cause morphological changes in the
network structure. Two modification patterns namely, a block modification pattern
(four consecutive modified monomer) and evenly-spaced were analyzed in detail. Similar,
network structures were obtained as for the evenly-spaced pattern for most conditions.
However, the formation of hydrophobic clusters in the network becomes more favorable,
so that the transition from the uniform network to the cluster-channel morphology sets
on earlier in butyl modified chitosan. Finally, the minimum (0%) and maximum (100%)
degree of acetylation cases were also studied. The 0% acetylation form a uniform network
structure while 100% acetylation forms fibrils. The fibril are twisted and 50% anti-
parallel aligned that resemble with the structure of αchitin.
78
Chapter 5
Effect of chitosan hydrogel properties
on drug diffusion
5.1 Introduction
Chitosan hydrogels have been extensively utilized as materials for tissue engineering,
wound healing, and drug delivery79, 80, 81 . Their optimization as drug delivery carri-
ers for clinical applications, however, is extremely laborious due to the variety of tun-
able parameters controlling the structure and transport properties of these materials.
Computational models capable of predicting the transport properties of drugs across
chitosan-based hydrogels have the potential of accelerating pre-clinical development and
translation. Critical towards accurate modeling of drug diffusion is the development of
a detailed model connecting physicochemical properties and architecture of the chain
network and drug properties. In this work, we adopted DOX and GEM, established
chemotherapeutic drugs, as model molecules to study diffusion through modified chi-
tosan hydrogels. GEM is a small (263.2 g/mol), hydrophilic and electrically neutral
molecule, whereas DOX is a larger (543.5 g/mol) amphiphilic molecule. CG models of
DOX and GEM were initially prepared through the mapping shown in chapter 3, and
drug-drug, drug-solvent, and drug-chitosan interactions were obtained from initial all-
atom simulations as presented in chapter 3.
In this chapter, first an atomistic analysis was conducted to study the type of interaction
between the chitosan polymer-drugs. The CG model of chitosan introduced in chapter 2
& 3 was designed in a way that can afford the delivery for any desired drug combination
with a synergistic molar ratio and kinetics. As an example for the significance of hydro-
79
gel morphology, the diffusion of two example drug molecules GEM and DOX through
the different gel morphologies was also simulated. Diffusion trends were obtained for
both the drugs. To study the effect of multiple drugs interacting with each other on
diffusion trends, dual drug migration through the network was also studied.
5.2 Analysis of drug-chitosan interactions at all-atom
resolution
All-atom simulations of the drug-loaded chitosan hydrogels were analyzed to obtain a
molecular-level understanding of (i) the interactions between the drug molecules and the
modified chitosan chains and (ii) the dependence of these interactions on the type and
degree (χ) of modification.
Atomistic systems have a high density of chitosan monomers and low water content,
which were necessary to obtain a sufficient sampling of the forces needed to implement
the CG procedure as described in Chapter 2. This high density, however, results in
many non-bonded interactions between the polymer and drug molecules. It must be
also noted that the frequency of these interactions in real hydrogels, where the water
content is significantly higher, is much lower. Nonetheless, the analysis of these inter-
actions at the atomistic scale provides insight into the relative contribution of different
types of interaction between the drug and polymer. The interaction energies among the
drug molecules, the chitosan backbone, and the modification groups were separated into
electrostatic and Lennard-Jones (LJ) contributions and the ability to form hydrogen-
bond was analyzed. The resulting interactions between DOX and modification groups
and between DOX and the chitosan backbone as a function of χare reported in Fig-
ure 5.2. The analogous results for GEM are shown in Figure 5.3. The hydrogen-bond
interactions for DOX modification and DOX backbone as a function of χare reported in
Figure 5.4. The hydrogen-bond interaction for GEM modification and GEM backbone
as compared to DOX is negligible.
We observed that the interactions between the drug molecules and chitosan were con-
sistently dominated by the LJ-type, particularly for DOX as compared to GEM. DOX
interacts both with the chitosan backbone and the modifications; The daunosamine moi-
ety in DOX (ring structure shown in Figure 5.1) aligns with the pyranose rings in the
chitosan backbone resulting in the formation of multiple hydrogen bonds (Figure 5.1a
and b). However no alignment were observed between GEM and the chitosan back-
80
Figure 5.1: All-atom simulations depicting the interactions between the interactions
between DOX and GEM ( red for oxygen, blue for nitrogen, fluorine for pink and white
for hydrogen) and chitosan (grey:Glucosamine and yellow: N-acyl-glucosamine) for (a)
acetyl modified chitosan with DOX, (b) butanoyl modified chitosans with DOX, (c)
acetyl modified chitosan with GEM, and (d) butanoyl modified chitosans with GEM.
bone or modification as shown in Figure 5.1c and d) and doesnot result in multiple
hydrogen bond as well. As the χof acetylation increases from 16% to 50%, the number
of contacts between the DOX and the modifications increases from 55 to 170. With
butanoyl-modified chitosan, the number of non-bonded interactions between DOX and
the modifications increases further to 225 at χ= 16% and 550 at χ= 32 %. The increase
81
Figure 5.2: Lennard-Jones and coulombic contribution to the DOX-modification group
interaction energy for (a) acetyl-chitosan and (b) butanoyl-chitosan; and Lennard-Jones
and coulombic contribution to the DOX-backbone interaction energy for (c) acetyl-
chitosan and (d) butanoyl-chitosans at different χ.
in the number of non-bonded interactions is reflected by the increase of interaction en-
ergy shown Figure 5.2. On the other hand, the dependence of the backbone interactions
with χis less clear. In particular, the number of non-bonded interactions decreases with
χfor acetyl modifications and increases slightly for butanoyl modifications. It should be
noted, however, that DOX molecules tend to aggregate into clusters(Figure 5.1a and b),
thereby reducing their ability to form interactions with the backbone, especially within
the timescale of the all-atom simulations, where equilibration of the aggregate size is not
82
Figure 5.3: Lennard-Jones and coulombic contribution to the GEM-modification group
interaction energy for (a) acetyl-chitosan and (b) butanoyl-chitosan; and Lennard-Jones
and coulombic contribution to the GEM-backbone interaction energy for (c) acetyl-
chitosan and (d) butanoyl-chitosans at different χ.
accessible. GEM shows overall weaker interactions with the chitosan chains as shown in
Figure 5.3. The LJ contribution to the interaction energy, in fact, is about 50% of that
observed with DOX, whereas the electrostatic contribution, is approximately the same
for both the drugs.
83
Figure 5.4: Hydrogen bond contacts for between DOX and modification groups in a (a)
acetyl- modified and (b) butanoyl-modified; and hydrogen bond contacts for between
DOX and the backbone in a (c) acetyl-modified and (d) butanoyl-modified chitosans at
different degrees of modification.
5.3 Dynamic properties of the CG system
To be able to reliably characterize their dynamics in the CG system, the diffusion co-
efficients of the drug molecules were calculated from the slope of the mean-squared
displacement (MSD) of the drug molecules with respect to time and compared to the
all-atom dynamics as listed in Table 5.1. As expected, the dynamics in the CG system
are significantly faster than in the all atom models, which is one of the factors leading
to the great efficiency of CG models in general. Comparing the diffusion coefficients
84
for water at atomistic and at CG resolution, Datom and DCG, respectively, the ratio is
τD=DCG
Datom = 6.4,i.e. dynamics in the CG system are 6 to 7 times faster than those
in the atomistic model. The value Datom = 2.74 ×10−5cm2/sis in good agreement with
diffusion constants reported in the literature for the TIP5P water model82 . Using the
same factor τDto scale the diffusion constants observed for DOX and GEM, the value
obtained for DOX Datom ≈0.1×10−5cm2/sis in reasonable agreement with the exper-
imentally determined value of Datom ≈0.21 ×10−5cm2/s83 . The atomistic values for
DOX and GEM on the other hand are significantly lower. It is possible that the small
box size and the high concentration of drug molecules contribute to this larger difference
and reduce the dynamics of DOX and GEM in the atomistic systems. Thus, especially
when the expected speedup in dynamics is taken into account, the CG model provides
a reasonable prediction of the diffusion coefficients of the two model drugs.
system D
water atomistic 2.741 ±0.164
water CG 17.530 ±0.733
DOX atomistic 0.056 ±0.0167
DOX CG 0.625 ±0.058
GEM atomistic 0.143 ±0.005
GEM CG 3.683 ±0.626
Table 5.1: Diffusion coefficients (10−5cm2/s) comparison for DOX and GEM with pure
water in atomistic and CG simulation respectively.
5.4 Simulation of single drug migration through modi-
fied chitosan hydrogels: Lower drug concentration
The all-atom simulations described in section 5.1 were performed to evaluate the molecular-
level interactions and elucidate the thermodynamic mechanisms by which the drugs
interact with the chitosan backbone and the modification groups.
The equilibrated chitosan hydrogel structures obtained in Chapter 4, were loaded
with 10 drug molecules of either DOX or GEM initially placed at random positions.
The concentration of drug molecules was chosen low to avoid strong drug-drug interac-
tion. In choosing the number of drug molecules, we also considered the trade-off between
85
ensuring reproducible simulations (higher drug loading) and avoiding the formation of
aggregates (lower drug loading). While, in fact, an insufficient number of drug molecules
can lead to poor statistical significance, especially when diffusing through structurally
non-homogeneous systems, excessive drug loading results in the formation of aggregates,
which would distort drug migration through the polymer network. However, similar
result were obtained for loading the system with 20 drug molecules as described be-
low in next section. Regarding the choice of distributing the drug molecules randomly
across the network, we note that the diffusion constants obtained by arranging the drug
molecules into different initial distributions showed the same trends and, in most cases,
agreed within the fitting error.
CG simulations of 10 ns were performed and the diffusion of the drug molecules
through the networks was monitored. The diffusion constants were calculated and are
reported in Figure 5.5. The error bars correspond to the fitting error due to the non-
linearity of MSD fitting. It should be noted that the diffusion constants did not signifi-
cantly change with simulations performed at longer time scales.
As anticipated, the differences in hydrogel morphology (homogeneous vs clusters/
channels) and physicochemical properties of the drugs and the modification groups result
in different trends of drug diffusion vs χ. GEM migrates through all networks with a
diffusion constant similar to that of free GEM in water, indicating that there is no
effect of the polymer on its diffusion (Figure 5.5 a & c). GEM molecules, which are
small and hydrophilic, form a relatively low number of nonbonded interactions with the
chitosan chains, irrespective of the type of modification and χ, as listed in Table 5.6.
This is also consistent with the size of GEM, which has a radius of gyration of 0.33
nm. The predominant pore diameters in the uniform hydrogels is on the order of 1.2
nm, so that GEM molecules can travel easily through the pore network provided that
there are no strong interactions with the polymers. Owing to the limited interaction
with the polymer chains, GEM travels easily through the water-filled channels in both
homogeneous and cluster/channel hydrogel morphologies. Only a slight reduction in
GEM diffusion is observed in butanoyl-chitosan at higher χ.
On the other hand, DOX molecules, which are larger and more hydrophobic, show
markedly different values of diffusion coefficient, and most notably, an inversion in the
diffusion trend with χbetween acetyl- and butanoyl- modifications. In acetylated-
chitosan networks, DOX diffusion decreases at higher χ, independently of the modi-
fication pattern (Figure 5.5 b). However, the network structure and pore size distri-
86
Figure 5.5: Drug diffusion constants vs. χfor single-drug migration across different
chitosan networks for evenly-spaced (black) and blocky (red) modification patterns: (a)
GEM and (b) DOX in acetyl-chitosan, and (c) GEM and (d) DOX in butanoyl-chitosan.
butions were found to remain unaltered with χ, both for the evenly-spaced and blocky
acetylated-chitosan networks , indicating that network morphology is unlikely to be the
cause of DOX decreased mobility. Rather, the decreased diffusion coefficients of DOX
at larger χcan be attributed to increased number of non-bonded interactions between
DOX and the modified chitosan molecules. As reported in Table 5.6, in fact, the num-
ber of non-bonded interactions between DOX and the acetyl moieties increases with χ
for both modification patterns, although the number of non-bonded interactions with
87
Drug Backbone, Low Mod, Low Backbone, High Mod, High
DOX
Acetyl
(evenly spaced) 217 ± 14 (93.5%) 15 ± 2 (6.5%) 239 ± 11 (83.9%) 46 ± 4 (16.1%)
Acetyl
(blocky) 225 ± 12 (84.9%) 40 ± 3 (15.1%) 221 ± 11 (79.8%) 56 ± 4 (20.2%)
Butanoyl
(evenly spaced) 219 ± 15 (78.5%) 60 ± 7 (21.5%) 316 ± 23 (83.6%) 62 ± 6 (16.4%)
Butanoyl
(blocky) 276 ± 13 (83.1%) 56 ± 7 (16.9%) 371 ± 18 (88.1%) 50 ± 7 (11.9%)
GEM
Acetyl
(evenly spaced) 26 ± 10 (89.7%) 3 ± 1 (10.3%) 33 ± 10 (80.5%) 8 ± 3 (19.5%)
Acetyl
(blocky) 30 ± 9 (88.2%) 4 ± 2 (11.8%) 41 ± 11 (80.4%) 10 ± 4 (19.6%)
Butanoyl
(evenly spaced) 34 ± 10 (89.5%) 4 ± 2 (10.5%) 169 ± 25 (87.6%) 24 ± 7 (12.4%)
Butanoyl
(blocky) 61 ± 11 (81.3%) 14 ± 5 (18.7%) 163 ± 16 (80.3%) 40 ± 6 (19.7%)
Figure 5.6: Number and (percentage of total) of non-bonded interactions observed be-
tween drug molecules and chitosan chains (backbone and modifications) in the different
chitosan networks over the drug molecules trajectories during the simulation; where
DOX molecules are treated as a group so that any given chitosan or modification site
can only contribute one contact. In the table, low represents χ= 16% for both sys-
tems, and high represents χ= 50% and 32% for acetyl-and butanoyl- modified systems,
respectively
the backbone remains constant. For example, a DOX molecule forms on average of 15
±2 interactions with acetyl-chitosan at χ= 16%, and up to 60±4 interactions with
acetyl-chitosan at χ= 50% with evenly spaced pattern. This indicates that hydrophobic
non-bonded interactions are the main cause of the slowed diffusion.
88
Figure 5.7: Mean-squared displacement (MSD) plot vs. time for diffusion of DOX
through the blocky butyl-modified chitosan network at χ= 32%, where DOX-Captured
refers to DOX that becomes entrapped within a cluster and DOX-Free refers to DOX
that remains in the pores of the cluster/channel morphology during the simulation.
Inversely, DOX diffusion through butanoyl-chitosan gels increases with χ, indepen-
dently of modification patterns (Figure 5.5 d). As χincreases, the butanoyl-chitosan
network undergoes a transition from homogeneous to cluster/channel morphology. CG
simulations of DOX migration through butanoyl-chitosan systems show dramatic dif-
ferences to the fate of DOX molecules depending on the morphology of the network,
i.e., homogeneous or cluster/channel. The migration of DOX through the homogeneous
butanoyl-chitosan network (low χ, Figure 5.9 c) is identical to that of DOX through
homogeneous acetyl-chitosan network (low and high χ, Figure 5.9 a and b). Through
clustered butanoyl-chitosan, instead, DOX molecules either adsorb onto/within the chi-
tosan clusters or travel freely through the large pores (Figure 8 d). This indicates that,
as in acetyl-chitosan systems, interactions between DOX and butanoyl-chitosan at high
χoccur (see Table 5.6). However, the DOX molecules freely migrating through the
large pore are responsible for the increase in the overall diffusion coefficients; this is
corroborated by the comparison in mean-squared displacement MSD vs. time for both
adsorbed and freely migrating DOX molecules as shown in Figure 5.7. It is also crucial
to note that the chitosan vs. water ratio and consequently the channel vs. cluster ratio
are much higher in the experimental systems than in the simulated networks. We there-
89
Figure 5.8: Snapshot of DOX migration through: (a) acetyl-chitosan networks at low
χ(16%); (b) acetyl-chitosan networks at high χ(50%); (c) butanoyl-chitosan network
at low χ(16%); (d) butanoyl-chitosan networks at high χ(32%) where the chitosan
backbone is represented by red beads, modifications are represented by yellow beads,
and DOX is represented by black beads
fore expect the actual hydrogel to mirror the increase in the diffusion coefficient with χ
observed in the CG simulations. The diffusion and release of DOX from hydrogels with
such a cluster-pore morphology will be therefore closely related to the fraction of the
molecules diffusing in the pores. This is determined by the partition coefficient Kcof the
molecules to the hydrophobic clusters, which is defined as Kc=ccl
cw=Ncl
Nw
Vw
Vcl . Therefore,
90
Figure 5.9: Snapshot of GEM migration through: (a) acetyl-chitosan networks at low
χ(16%); (b) acetyl-chitosan networks at high χ(50%); (c) butanoyl-chitosan network
at low χ(16%); (d) butanoyl-chitosan networks at high χ(32%) where the chitosan
backbone is represented by red beads, modifications are represented by yellow beads,
and GEM is represented by green beads
the fraction of DOX molecules stuck to the clusters is Ncl
Nw=KcVcl
Vwand depends both
on the affinity to the hydrophobic clusters, and the volume ratio of clusters and water
phase. Using octanol-water partition data as a rough estimate, Ncl
Nw≈4Vcl
Vw, so that in
gels with high water contents as is generally the case84 diffusion will be dominated by
molecules diffusing in the pores. Overall, the diffusion trends shows the same behavior
91
on compared with the experimental diffusion trends21 .
The faster dynamics enabled by the smoother energy landscapes featured in CG sim-
ulations may introduce some error in the calculation of the diffusion constants; however,
the comparison among the migration of drug molecules in different networks provides
a reliable evaluation of the influence of the network morphology and physicochemical
properties on drug transport. In the cluster/channel morphology, in particular, the dif-
ference between diffusional pathways of single-drug molecules becomes very pronounced,
with drug molecules experiencing sharper differences in the morphology of the medium
through which they diffuse. Collectively, these single-drug simulations indicate that the
hydrophobicity driven morphing of the hydrogel structure affects drug diffusion by two
opposing mechanisms: first, the larger pores favor the migration of the drug molecules;
and second, the nesting of the modification groups within the core of the clusters is
responsible for strong adsorption of the drug molecules that embed in the chitosan clus-
ters; notably, our simulations indicate that while both DOX and GEM are affected by
the first mechanism, only DOX undergoes the second mechanism. We anticipate that
the high water/polymer ratio of the real chitosan hydrogels makes the first mechanism
dominant over the second.
5.4.1 Effect of higher drug concentration on diffusion trends
To evaluate how the single drug diffusion trends change with high drug concentration,
20 DOX or 20 GEM molecules were inserted at random locations in the equilibrated
network structures for 16% & 50% acetylation and 16% & 32% butylation and their
motion through the chitosan hydrogels was analyzed. The results indicate that GEM
shows same diffusion trend as observed in Section 5.3 i.e. with diffusion constants that
are only minimally affected by interactions with the chitosan gels.
Similarly, the MSD curves for DOX shown in Figure 5.10 depend noticeably on
the type and level of modification, and the diffusion is slowed down for acetylation as
observed for the smaller number of DOX molecules. At χ= 16% the MSD of DOX
in the acetylated and butylated chitosan networks looks very similar, leading also to
approximately the same diffusion constants (see Table 5.2). However, when the level of
modification is increased, the MSD curve for 50% acetylation falls below that in 32%
acetylation, meaning that DOX diffusion decreases for larger χ. Inversely the diffusion
rate through 32% butylated chitosan has significantly increased, corresponding to a
diffusion constant that is more than twice as large as for χ= 16%.
92
Figure 5.10: Mean-square Displacement of DOX (a) Acetylated chitosan network (b)
butylated chitosan network for a uniformly-spaced modification pattern.
For further analysis the interactions between DOX molecules and the chitosan net-
work were analyzed by counting the number of contacts within 0.6 nm. As summarized
in Table 5.3 the number of interactions between DOX and the M beads increases for
χ=16% have 27 ±3while χ=50% have 83 ±5average number of contacts, whereas the
number of contacts with the chitosan backbone remains approximately constant at 400-
420 contacts over the range of χvalues. Thus the slow down of DOX diffusion through
the acetylated gel was observed due to the increase in interactions with the hydrophobic
modifications at higher χ, as explained in the previous section. In comparison, GEM
interacts much less with the chitosan network with only about 30-40 contacts in total,
consistent with its unhampered diffusion through the gel.
DOX GEM
Evenly-spaced Blocky Evenly-spaced Blocky
Modification Diffusion Const. Diffusion Const. Diffusion Const. Diffusion Const.
16% Acetylation 0.268 ±0.100 0.213 ±0.034 3.878 ±0.764 3.787 ±0.422
50% Acetylation 0.161 ±0.003 0.191 ±0.000 3.572 ±0.470 3.905 ±0.150
16% Butylation 0.257 ±0.042 0.150 ±0.050 3.906 ±0.149 3.588 ±0.464
32% Butylation 0.827 ±0.039 0.550 ±0.082 3.078 ±0.126 3.16 ±0.30
Table 5.2: Diffusion constant (10−5cm2/s) for DOX and GEM for acetylation and buty-
lation at different degree of modification.
However, the contact data presented in Table 5.3 shows that DOX still forms a large
number of contacts with the chitosan network. Going from 16% to 32% butylation,
93
DOX GEM
Evenly-spaced Blocky Evenly-spaced Blocky
Modification Avg. Contacts Avg. Contacts Avg. Contacts Avg. Contacts
16% Acetylation 27 ±3 52 ±4 6 ±2 7 ±3
50% Acetylation 83 ±5 98 ±5 18 ±4 21 ±5
16% Butylation 50 ±5 96 ±19 7 ±3 34 ±8
32% Butylation 89 ±9 124 ±10 58 ±10 80 ±18
Table 5.3: Average number of contacts for DOX and GEM for acetylation and butylation
at different degree of modification.
especially the number of interactions with the polymer backbone increases, following
same trends as shown by the low drug concentration diffusion trends. Overall, diffusion
trends remain consistent with the previous study21 . The trends for both the drugs are
independent of the drug concentration.
5.5 Simulation of dual drug migration through modi-
fied chitosan hydrogels
We finally sought to understand how the combination of multiple drugs affects the molec-
ular interactions and diffusion through the networks for different types and degrees of
modification. Such understanding is needed to gain control over the kinetics and mo-
lar ratio of release when using synergistic drug combinations. While some in silico
models have been developed to simulate the migration of single drugs through polymer
hydrogel networks85, 86 models for multiple drugs are much less common. Combining
drugs, especially with different physio-chemical properties, in fact, introduces consider-
able complexity related to drug-drug interactions and their outcome on drug-polymer
interactions. The combination of DOX and GEM has been proven to be a highly effective
drug regimen. Both are FDA- approved chemotherapeutics and have been studied in a
variety of drug delivery vehicles, such as micellar nanoparticles, PDCs, polymersomes,
mesoporous silica nanoparticles, and nanostructured lipid carriers87, 88, 89, 90 .
Following the same procedure as for the single drug, we performed CG simulations
of equimolar GEM and DOX migration through the different chitosan networks, and
calculated the diffusion constants for both drugs. Notably, the values derived from dual-
drug release (Figure 5.11) show rather different trends from those obtained with single
94
Figure 5.11: Drug diffusion constants vs. χfor dual-drug migration across different
chitosan networks: (a) GEM in acetyl-chitosan, (b) DOX in acetyl-chitosan, (c) GEM
in butanoyl-chitosan, and (d) DOX in butanoyl-chitosan.
drugs (Figure 5.5). In particular, the diffusion of DOX now remains almost constant
across the entire range of χfor both acetyl- and butanoyl-chitosan. The diffusion of
GEM is also markedly different. For alternated modification patterns, GEM diffusion
decreases with χthrough both acetyl- and butanoyl-chitosan. For blocky modification
patterns, instead, GEM diffusion increases with χthrough acetyl-chitosan and is almost
constant with χthrough butanoyl-chitosan.
These differences indicate that drug-drug interactions affect significantly their inter-
95
Figure 5.12: Radial distribution functions of GEM around the center of mass of DOX
in acetyl-chitosan networks with (a) evenly-spaced and (b) blocky modification, and
butanoyl-chitosan networks with (c) evenly-spaced and (d) blocky modification.
actions with the network. This is confirmed by the radial distribution functions of GEM
around DOX molecules, which indicate a strong tendency of the two drug molecules to
aggregate. As can be seen from the double-peak shape of the curves in Figure 5.12,
in fact, clusters of GEM around DOX and values <1 out to large distances are fre-
quently formed for both modification types, the trends of the drug-drug interaction with
increasing χare reversed for the blocky pattern.
96
5.6 Conclusion
Hydrogels constructed with native or modified polysaccharides and loaded with chemother-
apeutic drugs have been extensively studied, and a conspicuous number of them have
entered the clinical pipeline through the last decade91, 92, 93, 94 A growing segment in this
field is represented by polymer conjugates and hydrogels that deliver synergistic combi-
nations of drugs. In developing these systems for a given drug combination, the choice
of the modification groups, degree of modification, and initial drug loading are crucial
to ensure the therapeutic efficacy of the formulation. Empirical exploration of such
wide design space, however, is cumbersome. In this chapter, we have described the
development of a computational model that could serve as a powerful guide to pharma-
ceutical chemists in the identification of the design parameters that afford a schedule
and a ratio of drug release that ensure a successful therapeutic outcome. The proposed
model has been validated by comparison to experimental data by closely corresponding
systems and managed to accurately predict complex phenomena, such as the different
microscale morphologies present in hydrogels constructed with different types and de-
grees of modification, and the migration of not only one, but also two drugs through
these modified polymer networks. While focusing on hydrophobically modified chitosan
hydrogels and the GEM-DOX drug pair, this method is applicable to other polymer
substrates, modification moiety, and therapeutic payload.
97
Chapter 6
Conclusion
In this thesis, self-assembly of glucose and chitosan-based polysaccharides were studied
using atomistic and coarse-grained molecular dynamics simulation. In chapter 2, all-
atom simulations showed the effect of linkage, type, and various substitutions on the
flexibility of the polymer. The effect of the modification pattern on the molecules flex-
ibility was shown. To study these polymer systems at larger length and time scales,
coarse-grained models for each molecules were developed. The CG models were used to
study polymer aggregation for cellulose and chitosan based oligomers. The CG model
was able to predict fibril formation for cellulose and chitin. For methylated and flu-
orinated cellulose molecules, vastly different structures were observed, depending on
the pattern of modification. Whereas for methylated cellulose with alternate pattern
and fluorinated cellulose with blockwise pattern showed polymer aggregation in water
while, methylated cellulose with blockwise pattern and fluorinated cellulose with alter-
nate pattern were soluble in water. These results were verified experimentally and could
be explained based on the changed flexibility of the molecules and interactions of the
modifications.
In chapter 3, a CG force-field was developed. The transferability of the obtained model
corresponding to the degree of polymerization, degree of modification, and different
solute concentrations were demonstrated. Chapter 5 showed the diffusion of two anti-
cancer drugs namely DOX and GEM through the different networks. GEM diffuses
through the polymer networks similarly as through water. Its diffusion is independent
of the type and degree of modification. In the case of DOX, however opposing diffusion
trends with respect to the degree of modification are found for acetyl and butyl. Finally,
simulation with both types of drug showed different irregular diffusion trends as com-
98
pared to the single drug diffusion trends. Overall, the proposed combination of all-atom
and CG simulation as used here has demonstrated good predictive power and repre-
sents a reliable and predictive toolbox for understanding and predicting the properties
of carbohydrate aggregates, as can be used for example in pharmaceutical applications.
99
List of Figures
1.1 Chemical structures of (a). β-D- glucose, (b). β-D-Acetyl-glucosamine . . 10
1.2 Intra- and inter- molecular hydrogen bonds between adjacent monomer
and polymer chains .............................. 11
1.3 Different approaches to study the molecular systems with various resolu-
tion depending on the properties of interest. In a bottom-up model devel-
opment, coarser resolution simulations are guided by detailed level studies
like quantum mechanics. In top-down approaches, macroscopic proper-
ties are used to guide finer-resolution simulations like classical atomic or
coarse grained molecular simulations. ................... 13
1.4 Atomistic representation of (a) β-D-glucosamine solution and the corre-
sponding (b) coarse-grained representation. ................. 14
1.5 Bonded interaction potentials include (a) bond, (b) angle, (c) dihedral,
and (d) improper dihedral .......................... 16
1.6 Demonstrating force matching32 procedure by showing set of atomistic
forces, ⟨fI⟩and its corresponding resultant CG force FIfor single water
molecule. .................................... 18
2.1 Example of modified cellulose structures (a) Methylation of alternating
monomers (b) Alternative fluorination of alternating monomers increases
the hydrophilicity. (c) Chitosan containing amine groups increases the
positive charge in the polymer (d) Chitin, acetylation of all monomers
increases the hydrophobicity. ......................... 23
2.2 Definition of the dihedral angels ψ(C1,O4,C4,H4) and ϕusing the atoms
(H1,C1,O4,C4) ................................ 24
2.3 Simulation snapshots in (a). β-D 1-4 linked Glucose and (b) for β-D 1-6
linked Glucose. ................................ 26
100
2.4 Chemical structure of (a). β-D- glucose methyl modified at C3 atom, (b).
β-D-glucose fluorine modified at C3 atom. ................. 26
2.5 Analysis of ψdistribution for (a) alternating methyl modified cellulose
(b) blockwise methyl modified cellulose, (c)alternating fluorine modified
cellulose, and (d) blockwise fluorine modified cellulose. The residues are
numbered from the nonreducing end to the reducing end. ........ 27
2.6 Analysis of end-to-end distance (a) alternating methyl modified cellulose
(b) blockwise methyl modified cellulose, (c)alternating fluorine modified
cellulose, and (d) blockwise fluorine modified cellulose. The end-to-end
distance was monitored over 100ns. The residues are numbered from the
nonreducing end to the reducing end. .................... 28
2.7 Simulation snapshots of system 1 and 3 with DP = 6 and 2000 water
molecules. The simulation snapshots show the carbon atom in gray, oxy-
gen in red, hydrogen in white, and fluorine in pink. The hydrophobic
modification are encircled in yellow and hydrophilic modification in blue
for (a) glucose (b) alternated methyl modified glucose (c) block methyl
modified glucose (d) alternated fluorine modified glucose, and (e) block
fluorine modified glucose. ........................... 29
2.8 Chemical structure of (a).β-D- glucosamine(NH2), (b) β-D- glucosamine(NH+
3),
and (c). β-D-Acetyl-glucosamine ....................... 30
2.9 Simulation snapshots of system 4 and 5 with DP = 6 and 2000 water
molecules. The simulation snapshots show the carbon atom in gray, oxy-
gen in red, hydrogen in white, and nitrogen in blue for (a) N-acetyl glu-
cosamine (b) neutral glucosamine (c) charged glucosamine. ........ 30
2.10 Analysis of end-to-end distance and conformational maps of ϕand ψ
of chitosan (a) N-acetyl glucosamine (b) neutral glucosamine(NH2) (c)
charged glucosamine(NH+
3) obtained by MD simulations. The end-to-end
distance was monitored over 100ns.The residues are numbered from the
nonreducing end to the reducing end. The dihedral angles ϕand ψare
shown on x- and y- axes, respectively. ................... 31
2.11 All atom and coarse-grained representation of (a). β-D- glucose, (b).
methyl modified glucose (c). fluorine modified glucose (d) N-acetyl glu-
cosamine .................................... 33
101
2.12 Comparsion of AA and CG RDFs for interaction involved with CG site
type A and B with itself and water. (a) A-A (b) A-WAT, (c) B-B, and
(d). B-WAT .................................. 34
2.13 Comparison of AA and CG RDFs and their corresponding polymer snap-
shot for interaction involved for (a) M-M (b) M-WAT (c) snapshot, for
alternate methyl modification and (d) M-M (e) M-WAT and (f) snapshot
for block methyl modification. ........................ 35
2.14 Comparison of AA and CG RDFs of neutral and charged glucosamine
and their corresponding polymer snapshot for interactions (a) A-A (b)
A-WAT (c) snapshot for uncharged system (d) A-A (e) A-WAT and (f)
snapshot for charged system. ........................ 36
2.15 Simulation snapshots of (a) cellulose and (b) chitin networks with 100
polymers of DP=12 and 32 water molecules/monomer.The simulation
snapshots show the polymer backbone (A,B,C beads) in red, modifica-
tions (M beads) in yellow and water molecules as blue dots ........ 39
2.16 Simulation snapshots and pore-size distributions of cellulose networks
with 100 polymers of DP=12 and 32 water molecules/monomer. The
simulation snapshots show the cellulose backbone (A,B,C beads) in red,
methyl modifications (M bead) in yellow, fluorine (F bead) modifications
in blue, and water molecules as blue dots for (a) alternately methylated
cellulose (b) blockwise methylated cellulose (d) alternately fluorinated
cellulose, and (d) blockwise fluorinated cellulose. .............. 41
2.17 Simulation snapshots of single cellulose polymers with DP=10 and 40000
water molecules/monomer. The snapshots show the cellulose backbone
(A,B,C beads) in red, methyl and acetyl modifications (M bead) in yellow,
fluorine modification F in blue and water molecules as blue dots for a)
pure cellulose (b) 50% methylated cellulose with alternating pattern (c)
50% methylated cellulose with block pattern (d) 50% fluorinated cellulose
with alternating pattern (e) 50% fluorinated cellulose with block pattern
(f) chitin. ................................... 43
102
2.18 Simulation snapshots and pore-size distributions of chitosan networks
with 100 polymers of DP=12 and 32 water molecules/monomer. The
simulation snapshots show the chitosan backbone (A,B,C beads) in red,
charged glucosamine A in blue, and water molecules as blue dots for (a)
Glucosamine(NH2) (b) Glucosamine(NH+
3). ................ 44
3.1 Chemical structure of (a). N-acetyl glucosamine, (b) N-butyl glucosamine,
and (c) N-heptyl glucosamine. ........................ 47
3.2 All-atom and coarse-grained representations of (a) DOX, (b) GEM, (c)
unmodified glucosamine monomer, (d) acetyl-glucosamine, (e) butanoyl-
glucosamine, and (f) heptanoyl-glucosamine. ................ 49
3.3 Comparison of CG and atomistic RDFs of the distances between CG (A,
B, C, M for χAc=16% . Note: A, B, and C beads map the GlcN monomers,
whereas M map the modification(acetyl) group. .............. 51
3.4 Comparison of CG and atomistic RDFs of the distances between CG
(A, B, C, M for χBut=16% . Note: A, B, and C beads map the GlcN
monomers, whereas MA and MB map the modification groups ...... 52
3.5 Comparison of CG and atomistic RDFs of the distances between CG (A,
B, C, M for χHep=8% . Note: A, B, and C beads map the GlcN monomers,
whereas MA, MB, and MC map the modification groups ......... 53
3.6 RDFs for interaction involved with CG site type WAT with other beads
for (a),(b),(c) and (d) 16 acetylation and (e), (f) for 16 % butylation. . . 54
3.7 (a-c) Angle distributions in the atomistic and CG models; (d-f) molecular
conformations corresponding to the three free energy minima of the ϕ−ψ
dihedral angles. All-atom models are drawn as grey sticks, CG molecules
as red (ABC) and yellow (M) beads. .................... 55
3.8 RDFs of the distances between modification beads: (a) M-M beads in
acetyl-chitosan with χAc = 32%; (b) MA-MA and (c) MB-MB beads
in butanoyl-chitosan with χBut = 32%; and (d) MA-MA, (e) MB-MB,
and (f) MC-MC beads in heptanoyl-chitosan with χHep = 16%. The
RDFs obtained from the atomistic, native CG, and CG with transferred
potential models are in black, red, and blue, respectively. ......... 57
103
3.9 RDFs for 16% acetylation comparing all-atom results with DP=16 and 32
waters/monomer (black), CG results with DP=16 and 32 waters/monomer
(red) and transferred CG interaction potential from DP=16 and 12 wa-
ters/monomer (blue). ............................ 59
3.10 RDFs for 16% acetylation comparing all-atom results with DP=16 and 12
waters/monomer (black), CG results with DP=16 and 12 waters/monomer
(red) and CG results for DP=50 and 32 waters/monomer (blue). . . . . 60
3.11 Simulation snapshots of DOX and GEM with 10 water molecules/molecules.
The simulation snapshots show DOX in purple, GEM in green and water
molecules as blue dots for a) DOX in water (b) GEM in water ...... 61
3.12 RDFs for interaction of DOX and GEM with water (a). DA bead from
DOX with water, (b) GA bead from GEM with water .......... 61
4.1 Simulation snapshots and pore-size distributions of chitosan networks
with 50 polymers of DP=50 and 32 water molecules/monomer with (a-c)
acetyl, (d-f) butyl and (g-i) heptyl modifications. The simulation snap-
shots show the chitosan backbone (A,B,C beads) in red, modifications
(M, MA, MB, MC beads) in yellow and water mlecules as blue dots for
(a) χ= 16% acetylation (b) χ= 50% acetylation, (d) χ= 16% butyla-
tion, (e) χ= 32% butylation, (g) χ= 8% heptylation and (h) χ= 24%
heptylation. ................................. 66
4.2 Simulation snapshots of single chitosan polymers with DP=50 and 100000
water molecules with (a-d) acetyl, (e-h) butyl and (i-l) heptyl modifi-
cations. The simulation snapshots show the chitosan backbone (A,B,C
beads) in red, modifications (M, MA, MB, MC beads) in yellow (a)
χ= 16% acetylation with evenly-spaced pattern (b) χ= 50% acetylation
with evenly-spaced pattern, (c) χ= 16% acetylation with blocky pattern
(d) χ= 50% acetylation with blocky pattern, (e) χ= 16% butylation
with evenly-spaced pattern, (f) χ= 32% butylation with evenly-spaced
pattern, (g) χ= 16% butylation with blocky pattern, (h) χ= 32% buty-
lation with blocky pattern, (i) χ= 8% heptylation with evenly-spaced
pattern and (j) χ= 24% heptylation with evenly-spaced pattern, (k)
χ= 8% heptylation with blocky pattern and (l) χ= 24% heptylation
with blocky pattern. ............................. 68
104
4.3 Simulation snapshots and pore size distributions of chitosan networks
with 20 polymers of DP=50 and 100 water molecules/monomer with (a-c)
acetyl, (d-f) butyl and (g-i) heptyl modifications. The simulation snap-
shots show the chitosan backbone (A,B,C beads) in red, modifications
(M, MA, MB, MC beads) in yellow and water molecules as blue dots for
(a) χ= 16% acetylation (b) χ= 50% acetylation, (d) χ= 16% butyla-
tion, (e) χ= 32% butylation, (g) χ= 8% heptylation and (h) χ= 24%
heptylation. .................................. 70
4.4 Effect of modification pattern: a) scheme of the evenly spaced and block-
wise modification pattern; b) pore size distribution for χ= 16% butyla-
tion with the two patterns; c) and d) simulation snapshots of for χ= 16%
butylation with c) evenly spaced and (d) blockwise modification. . . . . 72
4.5 Simulation snapshots and pore-size distributions of chitosan networks
with 50 polymers of DP=50 and 32 water molecules/monomer with (a-c)
acetyl, (d-f) butyl and (g-i) heptyl modifications with blockwise modi-
fication pattern. The simulation snapshots show the chitosan backbone
(A,B,C beads) in red, modifications (M, MA, MB, MC beads) in yellow
and water molecules as blue dots for (a) χ= 16% acetylation (b) χ= 50%
acetylation, (d) χ= 16% butylation, (e) χ= 32% butylation, (g) χ= 8%
heptylation and (h) χ= 24% heptylation. ................. 73
4.6 Simulation snapshots and pore-size distributions of chitosan networks
with 20 polymers of DP=50 and 100 water molecules/monomer with (a-c)
acetyl, (d-f) butyl and (g-i) heptyl modifications with blockwise modifi-
cation pattern. The simulation snapshots show the chitosan backbone
(A,B,C beads) in red, modifications (M, MA, MB, MC beads) in yellow
and water mlecules as blue dots for (a) χ= 16% acetylation (b) χ= 50%
acetylation, (d) χ= 16% butylation, (e) χ= 32% butylation, (g) χ= 8%
heptylation and (h) χ= 24% heptylation. .................. 74
4.7 (a) Simulation snapshot and (b) pore-size distribution of deacetylated chi-
tosan networks with 20 polymers of DP=50 and 100 water molecules/monomer.
The simulation snapshots show the chitosan backbone (A,B,C beads) in
red. ....................................... 76
105
4.8 (a) Simulation snapshot of 100% acetylated chitosan networks with 20
polymers of DP=50 and 100 water molecules/monomer showing the chi-
tosan backbone (A,B,C beads) in red, modifications M beads in yellow,
(b) alignement of polymer in blue as non-reducing end at the top and
orange as non-reducing end at the bottom while (c) shows the structure
of two antiparallely aligned polymers superimpose to α-chitin structure
from Ref61, 62 ................................. 77
5.1 All-atom simulations depicting the interactions between the interactions
between DOX and GEM ( red for oxygen, blue for nitrogen, fluorine
for pink and white for hydrogen) and chitosan (grey:Glucosamine and
yellow: N-acyl-glucosamine) for (a) acetyl modified chitosan with DOX,
(b) butanoyl modified chitosans with DOX, (c) acetyl modified chitosan
with GEM, and (d) butanoyl modified chitosans with GEM. ....... 81
5.2 Lennard-Jones and coulombic contribution to the DOX-modification group
interaction energy for (a) acetyl-chitosan and (b) butanoyl-chitosan; and
Lennard-Jones and coulombic contribution to the DOX-backbone interac-
tion energy for (c) acetyl-chitosan and (d) butanoyl-chitosans at different
χ......................................... 82
5.3 Lennard-Jones and coulombic contribution to the GEM-modification group
interaction energy for (a) acetyl-chitosan and (b) butanoyl-chitosan; and
Lennard-Jones and coulombic contribution to the GEM-backbone interac-
tion energy for (c) acetyl-chitosan and (d) butanoyl-chitosans at different
χ......................................... 83
5.4 Hydrogen bond contacts for between DOX and modification groups in
a (a) acetyl- modified and (b) butanoyl-modified; and hydrogen bond
contacts for between DOX and the backbone in a (c) acetyl-modified and
(d) butanoyl-modified chitosans at different degrees of modification. . . . 84
5.5 Drug diffusion constants vs. χfor single-drug migration across different
chitosan networks for evenly-spaced (black) and blocky (red) modification
patterns: (a) GEM and (b) DOX in acetyl-chitosan, and (c) GEM and
(d) DOX in butanoyl-chitosan. ........................ 87
106
5.6 Number and (percentage of total) of non-bonded interactions observed be-
tween drug molecules and chitosan chains (backbone and modifications)
in the different chitosan networks over the drug molecules trajectories
during the simulation; where DOX molecules are treated as a group so
that any given chitosan or modification site can only contribute one con-
tact. In the table, low represents χ= 16% for both systems, and high
represents χ= 50% and 32% for acetyl-and butanoyl- modified systems,
respectively .................................. 88
5.7 Mean-squared displacement (MSD) plot vs. time for diffusion of DOX
through the blocky butyl-modified chitosan network at χ= 32%, where
DOX-Captured refers to DOX that becomes entrapped within a clus-
ter and DOX-Free refers to DOX that remains in the pores of the clus-
ter/channel morphology during the simulation. ............... 89
5.8 Snapshot of DOX migration through: (a) acetyl-chitosan networks at low
χ(16%); (b) acetyl-chitosan networks at high χ(50%); (c) butanoyl-
chitosan network at low χ(16%); (d) butanoyl-chitosan networks at high
χ(32%) where the chitosan backbone is represented by red beads, modifi-
cations are represented by yellow beads, and DOX is represented by black
beads ...................................... 90
5.9 Snapshot of GEM migration through: (a) acetyl-chitosan networks at low
χ(16%); (b) acetyl-chitosan networks at high χ(50%); (c) butanoyl-
chitosan network at low χ(16%); (d) butanoyl-chitosan networks at high
χ(32%) where the chitosan backbone is represented by red beads, mod-
ifications are represented by yellow beads, and GEM is represented by
green beads .................................. 91
5.10 Mean-square Displacement of DOX (a) Acetylated chitosan network (b)
butylated chitosan network for a uniformly-spaced modification pattern. 93
5.11 Drug diffusion constants vs. χfor dual-drug migration across different chi-
tosan networks: (a) GEM in acetyl-chitosan, (b) DOX in acetyl-chitosan,
(c) GEM in butanoyl-chitosan, and (d) DOX in butanoyl-chitosan. . . . . 95
5.12 Radial distribution functions of GEM around the center of mass of DOX
in acetyl-chitosan networks with (a) evenly-spaced and (b) blocky modi-
fication, and butanoyl-chitosan networks with (c) evenly-spaced and (d)
blocky modification. ............................. 96
List of Tables
2.1 Comparison between atomistic and CG end-to-end distance and radius of
gyration for cellulose and chitosan and their derivatives. Here A and B
correspond to alternating and blockwise modification respectively . . . . 37
2.2 Polymer end-to-end distance for single polymers in solution. The errors
represent one standard deviation. Here A and B correspond to alternating
and blockwise modification respectively ................... 42
3.1 Table showing comparison between atomistic and CG end-to-end distance
and radius of gyration for (a). 16% acetylation ,(b). 16% Butylation and
(c). 16% Heptylation ............................. 56
4.1 Number of contacts between modification beads formed in the chitosan
networks ................................... 67
4.2 Polymer end-to-end distance for single polymers in solution and in the
network. The errors represent one standard deviation. ........... 67
4.3 Number of contacts between modification beads formed in the chitosan
networks with modifications grouped in blocks of four ........... 75
5.1 Diffusion coefficients (10−5cm2/s) comparison for DOX and GEM with
pure water in atomistic and CG simulation respectively. .......... 85
5.2 Diffusion constant (10−5cm2/s) for DOX and GEM for acetylation and
butylation at different degree of modification. ............... 93
5.3 Average number of contacts for DOX and GEM for acetylation and buty-
lation at different degree of modification. .................. 94
Bibliography
1Lois M Crowe, Robert Mouradian, John H Crowe, Susan A Jackson, and Christopher
Womersley. Effects of carbohydrates on membrane stability at low water activities.
Biochimica et Biophysica Acta (BBA)-Biomembranes, 769(1):141–150, 1984.
2R James Stubbs, AM Prentice, and WP James. Carbohydrates and energy balance.
Annals of the New York Academy of Sciences, 819:44–69, 1997.
3Max Kleiber et al. The fire of life. an introduction to animal energetics. The fire of
life. An introduction to animal energetics., 1961.
4Paul R Crocker and Ten Feizi. Carbohydrate recognition systems: functional triads
in cell—cell interactions. Current opinion in structural biology, 6(5):679–691, 1996.
5Akiko Takada, Katsuyuki Ohmori, Naofumi Takahashi, Kiyotaka Tsuyuoka, Akihiro
Yago, Koichi Zenita, Akira Hasegawa, and Reiji Kannagi. Adhesion of human can-
cer cells to vascular endothelium mediated by a carbohydrate antigen, sialyl lewis a.
Biochemical and biophysical research communications, 179(2):713–719, 1991.
6Klaus N Englyst and Hans N Englyst. Carbohydrate bioavailability. British Journal
of Nutrition, 94(1):1–11, 2005.
7Artur Cavaco-Paulo. Mechanism of cellulase action in textile processes. Carbohydrate
Polymers, 37(3):273–277, 1998.
8Zhijun Shi, Yue Zhang, Glyn O Phillips, and Guang Yang. Utilization of bacterial
cellulose in food. Food hydrocolloids, 35:539–545, 2014.
9Robert L Barcus and David W Bjorkquist. Fibers and pulps for papermaking based
on chemical combination of poly (acrylate-co-itaconate), polyol and cellulosic fiber,
November 3 1992. US Patent 5,160,789.
109
10 Hongli Zhu, Wei Luo, Peter N Ciesielski, Zhiqiang Fang, JY Zhu, Gunnar Henriksson,
Michael E Himmel, and Liangbing Hu. Wood-derived materials for green electronics,
biological devices, and energy applications. Chemical Reviews, 116(16):9305–9374,
2016.
11 Yang Yu, Theodore Tyrikos-Ergas, Yuntao Zhu, Giulio Fittolani, Vittorio Bordoni,
Ankush Singhal, Richard J Fair, Andrea Grafmüller, Peter H Seeberger, and Mar-
tina Delbianco. Systematic hydrogen-bond manipulations to establish polysaccharide
structure–property correlations. Angewandte Chemie, 131(37):13261–13266, 2019.
12 Stephen C Fry. The structure and functions of xyloglucan. Journal of Experimental
Botany, 40(1):1–11, 1989.
13 Hale Cigdem Arca, Laura I Mosquera-Giraldo, Vivian Bi, Daiqiang Xu, Lynne S Tay-
lor, and Kevin J Edgar. Pharmaceutical applications of cellulose ethers and cellulose
ether esters. Biomacromolecules, 19(7):2351–2376, 2018.
14 Per M Claesson and Barry W Ninham. ph-dependent interactions between adsorbed
chitosan layers. Langmuir, 8(5):1406–1412, 1992.
15 Levan Tsereteli and Andrea Grafmüller. An accurate coarse-grained model for chitosan
polysaccharides in aqueous solution. PloS one, 12(7):e0180938, 2017.
16 Claire Chatelet, Odile Damour, and Alain Domard. Influence of the degree of acety-
lation on some biological properties of chitosan films. Biomaterials, 22(3):261–268,
2001.
17 Steven W Benner, Vijay T John, and Carol K Hall. Simulation study of hydrophobi-
cally modified chitosan as an oil dispersant additive. The Journal of Physical Chem-
istry B, 119(23):6979–6990, 2015.
18 Yunqiang Chen, Libin Chen, Hua Bai, and Lei Li. Graphene oxide–chitosan composite
hydrogels as broad-spectrum adsorbents for water purification. Journal of Materials
Chemistry A, 1(6):1992–2001, 2013.
19 Hiroshi Ueno, Takashi Mori, and Toru Fujinaga. Topical formulations and wound
healing applications of chitosan. Advanced drug delivery reviews, 52(2):105–115, 2001.
20 Narayan Bhattarai, Jonathan Gunn, and Miqin Zhang. Chitosan-based hydrogels for
controlled, localized drug delivery. Advanced drug delivery reviews, 62(1):83–99, 2010.
110
21 John D Schneible, Ankush Singhal, Radina L Lilova, Carol K Hall, Andrea Grafmuller,
and Stefano Menegatti. Tailoring the chemical modification of chitosan hydrogels to
fine-tune the release of a synergistic combination of chemotherapeutics. Biomacro-
molecules, 20(8):3126–3141, 2019.
22 Todd R Hoare and Daniel S Kohane. Hydrogels in drug delivery: Progress and chal-
lenges. Polymer, 49(8):1993–2007, 2008.
23 Yong Qiu and Kinam Park. Environment-sensitive hydrogels for drug delivery. Ad-
vanced drug delivery reviews, 53(3):321–339, 2001.
24 Paul Tangney and Sandro Scandolo. An ab initio parametrized interatomic force field
for silica. The Journal of chemical physics, 117(19):8898–8904, 2002.
25 Luca Monticelli, Senthil K Kandasamy, Xavier Periole, Ronald G Larson, D Peter
Tieleman, and Siewert-Jan Marrink. The martini coarse-grained force field: extension
to proteins. Journal of chemical theory and computation, 4(5):819–834, 2008.
26 Karl N Kirschner, Austin B Yongye, Sarah M Tschampel, Jorge González-Outeiriño,
Charlisa R Daniels, B Lachele Foley, and Robert J Woods. Glycam06: a general-
izable biomolecular force field. carbohydrates. Journal of computational chemistry,
29(4):622–655, 2008.
27 Jörg Sauter and Andrea Grafmüller. Procedure for transferable coarse-grained models
of aqueous polysaccharides. Journal of chemical theory and computation, 13(1):223–
236, 2016.
28 Martin Lısal, Jiřı Kolafa, and Ivo Nezbeda. An examination of the five-site potential
(tip5p) for water. The Journal of chemical physics, 117(19):8892–8897, 2002.
29 Jorg Sauter and Andrea Grafmuller. Solution properties of hemicellulose polysaccha-
rides with four common carbohydrate force fields. Journal of chemical theory and
computation, 11(4):1765–1774, 2015.
30 Sergei Izvekov and Gregory A Voth. A multiscale coarse-graining method for biomolec-
ular systems. The Journal of Physical Chemistry B, 109(7):2469–2473, 2005.
31 Dirk Reith, Mathias Pütz, and Florian Müller-Plathe. Deriving effective mesoscale po-
tentials from atomistic simulations. Journal of computational chemistry, 24(13):1624–
1636, 2003.
111
32 William George Noid, Jhih-Wei Chu, Gary S Ayton, Vinod Krishna, Sergei Izvekov,
Gregory A Voth, Avisek Das, and Hans C Andersen. The multiscale coarse-graining
method. i. a rigorous bridge between atomistic and coarse-grained models. The Journal
of chemical physics, 128(24):244114, 2008.
33 Siewert J Marrink, Alex H De Vries, and Alan E Mark. Coarse grained model for
semiquantitative lipid simulations. The Journal of Physical Chemistry B, 108(2):750–
760, 2004.
34 Mark Bathe, Gregory C Rutledge, Alan J Grodzinsky, and Bruce Tidor. A coarse-
grained molecular model for glycosaminoglycans: application to chondroitin, chon-
droitin sulfate, and hyaluronic acid. Biophysical journal, 88(6):3870–3887, 2005.
35 Karim Mazeau, Serge Pérez, and Marguerite Rinaudo. Predicted influence of n-acetyl
group content on the conformational extension of chitin and chitosan chains. 2000.
36 Pu Liu, Sergei Izvekov, and Gregory A Voth. Multiscale coarse-graining of monosac-
charides. The Journal of Physical Chemistry B, 111(39):11566–11575, 2007.
37 Jörg Sauter and Andrea Grafmüller. Efficient osmotic pressure calculations using
coarse-grained molecular simulations. Journal of chemical theory and computation,
14(3):1171–1176, 2018.
38 Steven W Benner and Carol K Hall. Development of a coarse-grained model of chitosan
for predicting solution behavior. The Journal of Physical Chemistry B, 120(29):7253–
7264, 2016.
39 Peter H Seeberger. The logic of automated glycan assembly. Accounts of chemical
research, 48(5):1450–1463, 2015.
40 Oliviana Calin, Steffen Eller, and Peter H Seeberger. Automated polysaccharide syn-
thesis: assembly of a 30mer mannoside. Angewandte Chemie International Edition,
52(22):5862–5865, 2013.
41 MJ Abraham, D van der Spoel, E Lindahl, and B Hess. Gromacs user manual, version
5.1. 2; gromacs development team, 2016. There is no corresponding record for this
reference.
42 DA Case, TA Darden, TE Cheatham III, CL Simmerling, J Wang, RE Duke, R Luo,
RC Walker, W Zhang, KM Merz, et al. Amber 12; university of california: San
112
francisco, 2012. There is no corresponding record for this reference, pages 1–826,
2010.
43 Eric J Sorin and Vijay S Pande. Exploring the helix-coil transition via all-atom
equilibrium ensemble simulations. Biophysical journal, 88(4):2472–2493, 2005.
44 Marko Wehle, Ivan Vilotijevic, Reinhard Lipowsky, Peter H Seeberger, Daniel
Varon Silva, and Mark Santer. Mechanical compressibility of the glycosylphos-
phatidylinositol (gpi) anchor backbone governed by independent glycosidic linkages.
Journal of the American Chemical Society, 134(46):18964–18972, 2012.
45 Sander Pronk, Szilárd Páll, Roland Schulz, Per Larsson, Pär Bjelkmar, Rossen Apos-
tolov, Michael R Shirts, Jeremy C Smith, Peter M Kasson, David Van Der Spoel, et al.
Gromacs 4.5: a high-throughput and highly parallel open source molecular simulation
toolkit. Bioinformatics, 29(7):845–854, 2013.
46 Junmei Wang, Romain M Wolf, James W Caldwell, Peter A Kollman, and David A
Case. Development and testing of a general amber force field. Journal of computational
chemistry, 25(9):1157–1174, 2004.
47 Christopher I Bayly, Piotr Cieplak, Wendy Cornell, and Peter A Kollman. A well-
behaved electrostatic potential based method using charge restraints for deriving
atomic charges: the resp model. The Journal of Physical Chemistry, 97(40):10269–
10280, 1993.
48 Tom Darden, Darrin York, and Lee Pedersen. Particle mesh ewald: An n log
(n) method for ewald sums in large systems. The Journal of chemical physics,
98(12):10089–10092, 1993.
49 Berk Hess, Henk Bekker, Herman JC Berendsen, and Johannes GEM Fraaije. Lincs: a
linear constraint solver for molecular simulations. Journal of computational chemistry,
18(12):1463–1472, 1997.
50 Shuichi Miyamoto and Peter A Kollman. Settle: An analytical version of the shake
and rattle algorithm for rigid water models. Journal of computational chemistry,
13(8):952–962, 1992.
51 Shuichi Nosé. S. nosé, j. chem. phys. 81, 511 (1984). J. Chem. Phys., 81:511, 1984.
113
52 William G Hoover. Canonical dynamics: Equilibrium phase-space distributions. Phys-
ical review A, 31(3):1695, 1985.
53 Michele Parrinello and Aneesur Rahman. Polymorphic transitions in single crystals:
A new molecular dynamics method. Journal of Applied physics, 52(12):7182–7190,
1981.
54 Shuichi Nosé and ML Klein. Constant pressure molecular dynamics for molecular
systems. Molecular Physics, 50(5):1055–1076, 1983.
55 Ilario G Tironi, René Sperb, Paul E Smith, and Wilfred F van Gunsteren. A gen-
eralized reaction field method for molecular dynamics simulations. The Journal of
chemical physics, 102(13):5451–5459, 1995.
56 Wilfred F Van Gunsteren and Herman JC Berendsen. A leap-frog algorithm for
stochastic dynamics. Molecular Simulation, 1(3):173–185, 1988.
57 Supriyo Bhattacharya and Keith E Gubbins. Fast method for computing pore size
distributions of model materials. Langmuir, 22(18):7726–7731, 2006.
58 Franz Kappel and Alexei V Kuntsevich. An implementation of shor’s r-algorithm.
Computational Optimization and Applications, 15(2):193–205, 2000.
59 Victor Ruhle, Christoph Junghans, Alexander Lukyanov, Kurt Kremer, and Denis
Andrienko. Versatile object-oriented toolkit for coarse-graining applications. Journal
of chemical theory and computation, 5(12):3211–3223, 2009.
60 SY Mashayak, Mara N Jochum, Konstantin Koschke, Narayana R Aluru, Victor Rühle,
and Christoph Junghans. Relative entropy and optimization-driven coarse-graining
methods in votca. PLoS one, 10(7):e0131754, 2015.
61 Gregg T Beckham and Michael F Crowley. Examination of the α-chitin structure and
decrystallization thermodynamics at the nanoscale. The Journal of Physical Chemistry
B, 115(15):4516–4522, 2011.
62 Eduardo F Franca, Roberto D Lins, Luiz CG Freitas, and TP Straatsma. Charac-
terization of chitin and chitosan molecular structure in aqueous solution. Journal of
Chemical Theory and Computation, 4(12):2141–2149, 2008.
114
63 Wesley K Lay, Mark S Miller, and Adrian H Elcock. Optimizing solute–solute interac-
tions in the glycam06 and charmm36 carbohydrate force fields using osmotic pressure
measurements. Journal of chemical theory and computation, 12(4):1401–1407, 2016.
64 Jodi A Hadden, Alfred D French, and Robert J Woods. Unraveling cellulose microfib-
rils: a twisted tale. Biopolymers, 99(10):746–756, 2013.
65 Liming Hu, Yun Sun, and Yan Wu. Advances in chitosan-based drug delivery vehicles.
Nanoscale, 5(8):3103–3111, 2013.
66 Raphaël Riva, Heloise Ragelle, Anne des Rieux, Nicolas Duhem, Christine Jérôme,
and Véronique Préat. Chitosan and chitosan derivatives in drug delivery and tissue
engineering. In Chitosan for biomaterials II, pages 19–44. Springer, 2011.
67 Mokhamad Nur and Todor Vasiljevic. Can natural polymers assist in delivering insulin
orally? International journal of biological macromolecules, 103:889–901, 2017.
68 Seong-Chul Hong, Seung-Yup Yoo, Hyeongmin Kim, and Jaehwi Lee. Chitosan-based
multifunctional platforms for local delivery of therapeutics. Marine drugs, 15(3):60,
2017.
69 PL Lam and Roberto Gambari. Advanced progress of microencapsulation technologies:
In vivo and in vitro models for studying oral and transdermal drug deliveries. Journal
of controlled release, 178:25–45, 2014.
70 Jorg Sauter and Andrea Grafmuller. Predicting the chemical potential and osmotic
pressure of polysaccharide solutions by molecular simulations. Journal of chemical
theory and computation, 12(9):4375–4384, 2016.
71 Ornella Ortona, Gerardino D’Errico, Gaetano Mangiapia, and Donato Ciccarelli. The
aggregative behavior of hydrophobically modified chitosans with high substitution
degree in aqueous solution. Carbohydrate polymers, 74(1):16–22, 2008.
72 Gang-Biao Jiang, Daping Quan, Kairong Liao, and Haihua Wang. Preparation of
polymeric micelles based on chitosan bearing a small amount of highly hydrophobic
groups. Carbohydrate Polymers, 66(4):514–520, 2006.
73 M Rinaudo, R Auzely, C Vallin, and I Mullagaliev. Specific interactions in modified
chitosan systems. Biomacromolecules, 6(5):2396–2407, 2005.
115
74 Sei-ichi Aiba. Studies on chitosan: 3. evidence for the presence of random and block
copolymer structures in partially n-acetylated chitosans. International journal of bio-
logical macromolecules, 13(1):40–44, 1991.
75 Keisuke Kurita, Takanori Sannan, and Yoshio Iwakura. Studies on chitin, 4. evi-
dence for formation of block and random copolymers of n-acetyl-d-glucosamine and
d-glucosamine by hetero-and homogeneous hydrolyses. Die Makromolekulare Chemie:
Macromolecular Chemistry and Physics, 178(12):3197–3202, 1977.
76 Mette H Ottey, Kjell M Vårum, and Olav Smidsrød. Compositional heterogeneity of
heterogeneously deacetylated chitosans. Carbohydrate Polymers, 29(1):17–24, 1996.
77 Jérôme Berger, Marianne Reist, Joachim M Mayer, Olivia Felt, NA Peppas, and
Robert Gurny. Structure and interactions in covalently and ionically crosslinked chi-
tosan hydrogels for biomedical applications. European Journal of Pharmaceutics and
Biopharmaceutics, 57(1):19–34, 2004.
78 Simon M Iveson, James D Litster, Karen Hapgood, and Bryan J Ennis. Nucleation,
growth and breakage phenomena in agitated wet granulation processes: a review.
Powder technology, 117(1-2):3–39, 2001.
79 Min Mu, Xiaoling Li, Aiping Tong, and Gang Guo. Multi-functional chitosan-based
smart hydrogels mediated biomedical application. Expert opinion on drug delivery,
16(3):239–250, 2019.
80 Michelly CG Pellá, Michele K Lima-Tenório, Ernandes T Tenório-Neto, Marcos R
Guilherme, Edvani C Muniz, and Adley F Rubira. Chitosan-based hydrogels: From
preparation to biomedical applications. Carbohydrate Polymers, 196:233–245, 2018.
81 Ecaterina Stela Dragan and Maria Valentina Dinu. Polysaccharides constructed hy-
drogels as vehicles for proteins and peptides. a review. Carbohydrate polymers, page
115210, 2019.
82 Michael W Mahoney and William L Jorgensen. Diffusion constant of the tip5p model
of liquid water. The Journal of Chemical Physics, 114(1):363–366, 2001.
83 Marco Biondi, Sabato Fusco, Andrew L Lewis, and Paolo A Netti. New insights into
the mechanisms of the interactions between doxorubicin and the ion-exchange hydrogel
116
dc beadTM for use in transarterial chemoembolization (tace). Journal of Biomaterials
Science, Polymer Edition, 23(1-4):333–354, 2012.
84 Cary T Chiou, David W Schmedding, and Milton Manes. Partitioning of organic
compounds in octanol-water systems. Environmental Science & Technology, 16(1):4–
10, 1982.
85 Suvakanta Dash, Padala Narasimha Murthy, Lilakanta Nath, Prasanta Chowdhury,
et al. Kinetic modeling on drug release from controlled drug delivery systems. Acta
Pol Pharm, 67(3):217–23, 2010.
86 T Higuchi. Mechanism of sustained-action medication. theoretical analysis of rate of
release of solid drugs dispersed in solid matrices. Journal of pharmaceutical sciences,
52(12):1145–1149, 1963.
87 Kathryn M Camacho, Stefano Menegatti, Douglas R Vogus, Anusha Pusuluri, Zoë
Fuchs, Maria Jarvis, Michael Zakrewsky, Michael A Evans, Renwei Chen, and Samir
Mitragotri. Dafodil: a novel liposome-encapsulated synergistic combination of doxoru-
bicin and 5fu for low dose chemotherapy. Journal of Controlled Release, 229:154–162,
2016.
88 Dechun Liu, Yanbin Chen, Xiaoshan Feng, Miao Deng, Gangqiang Xie, Jianguang
Wang, Like Zhang, Qipeng Liu, and Pengfei Yuan. Micellar nanoparticles loaded
with gemcitabine and doxorubicin showed synergistic effect. Colloids and Surfaces B:
Biointerfaces, 113:158–168, 2014.
89 Twan Lammers, Vladimir Subr, Karel Ulbrich, Peter Peschke, Peter E Huber, Wim E
Hennink, and Gert Storm. Simultaneous delivery of doxorubicin and gemcitabine to
tumors in vivo using prototypic polymeric drug carriers. Biomaterials, 30(20):3466–
3475, 2009.
90 Shuqin Ni, Lei Qiu, Guodong Zhang, Haiyan Zhou, and Yong Han. Lymph cancer
chemotherapy: delivery of doxorubicin–gemcitabine prodrug and vincristine by nanos-
tructured lipid carriers. International journal of nanomedicine, 12:1565, 2017.
91 Sonke Svenson, Marc Wolfgang, Jungyeon Hwang, John Ryan, and Scott Eliasof. Pre-
clinical to clinical development of the novel camptothecin nanopharmaceutical crlx101.
Journal of controlled release, 153(1):49–55, 2011.
117
92 Thomas Schluep, Jungyeon Hwang, Isabel J Hildebrandt, Johannes Czernin, Chung
Hang J Choi, Christopher A Alabi, Brendan C Mack, and Mark E Davis. Phar-
macokinetics and tumor dynamics of the nanoparticle it-101 from pet imaging and
tumor histological measurements. Proceedings of the National Academy of Sciences,
106(27):11394–11399, 2009.
93 Glen J Weiss, Joseph Chao, Jeffrey D Neidhart, Ramesh K Ramanathan, Dawn Bas-
sett, James A Neidhart, Chung Hang J Choi, Warren Chow, Vincent Chung, Stephen J
Forman, et al. First-in-human phase 1/2a trial of crlx101, a cyclodextrin-containing
polymer-camptothecin nanopharmaceutical in patients with advanced solid tumor ma-
lignancies. Investigational new drugs, 31(4):986–1000, 2013.
94 Andrew J Clark, Devin T Wiley, Jonathan E Zuckerman, Paul Webster, Joseph Chao,
James Lin, Yun Yen, and Mark E Davis. Crlx101 nanoparticles localize in human
tumors and not in adjacent, nonneoplastic tissue after intravenous dosing. Proceedings
of the National Academy of Sciences, 113(14):3850–3854, 2016.
118
