Molecular dynamics studies of amylose plasticized with Brazilian Cerrado oils : part I

1Laboratório de Estudos Estruturais Moleculares – LEEM, Instituto de Química – IQ, Universidade de Brasília – UnB, Campus Darcy Ribeiro, Brasília, DF, Brasil 2Laboratório de Pesquisa em Polímeros e Nanomateriais – LabPolN, Instituto de Química – IQ, Universidade de Brasília – UnB, Campus Darcy Ribeiro, Brasília, DF, Brasil 3Laboratoire d’Ingénierie des Biomolécules – LIBio, Ecole Nationale Supérieure d’Agronomie et des Industries Alimentaires – ENSAIA, Institut National Polytechnique de Lorraine – INPL, Université de Lorraine – UL, Vandœuvre-lès-Nancy, France


Introduction
The use of starch-based polymers in films and wrapping, bags, personal hygiene items, cups, straws, strings, and many other materials is recognized for being a cost-effective green technology in comparison to synthetic polymeric materials derived from petroleum [1] .However, all these products demand for starch with some level of chemical manipulation, because in its natural form starch presents low film processability and poor mechanical properties [2][3][4][5] .
Starch is a semi-crystalline biopolymer [5] mainly formed by two components with repeating α-D-glucopyranosyl units, although other components such as proteins and lipids are present in small ammouts.These two main components are amylose, a mostly linear polysaccharide and amylopectin, a highly branched one.In the amylose linear main chain, the repeating units are connected by α(1-4)-glycosidic bonds (Figure 1), while in longer chains some branches may be present.It has a tendency to acquire an stable helical shape when complexed with several ligands [5] .The quantity of amylose in a starch can make it form stronger films, the amount of amylose in normal starches (i.e.maize, potato, rice, etc) is about 20-30% [5,6] .
A common approach for improving starch film ability is plasticization, which can be efficiently and environmentally friendly achieved with small amounts of vegetable oils [9] and its fatty acids [2,7,[10][11][12] .
Our group has developed starch films produced by casting of aqueous starch suspensions plasticized with vegetable oils from Cerrado, the second largest biome in South America [13] , as well as from Amazonia.More specifically, with pequi (Caryocar brasiliense) and buriti (Mauritia flexuosa L.) fruits, which are used in the region with different purposes such as human and animal feeding [14,15] , cosmetic and medicinal uses [16] .The fundamental importance in its composition is the abundance in unsaturated fatty acids, with the predominance of oleic acid, between 50-80% in mass, and also, the presence of tocopherols and carothenoids [11,14,15] .They have proved to efficiently plasticize starch and enable fabrication of films by casting [2,3,7,11] .A step further consisted on the preparation of bionanocomposites of plasticized starch and lamellar oxides, such as montimorillonite clay, for improved thermal and barrier properties [1,2,7,8,10,12] .Therefore, plasticization is an essential step to allow the use of starch in advanced biodegradable materials [2,3,7,12] .
The plasticization of starch with vegetable oils opens-up a venue for new biomaterials and processing as well as for theoretical and molecular models of chemical interactions.In particular, molecular dynamics simulation (MD) has been performed to understand, at the molecular level, the interaction between starch components, such as amylose, and common solvents (water and DMSO) as well as to evaluate current and future candidates for possible plasticization.This approach should be advantageous before any expensive and time-consuming real processing is carried out.For example, Yang et al. [17] performed MD of a single amylose chain (36 monomeric units) combined with glycerol as plasticizer.They found that strong hydrogen bonding was essential for cohesion of the amylose/glycerol model, which was further corroborated by experimental data (X-ray diffraction, tensile strength and water uptake), probed with starch/glycerol cast films.Tusch et al. [18] observed that at certain simulation times the helix character of a single helical amylose chain (55 monomeric units) in water/DMSO mixtures was stabilized after increasing the DMSO concentration, since DMSO helps to preserve intramolecular hydrogen bonds in starch.Their findings were in consonance with NMR data.López et al. [19] verified by MD that the addition of lipids, such as dipalmitoyl-phosphatidylcholine and glycerol monooleate, induces amylose fragments (13 monomeric units) to folding into a helical structure.The latter matches the V-amylose structure, which is a type of amylose polymorph [5,20] , observed in X-ray diffraction.Finally, Feng et al. [21] observed by MD similar V-amylose structures when in the presence of linoleic acid.
As observed in experimental data and corroborated further by MD studies, inter-and intra-molecular interactions play a major role on the molecular behavior of complex systems comprised by natural plasticizers and polysaccharides [12,[17][18][19]21] . In fct, they are essential for the rational design of high performance biomaterials as well as to predict its properties.Most of MD studies were limited to a single and short poly-glucose chains, mainly because larger and more complex systems would demand for longer simulation times and even more expensive computational instrumentation [17][18][19]21] .
In this regard, the present contribution provides a systematic study performed by Molecular Mechanics (MM) and Molecular Dynamics methods with the Polymer Consistent Force Field (PCFF) on the assembly of multiple amylose oligomers (25 to 120 monomeric units) in aqueous medium.The oligomer is complexed with a mixture of fatty acids (oleic, palmitic, and stearic) set in proportions that mimic those found in pequi oil.The molecular assemblies were confined into cells and studied under periodic boundary condition (PBC) and constant number of particles, temperature and volume (NVT) ensemble, at 363 K, at variable elapsed times.Its main goal was directed toward observation of the molecules' motion when evolving together, to the bonded and non-bonded intermolecular forces, and to the respective structural and behavioral correlations.These changes are known for reflecting the strength of intra-and intermolecular interactions that act on the molecular system.

Materials and Methods
All calculations were performed with BIOVIA Materials Studio TM suite [22] .The atomistic coordinates for all molecules and simulated cells were built using the Visualizer, Builder and Amorphous Cell modules [22] .The MM and MD calculations were carried out with PCFF [23,24] interfaced to the Forcite module [22] .The PCFF is a member of the consistent force field family.It was developed with a base on CFF91 and is intended for application to polymers, organic materials, metals, and zeolites [22][23][24][25][26] , also being suitable for carbohydrates [3,27] and water solvation [28] .The dynamic trajectories were carefully analyzed using data correlation graphs implemented in Analysis/Forcite [22] .

Specificities for the calculation protocol of amylosefatty acid -water systems: studies by NVT ensemble
The atomistic simulations essentially followed the calculation protocols hereafter described.
(i) The molecular systems were composed for a better understanding about the intra-and intermolecular interactions and the structural behavior when the three chemical species evolve together.These simulations were carefully made with different systems, which varied, fundamentally, in the proportion of fatty acids; in the initial form and number of the amylose oligomers, in which the sizes varied from 25 to 120 monomers.These two chemical species were created from features implemented in Visualizer and Builder programs, such as Fragment browser.The amylose chains were built by "Build Polymers" feature as isotactic homopolymers from the carbohydrates database of the BIOVIA Materials StudioTM program, using 1,4-α-D-glucose repeat units, with chains of variable lengths.After each construction, geometries and energies were fully optimized by MM and the results stored for future usage.Periodic cells were created with the Build Crystals functionality.The unit cells dimensions varied according to the number and length of amylose chains.The fatty acids were manually inserted into the cells, near the amylose oligomers, in the proportions as close as possible to those found in the buriti and pequi oils.The unit cells were filled by water molecules ranging around 2000 to 4000 molecules in explicit solvent model, to simulate the experimental conditions carried out in previous studies [2,7,10] , where the components of the system were dispersed in distilled water.
(i.1) Intermediate energy minimizations allowed corrections and eliminated steric hindrance with the steepest descent algorithm.The number of iterations limited to coarse accuracy of the gradient was set to 2.0 x 10 -3 kcal mol -1 .Optimizations continued throughout the conjugate gradient algorithm, with the convergence of the gradient limited to 1.0 x 10 -3 kcal mol -1 Å -1 .In general terms, these limit values indicate medium precision.The molecular dynamics method of the Forcite program follows the classical Newtonian equations of motion for the integrated Verlet algorithm [29] .The initial injected speeds were random and assigned to the atoms according to the Maxwell-Boltzmann distribution.
The observation times ranged from 100 ps to 1 ns, according to the systems' sizes.For all simulations, the time step was set to 1 fs.The coordinates of the molecular systems corresponding to the dynamics fluctuations (frames) were saved in periodic intervals of at least 1000 frames, which allowed posterior analysis via control graphs and molecular coordinates, across the Analysis program features.All dynamic simulations ran under the canonical ensemble NVT.The temperature control in Forcite is based on the Nosé-Hoover chain (NHC) thermostat [30][31][32] .
All simulations were conducted under PBC.
After a manual and progressive addition of significant amounts of water, while being intermediated by some energy optimization steps, systems were submitted to fast dynamic trajectories at 500 K, during 40 ps, for their rearrangement.In this case, all other molecules had their atomic coordinates fixed under constraints.This procedure allowed for rapid movement of the solvent around other molecules and prevents water molecules, initially misplaced, to assume unrealistic positions and, consequently, to conduce for not productive results.Solvent molecules out of the unit cell were erased or manually repositioned.The procedure restarted until coarse system stabilization has been reached.
(i.2) After system optimization, constraints were removed and a new optimization was performed.Then, the system was submitted to a dynamic trajectory, at 363 K until reorganization.This is the temperature at which systems were experimentally studied. [2,7,10]imulation times depend mainly on the size of the system under observation.When equilibrium state was not critical, i.e., when systems were designed primarily to motion analysis and conformational modification, tolerances of gradient convergence (rms) were set to 1.0 x 10 -3 kcal mol -1 , for energy; 5.0 x 10 -1 kcal mol -1 Å -1 , for maximum strength; and the precision of the electrostatic and van der Waals terms to 1.0 x 10 -3 kcal mol -1 (medium accuracy).When the system reached situations in which the organization had already occurred, they were used as an input file for subsequent simulations; the parameters for the convergence tolerance were limited to 1.0 x 10 -4 kcal mol -1 for energy; 5.0 x10 -3 kcal mol -1 Å -1 , for force; 5.0 x 10 -5 Å, for displacement; and the precision of the nonbonding interactions terms, 1.0x10 -4 kcal mol -1 (fine precision).The nonbonding interactions of those periodic systems were computed by the Ewald summation method [33,34] .Under NVT ensemble, the automatic correction of the simulated cells dimensions was not used.
In addition to the general understanding about the behavior of these different chemical species working together, these models allowed one to verify the permissible size of the oligomeric chain and the number fatty acids and water molecules that could be added to the system.All together, they enabled the calculations in PBC, according to our computational capacity; a CPU time that would be logical and reliable to make the dynamic trajectories in different situations; and if the number and the length of the amylose chains participating to the systems could lead to different behavior.

First system modeled from four amylose short chains of 40 monomeric units and 49 fatty acid molecules
In the first simulated system, four amylose oligomers with 40 monomeric units each (842 atoms/oligomer; 3,368 atoms) were confined into an orthorhombic cell with dimensions of The fatty acids were placed perpendicular to the amylose chains' medial plane, in a proportion of 28 oleic acids, 20 palmitic acids and 1 stearic acid molecules (2,568 atoms).Finally, 2,668 molecules of water were gradually added that resulted in a system composed of 13,940 atoms in the asymmetric unit.
The starting (t = 0 s) and subsequent molecular assembly attained after different simulated elapsed times as well as the time dependence of dynamics energies (kinetic, potential, non-bond, and total energies) and temperature are depicted in Figure 2. As seen in the first image from Figure 2 (clockwise), the frontal view of the initial system comprises four elongated amylose chains lined two by two with fatty acids, while water molecules are randomly distributed around them.The subsequent image provides side view of the initial system.In analogy to two comets, heads and tails are composed of amylose and fatty acids, respectively, in which two are face to face while other two are distant of about 180 degrees.Then, solvent molecules were left moving around during a short dynamic process at 500 K, while the main molecules are kept fixed, under constraints.After water reorganization and energetic optimization, the constraints were removed and the free system was optimized again, before being submitted to a dynamic process at 363 K, for 200 ps (middle images).
The thermalization period took about 5 ps.At 10 ps, amylose chains remained elongated.Between 20 to 80 ps, amylose chains progressively coiled.The two "comets" (chains in lilac and light orange) that were farther from the central assemblies (chains in purple and dark orange) were strongly attracted and then distorted until they approached the core, leading to a new assembly in which the amylose chains were strongly aggregated (middle images).The nonpolar fatty acids tails were driven away from amylose while their polar heads remained near to each other, probably closer to the electronegative amylose groups via hydrogen bonds.Meanwhile, water molecules moved out of the cavities formed by amylose chains which, in turn, approached and coiled.The elongated fatty acids covered the amylose coil in an unusual rearrangement.This behavior can be monitored throughout frames from 40 to 200 ps.This simulation allowed us to observe that the key factor is not the stabilization time.The complete system folding was reached at 80 ps, and during the remaining 120 ps, only minor positions were changed.The last trajectory frame is subjected to geometric and energy correction (Figure 3), after repositioning of the few water molecules that may have been dispersed into the cell.The computation time (CPU) required to perform the dynamic path detailed above, at 363 K, in PBC, during 200 ps, was 88 h.This is quite reasonable to treat the 13,940 atoms that make up the system, in an 8 Gb RAM PC and 4 processors.

Second system modeled from one amylose chain of 120 monomeric units and 32 fatty acid molecules
In the second system, one amylose oligomeric chain containing 120 monomeric units (2,522 atoms) was built with the Amorphous Cell program.Fatty acid molecules were randomly distributed around the amylose chain, in a proportion of 18 oleic acids, 13 palmitic acids and 1 stearic acid molecules (1,678 atoms), in a total of 4,200 atoms.The original cell containing the oligosaccharide chain was enlarged to a cubic cell of a = b = c = 100 Å; α = β = γ = 90 o (Figure 4).Water molecules were gradually added, resulting in a system composed by 8,043 atoms into the cell.Energetic optimizations and dynamics trajectories followed protocols described above.Solvent molecules were progressively added and, after each significant quantity, the system was reorganized by fast molecular dynamics trajectories of 20 ps at 500 K, reoptimized to avoid steric hindrance, and the atomic coordinates of all organic molecules were placed under constraints.Water molecules dispersed into the cell were repositioned or erased, and another ensemble was added to the system.Then the constraints were deleted and dynamic trajectories by 30 ps, at 363 K were carried out until system stabilization was reached.It is important to note that during dynamics, some water molecules were always dispersed.Those water molecules were then deleted, the system was reoptimized, and the procedure was restarted, as frequently as necessary.
The analyses of those results are essentially the same, regardless the shape and length of the amylose chain.The starting radius of gyration and volume varied from 27.80 Å and 17,941.67Å 3 to 18.14 Å and 18,620.02Å 3 , respectively, at the end of the process.The amylose chain, progressively, became compacted until being coiled.The polar head of fatty acids stabilized closer to the electronegative amylose groups and formed a hydrogen bond network.Their nonpolar portions rejected water molecules and covered the oligomer "protecting" it.The water molecules covered every other surface of the oligomer with no fatty acids and stabilized by a very dense hydrogen bond network.The unexpected rearrangement comprised by elongated fatty acids encapsulating the amylose chains occurred again.Those results permit one to suppose that once the number of fatty acid molecules increases they will encapsulate the amylose chain and most of water molecules would be rejected.Only few water molecules would stay into the ensemble, due to their orientation close to oxygenated groups of the organic compounds.

Third system modeled from one amylose chain of 120 monomeric units and 96 fatty acid molecules
In order to verify those considerations, a third and last system among the many different systems studied will be described.A different conformer of amylose oligomer comprised by 120 monomeric units as built with the Amorphous Cells program was involved by randomly distributed fatty acid molecules.A higher concentration of fatty acids was considered, in a proportion of 54 oleic acid, 40 palmitic acid and 2 stearic acid molecules (5,028 atoms).Compounds were confined in a starting cell of a = 140 Å, b =120 Å, c = 100 Å; α = β = γ = 90 o (Figure 5).Since a fast optimization was accomplished, an amount of 2,207 water molecules were distributed into the asymmetric unit around the organic ensemble.The entire system had 13,342 atoms.Energetic optimizations and dynamics trajectories following the described protocols were carried out, firstly to correct water positions, leaving the atomic coordinates of organic molecules constrained.The misplaced water molecules were progressively deleted or correctly repositioned into the cell.The asymmetric unit axis were increased to a = 140 Å, b =140 Å, c = 130 Å and then, restrictions on the molecular freedom were withdrawn.The full procedure involved many gradient optimizations, intermediated by a few and fast dynamics trajectories at 363 K.When the system was considered stabilized enough and reached its final form, 1,406 water molecules were conserved inside it.The amylose radius of gyration and volume in the starting conformation were 32.37 Å and 17,945.86Å 3 , while the final conformation reached 20.52 Å and 18,326.74Å 3 , respectively.It was possible to confirm the precedent analysis and interpretations.Even randomly distributed, the majority of water molecules flow thru the organic system.The molecules farther from the central system forces stabilized among themselves.As the acids cover the most part of amylose surface, a remaining small portion could be covered by water.The electrostatic force of water was not strong enough to overcome the repulsion force imposed by nonpolar tails of fatty acids.A significant number of water molecules remain into the cell since they were stabilized by hydrogen bonds (Figure 5).However, if dynamics trajectories were longer, these water clusters would be progressively driven away from the organic assembly and, then, out of the unit cell.Among the 1,406 water molecules inside the cell, as shown in Figure 5, only 510 were estimated to be located around amylose and fatty acids polar heads to contribute for the system stabilization.Also, it should be noted that the gradient convergence is limited because of the water dispersion.The system converges following medium precision, because the water movement is slow but continuous.In spite of the strong water repulsion by the amylose-fatty acid complex, this solvent is absolutely essential to have a consistent system.
As fatty acids easily cover the polysaccharide chain we could interpret the behavior of those promising compounds when acting as plasticizers, which was already stablished based in results achieved by our research group [7,11] .

Conclusions
A systematic study was performed by Molecular Mechanics and Molecular Dynamics methods with PCFF on the assembly of multiple amylose oligomers (25 to 120 monomeric units) in aqueous medium while complexed with a mixture of fatty acids (oleic, palmitic, and stearic) set in proportions that mimic those found in pequi oil.The systems were assembled in simulated boxes, composed by helix amylose chains, lined with fatty acids and immersed in water molecules randomly distributed.All these calculations series were carried in NVT ensemble and PBC conditions.For all of configuration simulated, the behavior was essentially convergent.
As the simulation time elapsed, the non-polar tails of fatty acid molecules encapsulated amylose chains while most of water molecules were driven away from the amylose-fatty acid complex.Nonetheless, remaining water molecules were located at amylose sites uncoated by fatty acids and formed a dense network of hydrogen bonds, which was responsible for the entire system stabilization.The rapid and strong association of fatty to amylose chains found in the simulations suggests the behavior of those plasticizers with amylose.These theoretical results are in good agreement with what was experimentally observed for the interaction between fatty acids and starch.For this reason, we think that the molecular dynamics method presented in this work is useful for predicting the behavior of the interaction between fatty acids and more complex polysaccharides, more similar to starch, as well as systems enriched by inclusion of clay-minerals.

Acknowledgements
Felipe Silva thanks the Brazilian agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), for the doctoral fellowship in Brazil, and the

Figure 2 .
Figure 2. Two different perspectives of the starting system, in which fatty acids interact with four amylose oligomers, surrounded by water.The system had not yet been led to equilibrium (0 ps).The top and side views of the cell (left and right) show the spatial arrangement of key molecules, such as four "heads and tails of comets".(Middle) Snapshots of dynamics trajectory identifying fluctuations at 40, 74 and 200 ps.The oligomer chains are highlighted by spheres (CPK) and colored in lilac and purple, light and dark orange.Fatty acids are highlighted by thick lines (stick) and water molecules are visible because of the red color of the oxygen atoms.(Bottom) Graphics of dynamic process control (Energy (kcal mol -1 ) vs. time (ps) and temperature (K) vs. time (ps)) to a 200 ps trajectory at 363 K.

Figure 3 .
Figure 3. (Top) Zoom of the system segment already energetically optimized.It indicates the elongated fatty acids (in CPK) and water molecules (in stick) around amylose chains.(Bottom) View of segment highlighting the water molecules (in ball-and-stick) involving densely amylose oligomers (in CPK) in the regions without interactions with the amylose oligomers.

Figure 4 .
Figure 4. (Top) Initial system structure submitted for dynamics.Waters molecules were undisplayed for clarity.On the left, atoms are represented in stick and colored by elements colors -C, in grey; O, in red; H, in white.On the center, the similar image highlight compounds in CPK and different colors (amylose, in light orange; oleic, in medium blue; palmitic, in green, stearic, in light blue).(Top, right; down, left) Snapshots of dynamics trajectory, intermediate frame corresponding to 40 ps; last frame, 100 ps, when system was already coiled and water molecules were fully stabilized.Image shows the energetic optimized system to highlight the positions of fatty acids covering amylose chain and the solvent molecules distribution.Oxygen atoms of water molecules are colored in red.(Down, right) Snapshots of control graphics for Energy (kcal mol -1 ) vs. time (ps) and temperature (K) vs. time (ps).

Figure 5 .
Figure 5. (Top) Initial system structure submitted for dynamics Starting form and molecular distribution for the last model which is composed by 54 oleic acids; 40 palmitic acids and 2 stearic acids disposed around 120 monomers of amylose chain.Waters molecules were undisplayed for clarity.Atoms are represented in stick and colored by elements colors -C, in grey; O, in red; H, in white.In the center, the similar image highlight compounds in CPK and different colors (amylose, in purple; oleic, in medium blue; palmitic, in green, stearic, in light blue).(Down) -Final result's snapshots of the third model.Fatty acids cover amylose chain and there is no flown of water to amylose.Only few solvent molecules stay stabilized over amylose surface, when that surface was enough exposed.When concentration of fatty acid is dense, water molecules are strongly reject far from organic molecules.