Cellulose Nanofibrils and Mechanism of their ... - ACS Publications

Dec 14, 2015 - 31, St. Petersburg, 199004 Russia. •S Supporting Information. ABSTRACT: Molecular dynamics (MD) simulation of a nanofibril of native ...
0 downloads 0 Views 811KB Size
Subscriber access provided by GAZI UNIV

Article

Cellulose Nanofibrils and Mechanism of their Mineralization in Biomimetic Synthesis of Hydroxyapatite/Native Bacterial Cellulose Nanocomposites. Molecular Dynamic Simulations. Natalia Viacheslavovna Lukasheva, and Dmitry Alekseevich Tolmachev Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.5b03953 • Publication Date (Web): 14 Dec 2015 Downloaded from http://pubs.acs.org on December 21, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Cellulose Nanofibrils and Mechanism of their Mineralization in Biomimetic Synthesis of Hydroxyapatite/Native Bacterial Cellulose Nanocomposites. Molecular Dynamics Simulations.

N.V. Lukasheva*, D.A. Tolmachev Institute of Macromolecular Compounds, Russian Academy of Sciences, Bol’shoi pr. 31, St. Petersburg, 199004 Russia *e_mail: [email protected]

Abstract Molecular dynamics (MD) simulation of a nanofibril of native bacterial cellulose (BC) in solutions of mineral ions is presented. The supersaturated calcium-phosphate (CP) solution with the ionic composition of hydroxyapatite and CaCl2 solutions with the concentrations below, equal, and above the solubility limits are simulated. The influence of solvation models (TIP3P and TIP4P-ew water models) on structural characteristics of the simulated nanofibril and on the crystal nucleation process is assessed. The structural characteristics of cellulose nanofibrils (in particular, of the surface layer) are found to be nearly independent of the solvation models used in the simulation and on the presence of ions in the solutions. It is shown that ionic clusters are formed in the solution rather than on the fibril surface. The cluster sizes are slightly different for the two water models. The effect of the ion-ion interaction parameters on the results is discussed. The main conclusion is that the activity of hydroxyl groups on the BC fibril surface is not high enough to cause adsorption of Ca2+ ions from the solution. Therefore, the nucleation of CP crystals takes place initially in solution, and then the crystallites formed can be adsorbed on BC nanofibril surfaces.

ACS Paragon Plus Environment

1

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 29

I. INTRODUCTION Hybrid organic-mineral composites are considered to be the most promising materials for bone implants [1]. The composites based on bacterial cellulose (BC) and calcium phosphates (CP) [2-4] represent a new class of biomaterials designed for tissue engineering and prosthetic surgery. These materials successfully combine the properties of both components. The choice of BC as an organic matrix was dictated not only by its extremely high biocompatibility and chemical purity, but also by its remarkable physical properties, such as a high tensile strength and a high elastic modulus which are comparable to those of collagen in bone tissue [5, 6, 7]. Moreover, BC fibers are morphologically similar to collagen fibers because, like bone collagen, they have a hierarchical supermolecular structure [8]. BC possesses unique sorption properties due to its network structure consisting of ribbons with nanoscale thicknesses. Every ribbon is formed by a large number (10−100) of crystalline nanofibrils separated by nanochannels, and all surfaces are covered with hydrophilic OH groups. Synthetic CPs are analogs of the mineral component of bone tissue [1]. To obtain composite materials with the properties similar to those of bones, they are often produced by mimicking the natural mineralization process (biomimetic route), i.e., by soaking an organic matrix in the solution containing mineral ions. The biomimetic mineralization process is regarded as a heterogeneous crystallization, that is, mineral ions from the solution precipitate on organic matrix molecules, thereby forming crystallites. To provide heterogeneous crystallization, the matrix surface must contain the functional groups that can cause heterogeneous nucleation. To prepare BC- and CP-based composite materials, both native and chemically modified BC is used [4, 9-13]. If the organic matrix surface contains active groups capable of ionizing in an aqueous solution, thus acquiring a positive or negative charge (as in the case of chemically modified BC), the mineral layer formation on such a surface is mainly governed by electrostatic interactions.

ACS Paragon Plus Environment

2

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

The use of BC as a polymer matrix in the native form is extremely attractive first of all because of its chemical purity. The composition and structure of synthesized (native BC)/(CP) composites were investigated by different methods, such as X-ray photoelectron spectroscopy [4], energy-dispersive X-ray analysis [9], atomic emission spectroscopy [10, 12], infrared spectroscopy [9-13], electron microscopy [9-13], X-ray diffraction [9, 10, 13], and electron diffraction [12]. All studies revealed the formation of apatite crystals on BC nanofibrils. Moreover, the IR spectroscopy [9, 13] showed the presence of the red shift (towards lower wave numbers) of the band attributed to the intramolecular hydrogen bonds. Such a shift points to a strong interaction between the OH groups (which cover the nanofibril surface of native BC) and apatite. Some authors believe that minerals are formed on cellulose fibrils due to adsorption of mineral ions from the solution onto their surfaces [9, 10, 12, 13]. Their opinion is based on the assumption [14, 15] that polar functional groups, such as hydroxyl groups in cellulose, trigger the apatite precipitation. According to [14, 15], calcium binds to cellulose more firmly than water molecules and calcium cations may form weak bonds with negative dipole moments produced by oxygen atoms in neighboring hydroxyl groups. Since the bonds are weak, calcium cations can form bonds with phosphate anions. However, other authors [11] think that cellulose hydroxyl groups do not have a high enough reactivity to grow CP crystals. They suppose that though oxygen atoms of hydroxyl groups have negative partial atomic charges, the interaction between them and calcium ions (ion−polar interactions) are very weak in water solutions and probably cannot compensate for entropy losses. The problem is that experimental studies deal with the final stage, when the mineral phase is already formed. In order to study the mechanism of the CP crystal nucleation as the process which is regulated by native BC and to answer the question whether this process is a heterogeneous or homogeneous crystallization in its nature, computer simulations are necessary. To accelerate the mineralization process, it is necessary to use an increased saturation of the calcium phosphate solution surrounding the organic matrix or matrix pre-incubation in

ACS Paragon Plus Environment

3

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

calcium-containing solutions (typically CaCl2) [4, 9-13]. The rationale for the pre-incubation procedure is the supposition that hydroxyl groups of cellulose molecules having a strong negative dipole can chelate free cations in the CaCl2 solution and form coordinated bonds, thereby creating centers for subsequent heterogeneous crystallization of CP minerals. Our goal was to gain insight into initial stages of the mineralization process. To this end, molecular dynamics simulations of a native cellulose fibril placed into ionic solutions, such as a supersaturated solution with the ionic composition of hydroxyapatite (HAP) and the CaCl2 solutions with different concentrations were carried out. Our task was to examine the interaction of ions with the cellulose nanofibril surface and also the interaction of ions with each other and to evaluate the competing processes occurring on the fibril surface and in the solutions. By using the TIP3P and TIP4P-ew water models, the influence of the solvation models on structural characteristics of the simulated nanofibril and the crystal nucleation process was investigated. This article is organized as follows. Section II describes simulation models of the BC nanofibril (IIA), ionic solutions (IIB), and force fields (IIC). The methods, software and simulation procedures are considered in Section IID. The results of simulations for the nanofibril in calcium-phosphate and CaCl2 solutions are given in Section III. Section IV is Discussion. The evaluation of the BC induction ability relying on the classical theory of heterogeneous crystallization and its comparison with the simulation results are considered in Section IVA. The effect of the solvation models and also the impact of the parameterization of the van der Waals interactions for ions are discussed in Section IVB. Conclusions are given in Section V. Some explanatory and additional information is given in the Supporting Information.

II. SIMULATION MODELS, POTENTIALS, METHOD, SOFTWARE and PROCEDURE A. Model of BC nanofibril The BC crystalline structure was determined as the polymorphic modification Iβ possessing a monoclinic crystalline cell and P21 symmetry. It is built via an alternation of two-

ACS Paragon Plus Environment

4

Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

dimensional layers of glucose residues bonded by hydrogen bonds (Figure 1a), and the alternating layers are shifted with respect to each other along the molecular axis (с axis) (Fig. 1b) by ±c/4 [16]. To build a full atom model of the BC nanofibril structure, crystallographic data on the cellulose Iβ with unit cell parameters a = 0.82 nm, b =0.78 nm, c = 1.038 nm, β = 90°, α = 90°, and γ = 96.6° and also the data on atomic coordinates of cellulose molecules in the crystal cell were used [16].

a) b) Figure 1. (a) Initial structure of the BC nanofibril (cross section) and (b) alternating layers of cellulose molecules.

The fibril model consisted of 64 cellulose oligomer molecules, each containing 8 glucose residues. The nanofibril was placed into a computational cell with periodic boundary conditions in all directions. The periodic boundary conditions in the direction of the long axes of the molecules were used in such a way that the cellulose chains were covalently bonded through the periodic cell border. This made it possible to simulate infinitely long molecules. The cell dimensions were 4.436 nm along the cellulose molecules and 10.298 nm in two transverse directions. To model the fibril behavior in an aqueous solution of mineral ions, the simulation box with the cellulose fibril placed at its center was filled with water molecules and with ions at a required concentration. To create this concentration of the solution, a fraction of randomly selected water molecules was replaced by ions.

ACS Paragon Plus Environment

5

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

B. Models of Ionic Solutions i. Solution with the HAP ionic composition To investigate the processes of formation and growth of apatite clusters on the native BC fibril, the ion concentration in our simulation system was taken to be equal to that used in the study of hydroxyapatite nucleation at a collagen template [17]. This concentration (C=2.4 mol/kg) is higher than the saturation limit but a typical way to accelerate the crystal nucleation and growth in simulation is to increase the ion concentration. The quantities of ions and water molecules corresponding to the simulated concentration are presented in Table S1 in Supporting Information.

ii. CaCl2 solutions To investigate the process of adsorption of calcium ions on the cellulose fibril from the CaCl2 solution, three concentrations were considered: below the solubility limit, С1 (1.4 mol/kg); near the solubility limit, С2 (6.8 mol/kg); and above the solubility limit, С3 (11.5 mol/kg). The quantities of ions and water molecules corresponding to these concentrations are given in Table S2 in Supporting Information.

C. Interaction Potentials i. For BC The interaction potentials CSFF [18] used in this study for cellulose modeling belong to the CHARMM set of the force fields widely used in the biopolymeric system modeling. For CSFF, the parameters obtained previously for sugars were mainly used [19]; however, the parameters for the dihedral angles of the primary and secondary hydroxyl groups were modified. The obtained force field made it possible to attain a correspondence of the calculated conformational energies and rotational frequencies of hydroxyl groups with the experimental data on the behavior of sugars in solution. It is known [20-28] that the force fields that correctly

ACS Paragon Plus Environment

6

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

simulate cellulose characteristics in an aqueous environment cannot simulate its crystalline structure accurately enough. However, the main task of our study was to investigate adsorption of ions on the nanofibril surface. To this end, it is necessary to correctly simulate the dynamic characteristics of hydroxyl groups, but we do not need a strict correspondence between the parameters of the simulated crystalline structure and the experimental data. It is believed that the CSFF (CHARMM 22) set of full-atom potentials [31] allows one to correctly simulate the mobility of hydroxyl groups located on the surface of a cellulose nanofibril immersed in an aqueous solution. There exist other force fields parameterized for cellulose modeling in an aqueous environment (CHARMM35 [29], GLYCAM06 [30], Gromos 45a4 [31], and OPLS [32]). However, their use results in a fibril twisting as early as after 10 ns of modeling [26, 3334]. Though the fibril twisting disappears, the transformation into the untwisted state takes a very long simulation time when these force fields are used (C35 - 400 ns, GLYCAM06 - 600 ns and Gr45a4 - more than 1 mks) [24]. The final crystal structure is stabilized by a three-dimensional network of intermolecular hydrogen bonds, in contrast to the two-dimensional network which is typical of the Iβ crystal structure. When CSFF is used, the structure of the same type is obtained, but a much shorter simulation time is needed (~10 ns) [20, 35]. This was an additional argument in favor of our choice of CSFF.

ii. For water For water, the CHARMM version of the TIP3P model (with the Lennard-Jones parameters for hydrogen and oxygen) [36] and the TIP4P-ew model [37] (parameterized for the use with the Ewald summation methods) were selected. It is known that the simulation results may depend on the use of one or another type of water models. The necessity to consider the systems with different water models was due to the fact that ions can be hydrated differently in such systems. It was shown in [38] that monovalent ions are hydrated differently when the TIP3P or TIP4P-ew water models (widely used for biomolecules) are applied, which can influence the

ACS Paragon Plus Environment

7

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29

formation of ion contacts with each other and with the BC fibril surface. In the case of supersaturated solutions of multivalent ions the interactions between ions prevail, and the role of water models in modeling such systems is less important. For this reason, only the TIP3P water model was used for modeling the solution with the HAP ionic composition. However, two water models (TIP3P and TIP4P-ew) were considered for the CaCl2 solutions. Though the CSFF force field is parameterized on the use of the TIP3P water model, it has been shown (for example, for short peptides [39]) that Charmm 22 is also compatible with the TIP4P-ew water model.

iii. For ions The parameters for the ions that make up hydroxyapatite ((PO4)3-, OH1- and Ca2+) were taken from [40]. It is important that the partial atomic charges for the (PO4)3- and OH1- ions were obtained by using the same quantum-chemical approach and the same set of the basis functions as for the cellulose molecule in the CSFF force field. The parameters for the ions in the CaCl2 solution were taken from the parameter set of the CHARMM 27 potentials [41].

D. Method, Software, and Procedure Modeling was performed via the method of molecular dynamics with the use of the program GROMACS 4.5.5 [42, 43]. The main modeling stages are briefly described in Supporting Information. All MD simulations were performed at atmospheric pressure in a NPT assembly at 300 K. To maintain constant temperature and pressure in the system, a Nose–Hoover thermostat and a barostat with temporal constants τT = 0.2 ps and τp = 0.5 ps, respectively, were used. Long-distance Coulomb interactions were calculated via the Ewald summation procedure. The C–H valence bond length in the cellulose molecule was limited with the LINCS algorithm [44]. Water molecules were maintained rigid throughout the process of modeling with the use of the SETTLE procedure [45]. The results were obtained by averaging over the trajectory in the interval 50–100 ns because, after modeling for 50 ns, the energy and the volume fluctuated about

ACS Paragon Plus Environment

8

Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

certain mean values. The scale of energy and volume fluctuations was ~0.1% and ~0.5% of the mean value, respectively. The dependence of the results on the starting configurations was examined by modeling three systems having the same number of water molecules and ions but different initial positions and velocities for each solution of interest. Discrepancies in the results did not exceed 3%.

III RESULTS A. Cellulose nanofibril structure and conformations. Hydration of the fibril surface. Effect of water model. As pointed out above, the use of any set of the potentials parameterized for cellulose in water does not lead to the simulated crystalline structure of cellulose with the parameters, which would completely correspond to the experimental values. A significant change in the fibril volume (“swelling”) was observed, and it was the same in the systems with the TIP3P and TIP4P-ew water models. This “swelling” was caused by the parameterization of the CSFF force field potentials directed toward the correct description of the behavior of the cellulose primary hydroxyl groups in a contact with water. Such a parameterization results in an equilibrium distribution of rotational isomers of primary hydroxyl groups which substantially differs from the distribution characteristic of the Iβ crystalline state of cellulose [18]. It was found that the fibril retains its integrity and a highly ordered structure, but it is stabilized by a three-dimensional network of H bonds instead of a two-dimensional one which exists in the Iβ structure. In our earlier work [35], changes in the unit cell volume were discussed in detail. Here we focus only on major results. Note that the obtained transverse fibril “swelling” results in some increase in the distances between hydroxyl groups of neighboring chains on the cellulose fibril surface. This increase is about 3% for those surfaces of the cellulose nanofibrils which form the nanoribbon surface. These changes in the simulated structure cannot affect the result in estimation of the ability of the hydroxyl groups located on native BC fibril surfaces to adsorb Ca ions. However, in the problems for which the degree of the structural conformity between two components is ACS Paragon Plus Environment

9

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

important, the changes in the distribution of primary hydroxyl groups should be taken into account. In order to evaluate the conformations of the cellulose molecules forming the fibril surface, the distribution functions of the dihedral angles averaged over the trajectory (50-100 ns) as well as over the array of molecules on the fibril surface were obtained. The dihedral angles defining the conformation of the primary hydroxyl groups and the conformation of the main chain of the cellulose molecule for two water models are presented in Supporting Information in Figure S1. We found that the conformational characteristics of the cellulose molecules for both water models were nearly the same. The characteristics of water interaction with the fibril surface (radial distribution functions g(r) for water molecules at hydroxyl groups, the degree of hydration of hydroxyl groups, and the potentials of mean force calculated from g(r)) are shown in Supporting Information in Figure S2. It is evident from this figure that the interaction of water molecules with the cellulose fibril surface and, hence, the degree of hydration have minor differences for both water models. Water molecules mostly interact with the primary hydroxyl groups. This is due to the fact that the primary hydroxyl groups are more accessible to water molecules owing to their greater distance from the glucose cycle. There is approximately one water molecule per primary hydroxyl group in the first coordination layer. Our results agree well with the findings of other authors [46, 47]. It is very important that the distribution functions for the dihedral angles of the cellulose molecules on the fibril surface in ionic solutions are the same as the functions obtained for the fibril in contact with pure water. Moreover, neither the distribution functions for water molecules nor the quantity of water molecules on the fibril surface undergo any changes due to the ions in the solutions.

ACS Paragon Plus Environment

1 10

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

B. Cellulose nanofibril in the solution with the HAP ionic composition Interactions between ions and the fibril surface can result in their adsorption on the surface, while their interactions with each other can lead to their aggregation (formation of clusters). An ionic cluster is defined as a complex of ions in which (i) each ion is bonded to at least one ion and two ions are considered bonded if they are at a distance shorter than a certain distance of contact, rc, and (ii) all ions are bonded to each other by a series of connections. The contact distance was selected from the pair radial distribution function for the ion pair. This distance determines the position of the first peak border and corresponds to the strict condition that only ions are included in the series of connections. The cluster size was determined by the quantity of incorporated ions. Figure 2 presents the cross section of the computational cell with the cellulose fibril and the ions in the initial state, after 12 ns, and after 100 ns of modeling (water is not shown in order not to complicate the pattern.)

a) b) c) Figure 2. Distributions of ions in the modeling cell: (a) in the initial state, (b) in 12 ns, and (c) at the end (100 ns) of simulation. Red balls – O atoms, cyan balls – Ca2+ ions, yellow balls– P atoms, white balls – H atoms, green balls – cellulose atoms.

Initially ions form small-size (3-5 ions) clusters which in the process of modeling form larger clusters (Figure 2b). After 100 ns of the simulation they form one cluster which has only a short-range order in the structure and which does not interact with the fibril. The fibril surface is surrounded only by water molecules (Figure 2c). (The pair distribution function for the Ca2+ and (PO4)3- ions and the mutual arrangement of ions in the final cluster are presented in Supporting Information, Figure S3). ACS Paragon Plus Environment

11

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

The results obtained in our simulation demonstrate that the clusters are not adsorbed on the cellulose fibril surface. As shown earlier, both experimentally and by computer modeling [48, 49], the preliminarily synthesized HAP crystallites can be adsorbed on BC fibrils. However, our simulation time was not long enough to get the transition from an amorphous to the crystalline structure. We also discovered that not only the clusters but also individual Ca2+ ions are not adsorbed on the cellulose fibril surface. This cannot be explained by the fact that the solution contains trivalent anions (PO4)3- and, obviously, hydroxyl groups on the fibril surface cannot compete with them for the formation of contacts with calcium ions. However, as the simulation results presented in the next section show, Ca2+ ions in the CaCl2 solutions (with monovalent anions) are not adsorbed on the cellulose fibril surface as well.

C. Cellulose nanofibril in CaCl2 solutions i. System evolution in time Figure 3 shows the cross section of a simulation box with the cellulose fibril and ions (concentration C3) at the initial moment (Figure 3a) and after 100 ns of the simulation (Figure 3b) (not to complicate the picture, water is not shown).

a) b) Figure 3. Cross-section of the simulation box with the cellulose fibril and ions: (a) at the beginning and (b) after 100 ns of MD simulations. Cyan balls – Ca2+ ions, grey balls – Cl1- ions, green balls – cellulose atoms.

Figure 4 shows the fractions of free ions as functions of time for salt concentrations C=C1, C=C2, and C=C3. The data were taken at each 100 ps of the trajectories from the ACS Paragon Plus Environment

1 12

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

simulations with the TIP3P and TIP4P-ew water models. The initial period (0 – 1 ns) is presented separately in Figure 4b.

a)

b)

Figure 4. Changes in quantities of free ions (in %) during the modeling process (a) for the total simulation time and (b) for the time interval 0-1 ns. The data are presented for the CaCl2 solution with concentrations C1, C2, and C3 for the simulations with the TIP3P and TIP4P-ew water models.

As one can see from Figure 4a, major changes occur at the very beginning of the simulations. In the interval 50–100 ns, where all the results were obtained with averaging, the standard deviations of the free ion quantities did not exceed 0.5%. At the concentrations corresponding to the solubility limit or above it, major changes occur within the first 100 ps and more than 80% of the ions are in a bound state, while at the low concentration the amount of free ions decreases sufficiently slowly. The rate of ion transition into a bound state depends on the water model. In the system with TIP4P-ew the ion-ion contacts are formed faster.

ii. Ca2+ ions on the nanofibril surface. Effect of water model. In order to find at what concentrations and how many Ca2+ ions are bonded to the fibril surface, the radial distribution functions (g (r)) of the Ca2+ ions with respect to oxygen atoms of surface hydroxyl groups were calculated for solutions with different concentrations for two models of water. The radial distribution functions (g (r)) and the number of the Ca2+ ions on the cellulose fibril surface are shown in Figure 5 for two water models. In Figure 5 and in the figures ACS Paragon Plus Environment

1 13

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 29

where three concentrations are presented the solid line corresponds to the concentration much lower than the solubility limit (C1); the dotted line is used for the concentration equal to the solubility limit (C2); and the dashed line is for the concentration exceeding the solubility limit (C3).

a)

b)

Figure 5. (a) Radial distribution functions (averaged over time) of Ca2+ ions at the oxygen atom of hydroxyl groups and (b) number of Ca2+ ions on the cellulose fibril surface for two water models.

The height of the first peak corresponding to the closest contact of Ca2+ ions with oxygen atoms of the hydroxyl groups increases with increasing concentration. Regardless of the water model chosen, no Ca2+ ions are observed on the fibril surface at the low concentration. As the salt concentration increases, calcium ions appear on the surface. There are some differences for two water models which indicate that in the case of TIP4P-ew the number of the closest contacts of the Ca2+ ions with the fibril surface are greater than for TIP3P. It should be emphasized that Ca2+ ions do not interact with secondary hydroxyl groups. The average number of Ca2+ ions contacting with one monomer (or with the primary hydroxyl group) of the cellulose molecule on the fibril surface is insignificant (Figure 5b). The evaluation shows that even at the highest concentration (C=C3) there are ≈2 (TIP3P) or ≈14 (TIP4P-ew) Ca2+ ions in the close contact with the entire surface of the simulated fibril (74.85 nm2), i.e., about 0.13% and 0.9% of the total amount of calcium ions in the system. In Section IV we discuss the reasons for different numbers of Ca2+ ions on the cellulose fibril surface obtained for ACS Paragon Plus Environment

1 14

Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

the systems with two water models. Though the number of the Ca2+ ions contacting with the fibril surface is very small, they could serve as centers for subsequent crystallization if they were adsorbed on the surface. The potentials of mean force for the interaction between Ca2+ and oxygen of the hydroxyl groups on the fibril surface calculated from g(r) and by using the umbrella sampling method [50] (see details in Supporting Information) are presented in Figure 6.

a)

b)

Figure 6. Potentials of mean force for interactions of Ca2+ ions with oxygen of hydroxyl groups calculated (a) from g(r) at the concentration C3 (for the TIP3P and TIP4P-ew water models) and (b) by using the umbrella sampling procedure and the CSFF (solid black line) and GLYCAM06 (dashed black line) force fields with the TIP3P water model.

The potentials of mean force have minor differences for the systems with the TIP3P and TIP4P-ew models as shown for the PMF calculated from g(r) (Figure 6a). However, the potentials calculated from g(r) (Figure 6a) and by using the umbrella sampling method (Figure 6b) are different. The PMF calculated from g(r) has the energy barrier that hinders the escape of Ca2+ from the fibril surface, and the PMF calculated by using umbrella sampling has no such a barrier. The g(r) function is a characteristic of the system structure. At saturated concentrations, after the ionic clusters are formed, the system viscosity increases dramatically, as can be seen from the root mean square displacement for Ca2+ ions in the systems with concentrations C1, C2 and C3 presented in Figure 7.

ACS Paragon Plus Environment

1 15

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 29

Figure 7. Mean-square displacement of Ca2+ ions as functions of time for CaCl2 solutions with concentrations C1, C2 and C3.

The diffusion coefficients obtained from the linear parts of the MSD(t) curves in the interval 50÷100 ns for different concentrations are as follows: 0.8×10-5 cm2/s for C=C1, 1.81 ×10-7 cm2/s for C=C2 and 3.14×10-8 cm2/s for C=C3. This means that at high concentrations, after the cluster formation is completed, we have a system in which mixing of its components cannot occur during our simulation because it requires much more time. At C=C3, for example, the calcium ion displacement at the distance of at least half the simulation box requires about 1.5 microseconds of the simulation time according to the estimated values of the diffusion coefficient. Thus, the duration of our simulations (after the ionic clusters were formed) was not enough for a further evolution of the system (for mixing of its components). As a result, we obtained averaging over the “frozen” structure, and g(r) characterized only the particle distribution in this structure. Thus, the presence of the first minimum in the PMF (Figure 6a) calculated from g(r) indicates that there are some ions entrapped on the surface and these ions cannot go away. The presence of the barrier at the distance corresponding to the position of the maximum on the PMF curve shows that there are no ions in this layer because it is the distance of the water-surface close contact (see Figure S2c in Supporting Information). We believe that long-time simulations will give the same result as that obtained in our modeling for the low concentration (C=C1), namely, the absence of Ca ions on the surface (see Figure 5b). The ACS Paragon Plus Environment

1 16

Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

dependence of PMF on the distance between the Ca2+ ion and the surface calculated by using the umbrella sampling algorithm (see Fig. 6b, solid line) demonstrates exactly this finding and shows that Ca2+ ions prefer to be in the bulk rather than on the surface. No doubt, the result can depend on the force field parameters and, first of all, on the partial charges on an oxygen atom of a primary hydroxyl group. Partial charges on oxygen atoms of primary hydroxyl groups are somewhat different in different force fields (See Table S3 in Supporting Information). In GLYCAM06, the largest negative charge is assigned to the oxygen atom of the primary hydroxyl group. By using the GLYCAM06 force field and the umbrella sampling method, we carried out the simulations to calculate the dependence of the PMF on the distance between Ca2+ and the cellulose fibril surface. We obtained (the dashed curve in Figure 6b) the dependence of the same type as by using the CSFF force field. The discrepancy in the PMF values at large distances (in the bulk) is due to the fact that the resulting curves are shifted to assign zero value for each system. The discrepancy is caused by the difference in the interaction of the Ca2+ ion with the primary hydroxyl group (~4.5 kJ/mol) in the simulations using the CSFF or GLYCAM06 force fields. To summarize, it can be concluded that hydroxyl groups on the native BC fibril surface cannot adsorb Ca2+ ions from solution. The formation of ionic clusters proved to be more favorable in solution than on the fibril surface.

iii. Ions in solution. Effect of water model. Since a decrease in the number of free ions in the system with time is not related to their adsorption on the fibril surface (as shown above), their transition into a bound state points to the formation of ionic clusters. Figure 8a presents the evolution of the average cluster size (averaged over all clusters and expressed through the quantity of ions contained in a cluster) with time for the solution concentration C=C1. To estimate the cluster size, the contact distance of 0.27 nm was determined from the Ca-Cl pair radial distribution function (see Figure S5a in Supporting

ACS Paragon Plus Environment

1 17

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 29

Information).

a)

b)

Figure 8. Evolution of the average cluster size in time (a) at C=C1 and (b) at C=C3 for the TIP3P (red line) and TIP4-ew (black line) water models.

As one can see from Figure 8a, the average size of the clusters which are formed in the system at C=C1 with TIP4P-ew is slightly larger than in the system with TIP3P (3.9 vs 3.3). This is due to the fact that in the system with the TIP4P-ew model the numbers of free ions and ionic pairs and triples are smaller, but the number of larger clusters is greater (see Figure S6 in Supporting Information). On the contrary, in the supersaturated solution (C=C3) the average cluster size in the system with the TIP4P-ew model is somewhat smaller (Figure 8b), but the distributions strongly overlap and thus the cluster sizes are essentially equal. We discuss the nature of the discrepancies obtained for two water models in Section IV.

IV. DISCUSSION i. Comparison with the classical theory of a heterogeneous crystallization The results obtained in our simulations suggest that in the presence of native BC fibrils the HAP nucleation process occurs in solution rather than on fibril surfaces (heterogeneously). According to the classical theory of heterogeneous crystallization, the function characterizing the induction ability of a foreign particle (when it can be considered to be a flat surface with respect ACS Paragon Plus Environment

1 18

Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

to the critical nucleus) is the function of surface tension energies between phases: f(m)=(1/4)*(23m+m3) where m = (γsf - γsc)/γcf and γsc is for BC-HAP, γsf is for BC-water, γcf is for HAP-water (the main formulas from [51, 52] that describe this process are presented in Supporting Information). We can estimate the effect of the substrate (such as BC) on the formation of the mineral phase (HAP) from solution by taking the experimental values γ of surface tension coefficients for BC, HAP, and water [53, 54] (Table S4 in Supporting Information). As a result, we obtain m = -0.51 and f = 0.85. These values show that the native BC still has some induction ability, but this ability is low. The resulting estimate is based on the main statement of the classical theory of crystallization that the germ of a new phase formed heterogeneously (on a substrate surface) has a crystalline structure and for this reason the surface free energy of a crystallite is considered in the theory. In our simulation, we obtained the amorphous clusters formed in solution rather than adsorbed on the cellulose fibril surface. Not only the clusters but also individual Ca2+ ions were not adsorbed on the cellulose fibril surface either in calcium phosphate or CaCl2 solutions.

ii. The nature of discrepancies between the results for the TIP3P and TIP4P-ew water models It was pointed out above that the discrepancies between the results obtained by using the TIP3P and TIP4P-ew water models concern the average cluster sizes (Figures 8a and 8b, at low and high salt concentrations, respectively), the rate of transition of ions into a bound state (Figure 4), and the amount of Ca2+ ions on the cellulose fibril surface (Figure 5b). The cluster formation depends on the competing interactions between all components of the system (ion-surface, surface-water, ion-ion, ion-water and water-water interactions). Since the systems considered differ only by the water models, the ion-ion and ion-surface interactions are the same for both systems. In addition, as shown in Figure S2, the interaction of water molecules with the fibril surface is nearly independent of the water model. Thus, it can be concluded that the discrepancies in the results are due to the differences in the ion-water and water-water

ACS Paragon Plus Environment

1 19

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 29

interactions. At low salt concentrations, the interactions between water molecules make the main contribution to the free energy of the system. The TIP3P water molecules are bonded weaker, and this is the main reason for the differences in the cluster sizes for two models of water. Stronger interactions of the chlorine ions with the TIP3P water molecules additionally promote the formation of the clusters having a smaller average size. As noted above, the rate of ion transition into a bound state at low concentrations depends on the water model. We have found that in the system with TIP4P-ew the ion-ion contacts are formed quicker mainly because the energy barrier for removal of a water molecule from a calcium ion (in order to form the ion pair with a chlorine ion) in the case of the TIP4P-ew water model is lower than in the case of the TIP3P model. (The potentials of mean force for the water-water and ion-water interactions are presented in Figure S7 in Supporting Information). At the high concentration (C=C3), the amount of calcium ions on the cellulose fibril surface in the system with TIP4P-ew is greater than in that with TIP3P (Figure 5b). Since the number of water molecules in the supersaturated solution is low, the above differences are due to the differences in the ion-water interactions. The calculations of the potential energies for the Cl1-water and Ca2+-water pairs showed that the potential energy minimum for the interaction of the calcium ion with a molecule of the TIP3P water model is deeper than for the interaction with a molecule of the TIP4P-ew water model. Thus, in the system with TIP4P-ew the calcium ions having a weaker interaction with water molecules can appear on the fibril surface in a somewhat greater amount than in the case of TIP3P (as our results demonstrated). Detailed information concerning the potential energies of the Cl1--water and Ca2+-water interactions is presented in Supporting Information, in Table S5).

iii. Average cluster size and force field parameters for ions The average cluster size in all the simulated systems turned out to be larger that that

ACS Paragon Plus Environment

2 20

Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

derived from experimental data [55] the interpretation of which speaks in favor of ionic pairs at low concentrations. The reason for the discrepancy between the simulated and experimental results is that the parameters for the ions incorporated into the CHARMM force field were optimized only for their interactions with biomolecules and were not additionally optimized for the correct description of their behavior in solution. Selection of correct parameters for ions is a difficult task. There are no direct experimental data for constraining the ion-ion interactions. As a result, there are substantial variations in ion-ion interaction parameters among various force fields. One way to overcome the arbitrariness in ion-ion interaction parameters is to appeal to the higher level ab initio calculations [56]. The ab initio PMF calculated in [56] was compared to the classical PMFs (See Figure S8 in Supporting Information) obtained by using the parameters of two commonly used classical force fields CHARMM (used in our study) and DS95 [57] (which is typically used to model salt solutions). It has been demonstrated that the CHARMM PMF reproduces the positions of minima quite well and matches the ab initio PMF, but somewhat overestimates the energy of the close contact between the ions. In the case of DS95 the energy of the close contact is underestimated but the contact distance is overestimated. It was attempted in [58] to improve the results by (i) using polarizable force fields or (ii) correcting the parameters so as to obtain the contact distance and energy corresponding to those derived from accurate electronic structure calculations. The first approach resulted in a much greater overestimation of the contact energy than by using CHARMM, but the second one strongly underestimated it (Figure 8s in Supporting Information). Thus, at present time CHARMM provides a more realistic parametrization for aqueous CaCl2. Of course, the parameterization should affect the Ca2+ interaction with the fibril surface. Nevertheless, even if we had the parameters which strictly reproduce the ab initio PMF, this would not change the results obtained by using CHARMM for the interaction of calcium ions with the cellulose fibril surface. In this case we will have not only a weaker Ca-Cl interaction, but also a weaker interaction between the calcium ions and the organic matrix (due to a less deep

ACS Paragon Plus Environment

2 21

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 29

potential well and larger contact distance). Thus, our main conclusion that the cellulose matrix does not adsorb calcium ions is still valid.

V. CONCLUSIONS We have simulated the cellulose nanofibril mineralization, for the first time to our knowledge. The simulations of the behavior of a BC fibril immersed in calcium phosphate and CaCl2 aqueous solutions for 100 ns at 300 K under atmospheric pressure were performed via the molecular dynamics method. Sets of the CSFF (CHARMM 22) and GLYCAM06 potentials for the cellulose molecules, CHARMM 27 for the ions, and the TIP3P and TIP4P-ew models for water molecules were used. The distribution functions for the dihedral angles of internal rotation of cellulose molecules on the surface of the fibril (placed into ionic solutions) were found to be the same as those obtained for the fibril in pure water. In other words, no effect of the ions present in the solutions on the cellulose molecules conformations was observed. Moreover, the distribution functions for water molecules and the quantity of water molecules on the fibril surface were unaffected by ions in the solutions. We have also found that neither any clusters nor individual Ca2+ ions are adsorbed on the cellulose fibril surface. Hydroxyl groups on the surface of the native BC fibril cannot adsorb Ca2+ ions from the solution. The formation of ionic clusters proves to be more favorable in solution than on the fibril surface. Therefore, it can be concluded that the pre-incubation of BC matrix in a CaCl2 solution cannot create centers for the subsequent heterogeneous crystallization of CP minerals on cellulose fibril surfaces. Consequently, the mineralization process cannot be caused by adsorption. It can be caused by absorption of ions (physical entrapping in the volume) due to the three dimensional lattice structure of the organic matrix. Thus, when native BC is used to produce BC-CP composite materials, the CP crystal formation should be regarded as the

ACS Paragon Plus Environment

2 22

Page 23 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

process of homogeneous crystallization (crystal nucleation and growth occur in solution), and then the crystallites formed can be adsorbed on the native BC fibril surface. Our goal was to study the ability of hydroxyl groups (located on surfaces of native BC fibrils) to adsorb Ca2+ ions. To solve this task, a strict compliance of the crystal structure characteristics of the simulated nanofibrils to the experimental data is not required. However, to ensure successful results in studying, for example, the interaction of native cellulose nanofibrils with polymeric (or biopolymeric) matrices in composite materials, it is important to reproduce correctly the experimentally determined nanofibril crystal structure. To this end, the force field for the cellulose nanofibril modeling should include two sets of parameters. One set, for the molecules in the fibril interior, must be chosen in such a way as to preserve the crystalline structure corresponding to the experimental data. The other set of the parameters should reproduce the conformational characteristics of the cellulose molecules of the nanofibril surface contacting with the environment. We are facing the same task in modeling of heterogeneous nucleation of mineral phases on the surfaces of chemically modified cellulose fibrils.

ACKNOWLEDGMENTS Computerized modeling was performed with the resources of SKIF Chebyshev, Lomonosov supercomputers of Moscow State University and Center of collective use "Complex of modeling and data research mega-class facilities" NRC "Kurchatov Institute". Unique identifier RFMEFI62114X0006. This study was supported by federal program of the Russian government, the agreement 14.613.21.0024, unique identifier RFMEFI61314X0024 The reported study was funded by RFBR according to the research project No. 14-03-32041 мол_a.

ACS Paragon Plus Environment

2 23

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 29

References: (1) Olszta, M.J.; Cheng, X.; Jee S.S.; Kumar R.; Kim Y.-Y.; Kaufman M.J.; Douglas E.P.; Gower L.B. Bone structure and formation: A new perspective. Materials Science & Engineering RReports. 2007, 58, 77–116. (2) Ogomi, D.; Serizawa, T.; Akashi, M. Bioinspired organic-inorganic composite materials prepared by an alternate soaking process as a tissue reconstitution matrix. J. Biomed. Mater. Res. 2003, A 67, 1360–1366. (3) Khripunov, A.K.; Baklagina, Yu.G.; Sinyaev, V.A.; Shustikova, E.S.; Paramonov, B.A.; Romanov, D.P.; Smislov, R.Ju.; Tkachenko, A.A. Investigation of nanocomposites based on hydrated calcium phosphates and cellulose. Glass Physics and Chemistry 2008, 34, 248−257. (4) Tazi, N.; Zhang, Z.; Messaddeq, Y.; Almeida-Lopes, L.; Zanardi, L.M.; Levinson, D.; Tazi, M.R. Hydroxyapatite bioactivated bacterial cellulose promotes osteoblast growth and the formation of bone nodules. AMB Express 2012, 2, 61–71. (5) Jonas, R.; Farah, L.F. Production and application of microbial cellulose. Polym. Degrad. Stab. 1998, 59, 101–106. (6) Yamanaka, S; Watanabe, K; Kitamura, N; Iguchi, M; Mitsuhashi, S; Nishi, Y; Uryu, M. The structure and mechanical properties of sheet prepared from bacterial cellulose. J. Mater. Sci. 1989, 24, 3141−3145. (7) Backdahl, H; Helenius, G; Bodin, A; Nannmark, U; Johansson, B.R.; Risberg, B.; Gatenholm, P. Mechanical properties of bacterial cellulose and interactions with smooth muscle cells. Biomaterials 2006, 27, 2141−2149. (8) Brown, R.M.Jr., The Biosynthesis of Cellulose. J. Macromol. Sci. A 1996, 33, 1345−1373. (9) Saska, S.; Barud, H.S.; Gaspar, M.M.; Marchetto, R.; Ribeiro J.L. Bacterial Cellulose – Hydroxyapatite Nanocomposites for Bone Regeneration Int. J. Biomat. 2011, 1−8. (10) Nge, T.T.; Sugiyama, J. Surface functional group dependent apatite formation on bacterial cellulose microfibrils network in a simulated body fluid. J. Biomed. Mater. Res. A 2007, 81, 124−134. (11) Wan, Y.Z.; Huang, Y.; Yuan, C.D.; Raman, S.; Zhu, Y.; Jiang, H.J.; He, F.; Gao, C. Biomimetic synthesis of hydroxyapatite/bacterial cellulose

nanocomposites for biomedical

applications Materials Science and Engineering C 2007, 27, 855–864. (12) Nge, T. T.; Sugiyama, J.; Bulone, V. Bacterial Cellulose-Based Biomimetic Composites

ACS Paragon Plus Environment

2 24

Page 25 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Biopolymers 2010, 346–368. (13) Hutchens, A.S.; Benson, S.R.; Evans, R.B.; Hugh, M.O.; Rawn J.C., Biomimetic synthesis of calcium-deficient hydroxyapatite in a natural hydrogel Biomaterials 2006, 27, 4661–4670. (14) Dumitru, S.; Vidal, P.F.; Chronet, E. Hydrogels based on polysaccharides. In Polysaccharides in Medicinal Applications. CRC Press: New York, 1996, 87–105. (15) Dobbins, R. J.; Solute-solvent interactions in polysaccharide systems. In Industrial gums, polysaccharides and their derivatives. 2nd ed. Academic Press: New York 1973, 19–25. (16) Nishiyama, Y.; Langan, P.; Chanzy, H. Crystal Structure and Hydrogen-Bonding System in Cellulose Iβ from Synchrotron X-ray and Neutron Fiber Diffraction J. Am. Chem. Soc. 2002, 124, 9074–9082. (17) Zhang, R; Ma, P.X. Biomimetic polymer/apatite composite scaffolds for mineralized tissue engineering. Macromol. Biosci. 2004, 4, 100–111. (18) Kuttel, M.; Brady, J.W.; Naidoo K.J. Carbohydrate solution simulations: producing a force field with experimentally consistent primary alcohol rotational frequencies and populations. J. Comput. Chem. 2002, 23, 1236–1243. (19) Palma, R. Molecular Mechanics Studies of Celluloses. J. Am. Chem. Soc. ACS Symp. Ser. 2000, 112–130. (20) Matthews, J.F.; Himmel, M. E.; Brady, J. W. Simulations of the Structure of Cellulose. Computational Modeling in Lignocellulosic Biofuel Production, J. Am. Chem. Soc. ACS Symp. Ser., Ch. 2 2010, 17–53. (21) Matthews, J.F. Molecular mechanics simulations of cellulose and cellobiose. Doctoral Thesis. Cornell University. 2008. (22) Maurer, R.J.; Sax, A.F.; Ribitsch, V. Molecular simulation of surface reorganization and wetting in crystalline cellulose I and II Cellulose 2013, 20, 25–42. (23) Chundawat, S.P.S.; Bellesia, G.; Uppugundla, N.; Da Costa Sousa, L.; Gao, D.; Cheh, A.M.; Agarwal, U.P.; Bianchetti, C.M.; Phillips, G.N.; Langan, P.; Balan, V.; Gnanakaran, S.; Dale, B.E. Restructuring the crystalline cellulose hydrogen bond network enhances its depolymerization rate. J. Am. Chem. Soc. 2011, 133, 11163–11174. (24) Matthews, J.F.; Beckham, G.T.; Bergenstrеhle-Wohlert, M.; Brady, J.W.; Himmel, M.E.; Crowley, M.F. Comparison of Cellulose Iβ Simulations with Three Carbohydrate Force Fields J. Chem. Theory Comput. 2012, 8, 735–748.

ACS Paragon Plus Environment

2 25

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 29

(25) Bergenstrеhle, M. Crystalline cellulose in bulk and at interfaces as studied by atomistic computer simulations. Doctoral Thesis. Stockholm. 2008. (26) Paavilainen, S.; Rog, T.; Vattulainen, I. Analysis of twisting of cellulose nanofibrils in atomistic molecular dynamics simulations. J. Phys. Chem. B. 2011, 115, 3747–3755. (27) Matthews, J.F.; Skopec, C.E.; Mason, P.E.; Zuccato, P.; Torget, R.W.; Sugiyama, J.; Himmel, M.E.; Brady, J.W. Computer simulation studies of microcrystalline cellulose Iβ. Carbohydr. Res. 2006, 341, 138–152. (28) Yui, T.; Nishimura, S.; Akiba, S.; Hayashi, S. Swelling behavior of the cellulose Iβ crystal models by molecular dynamics. Carbohydr. Res. 2006, 341, 2521–2530. (29) Guvench, O.; Greene, S.N.; Kamath, G.; Brady, J.W.; Venable, R.M.; Pastor, R.W.; Mackerell, A.D. Additive empirical force field for hexopyranose monosaccharides. J. Comput. Chem. 2008, 29, 2543–2564. (30) Kirschner, K.N.; Yongye, A.B.; Tschampel, S.M.; Gonzalez – Outeirino, J.; Daniels, C.R.; Foley, B.L.; Woods, R.J. GLYCAM06: a generalizable biomolecular force field. Carbohydrates. J. Comput. Chem. 2008, 29, 622–655. (31) Lins, R.D.; Hunenberger, P.H.A new GROMOS force field for hexopyranose –based carbohydrates. J. Comput. Chem. 2005, 26, 1400–1412. (32) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W.L. OPLS allbatom force field for carbohydrates. J. Comput. Chem. 1997, 18, 1955–1970. (33) Bergenstrеhle, M.; Matthews, J.; Crowley, M.; Brady, J. Cellulose crystal structure and force fields. In Proceedings of International Conference on Nanotechnology for the Forest Production Industry, Otanienmi, Espoo, Finland, 2010. (34) Hadden, J.A.; French, A.D.; Woods, R.J. Unraveling Cellulose Microfibrils: A Twisted Tale. Biopolymers 2013, 99, 746–756. (35) Tolmachev, D.A.; Lukasheva, N.V. Study of the Process of Mineralization of Nanofibrils of Native Bacterial Cellulose in Solutions of Mineral Ions: Modeling via the Method of Molecular Dynamics. Polymer Science, Ser. A 2014, 56, 545–557. (36) MacKerell, A.D.J.; Bashford, D.; Bellott, R.L.; Dunbrack, R.L.J.; Evanseck, J.D.; Field, M.J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F.T.K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D.T.; Prodhom, B.; Reiher, W.E., III; Roux, B.; Schlenkrich, M.; Smith, J.C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.;

ACS Paragon Plus Environment

2 26

Page 27 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. 1998, 102, 3586–3616. (37) Horn, H.W.; Swope, W.C.; Pitera, J.W.; Madura, J.D.; Dick, T.J.; Hura, G.L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665–9678. (38) Fennell, C.J.; Bizjak, A.; Vlachy, V.; Dill, K.A. Ion Pairing in Molecular Simulations of Aqueous Alkali Halide Solutions. J. Phys. Chem. B. 2009, 113, 6782–6791. (39) David, R.N.; Jeremy, C.S. Molecular Dynamics Simulations of Proteins: Can the Explicit Water Model Be Varied? J. Chem. Theory Comput. 2007, 3, 1550–1560. (40) Hauptmann, S.; Dufner, H.; Brickmann, J.; Kast, S.M.; Berry, R.S. Potential energy function for apatites. Phys. Chem. Chem. Phys. 2003, 5, 635–639. (41) MacKerell, J.A.D.; Banavali, N.; Foloppe, N. Development and current status of the CHARMM force field for nucleic acids. Biopolymers 2001, 56, 257–265. (42) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load –Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. (43) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A.E.; Berendsen, H.J.C. GROMACS: Fast, Flexible and Free. J. Comput. Chem. 2005, 26, 1701–1718. (44) Hess, B.; Bekker, H.; Berendsen, H.J.C.; Fraaije, J.G.E.M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463–1472. (45) Miyamoto, S.; Kollman, P.A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 1992, 13, 952–962. (46) Heiner, A.P.; Teleman, O. Interface between monoclinic crystalline cellulose and water: breakdown of the odd/even duplicity. Langmuir 1997, 13, 511–518. (47) Heiner, A.P.; Kuutti, L.; Teleman, O. Comparison of the interface between water and 4 surfaces of native crystalline cellulose by molecular–dynamics simulations. Carbohydr. Res. 1998, 306, 205–220. (48) Baklagina, Yu.G.; Lukasheva, N.V.; Khripunov, A.K.; Klechkovskaya, V.V.; Arkharova, N.A.; Romanov, D.P.; Tolmachev, D.A. Interaction between Nanosized Crystalline Components of a Composite Based on Acetobacter Xylinum Cellulose and Calcium Phosphates. Polymer Science, Ser. A 2010, 52, 419–429.

ACS Paragon Plus Environment

2 27

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 29

(49) Tolmachev, D.A.; Lukasheva, N.V Interactions Binding Mineral and Organic Phases in Nanocomposites Based on Bacterial Cellulose and Calcium Phosphates. Langmuir 2012, 28, 13473–13484. (50) Leach, A.R. Molecular Modelling: Principles and Applications. Pearson Education 2001, p. 580. (51) Liu, X.Y. Generic mechanism of heterogeneous nucleation and molecular interfacial effects. In Advances in Crystal Growth Research. Elseiver Science B: Amsterdam, 2001, 42–46. (52) Bhattacharya, S.N.; Gupta, R.K.; Kamal, M.R. Polymeric Nanocomposites: Theory and Practice. Carl Hanser:Munich, 2008, p 50. (53) Dourado, F.; Gama, F.M.; Chibowski, E.; Mote, M. Characterization of cellulose surface free energy. J. Adhes. Sci. Technol. 1998, 12, 1081−1090. (54) Wan, C.; Chen, B. Synthesis and characterization of biomimetic hydroxyapatite/sepiolite nanocomposites. R. Soc. Chem. Nanoscale 2011, 3, 693−700. (55) Chialvo, A.A.; Simonson J.M. The structure of CaCl2 aqueous solutions over a wide range of concentration. Interpretation of diffraction experiments via molecular simulation. J. Chem. Phys. 2003, 119, 8052–8061. (57) Dang, L.X.; Smith, D.E. Comment on “Mean Force Potential For The Calcium-Chloride Ion Pair In Water” J. Chem. Phys. 1995, 102, 3483–3484. (56) Timko, J.; De Castro, A.; Kuyucak, S. Ab initio calculation of the potential of mean force for dissociation of aqueous Ca–Cl, J. Chem. Phys. 2011, 134, 204510. (58) Dang, L.X.; Truong, T.B.; Ginovska-Pangovska, B. Note: Interionic potentials of mean force for Ca2+-Cl− in polarizable water. J. Chem. Phys. 2012, 136, 126101.

ACS Paragon Plus Environment

2 28

Page 29 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Graphical TOC Entry

ACS Paragon Plus Environment

2 29