Mineralization of phosphorylated cellulose: crucial role of surface structure and monovalent ions for optimizing calcium content

Natalia V. Lukasheva *a, Dmitry A. Tolmachev a and Mikko Karttunen abc
aInstitute of Macromolecular Compounds, Russian Academy of Sciences, Bolshoj pr. V.O., 31, 199004 St. Petersburg, Russia. E-mail: luk@imc.macro.ru; dm.tolmachev@yandex.ru
bDepartment of Chemistry, The University of Western Ontario, 1151 Richmond St, London, N6A 5B7, Ontario, Canada. E-mail: mkarttu@uwo.ca
cDepartment of Applied Mathematics, The University of Western Ontario, 1151 Richmond St, London, N6A 5B7, Ontario, Canada

Received 12th September 2018 , Accepted 26th November 2018

First published on 27th November 2018

Cellulose can be phosphorylated to produce organic matrices with highly adsorptive properties for, e.g., biocompatible materials for biomedical applications. We focus on the effects of phosphorylation of surfaces of crystalline nanocellulose and, in particular, on the competitive adsorption of mono- and divalent cations (Na+ and Ca2+) typically contained in mineralizing salt mixtures using all-atom molecular dynamics (MD) simulations. Phosphorylation was applied at 12% and 25% both in water and CaCl2 solutions. Our main result shows that Na+ and Ca2+ cations are concentrated in different interfacial layers with Na+ ions penetrating much closer to the surface. This behavior cannot be described by the Poisson–Boltzmann theory or implicit solvent simulations. Our analysis shows that the physical origin of this observation is due to a balance between the electrostatic interactions and hydration free energy associated with the ions. Adsorption levels of the different ions also respond differently to changes in the degree of phosphorylation. We show that the number of adsorbed Na+ ions per phosphate group increases whereas the number of adsorbed Ca2+ ions decreases with an increasing degree of phosphorylation (or when the number of binding sites increases). The decrease in the number of adsorbed Ca2+ ions can be explained by an increasing “charge–charge” repulsion between the Ca2+ ions attracted by the charged surface. Importantly, our results demonstrate the existence of an optimum degree of phosphorylation in terms of adsorbed Ca2+ ions and can be used as a guideline in materials design, for example, when choosing the cellulose matrix or with other similarly structured biomolecular and polymer surfaces.


Chemically modified bacterial cellulose (BC) is an attractive adsorptive material for biomacromolecules, metal ions and minerals.1–3 It is produced by replacing cellulose hydroxyl groups with functional groups that ionize in aqueous solutions. Due to its adsorptive properties, BC is used as gel-film in the production of hybrid organo-mineral materials such as nanocomposites based on cellulose and calcium phosphates, for example, for bone implants.4 For recent reviews on different aspects of BC see ref. 1 and 5–8.

The most common method for producing chemically modified BC to synthesize hybrid organo-mineral nanocomposites is phosphorylation2,9,10 during which the primary hydroxyl groups are mainly replaced by phosphate groups (Fig. S1, ESI). Structurally, BC gel-film is composed of crystalline nanocellulose, also called cellulose nanofibrils, which are aggregates of inter- and intramolecular hydrogen bonded (H-bonded) cellulose chains. The presence of a large number of charged groups on the surfaces of crystalline nanocellulose may change their morphology and structure, as well as mechanical properties and the structure of the mineral phase formed on their surfaces.1 Electron microscopy studies of phosphorylated BC have shown, however, that grafting of phosphate groups up to the degrees of phosphorylation of about 25%2 does not change the fibril morphology.10 Moreover, it has been shown that phosphorylation does not change the interlayer distance in the crystalline packing of cellulose molecules.10 It is known from literature that for other types of cellulose11–13 higher degrees of phosphorylation – up to 43–50% – are possible. However, such a high degree of phosphorylation leads to an increased swelling of fibers and an easier delamination of the fiber wall.

Despite the above, the effects of phosphorylation on the conformations of the individual cellulose molecules within the surface layer of crystalline nanocellulose remain largely unknown. In addition, due to the lack of experimental data, there is currently no information available about the structure of the interfacial layer between phosphorylated nanocellulose and water or salt solutions, yet such information is crucial for developing and controlling adsorption. Detailed knowledge of the adsorption of (mineral) ions is vital for understanding the mechanisms of mineral phase formation and in synthesis of organo-mineral composites; mineralization has also been shown to be strongly dependent on the nature of ions and electrostatic interactions.14,15

Salt mixtures used in organo-mineral composite synthesis contain, in particular, calcium and sodium cations. These cations have similar sizes but different valences and investigation of their concurrency for the adsorption sites on structured surfaces can be considered as a first step to understand the detailed mechanisms of the mineral phase formation. In addition, and importantly, ions can influence the adsorption of other molecules on mineral and phosphorylated surfaces – the role of valence is also very clear from the Poisson–Boltzmann type calculations from different systems. The importance of electrostatic interactions in a wide range of systems is discussed, e.g., in the review by Levin.16

Atomic-level MD simulations are a unique tool that can provide information about the details of the absorption mechanisms and experimentally testable predictions as has been previously shown in the context of biomineralization.17,18 Previous studies have focused on modeling of cellulose nanofibrils with surface modifications by grafted oligomeric molecules (typically poly(lactic acid)19), and grafted functional groups such as sulfonic and carboxyl ones.20,21

We use atomic level MD simulations to study the effects of phosphorylation on the surface structure of crystalline nanocellulose in both presence and absence of a calcium (CaCl2) salt solution as well as to study the competitive adsorption of Na+ and Ca2+ counterions. Our results show that the degree of phosphorylation has a strong effect on the conformations of the surface chains when CaCl2 is present. No such changes were observed in pure water. In addition, it was found that Na+ and Ca2+ cations concentrate in different regions of the interfacial layer with Na+ ions localized closer to the surface. We also found that adsorption levels (number per phosphate group) of the two counterions varied differently with the degree of phosphorylation. In the system with the higher degree of the substitution, the number of Ca2+ ions adsorbed per phosphate group becomes smaller while the Na+ ions show opposite behavior. The behavior of ions is a complex matter and cannot be fully captured by simplified models such as the classical Poisson–Boltzmann approximation16 or (implicit solvent) Monte Carlo simulations22,23 for they do not include, for example, correlations, the finite sizes of ions, solvation free energy or the atomistic nature of surfaces. Deviations with complex behavior have been reported in many related systems including lipid bilayers24 with charged lipids25 and in crystal growth,26 and biomineralization.27

Methods and models

System set up and force fields

It is known from experimental data that for BC the degree of phosphorylation typically does not exceed 25%.2 Motivated by this, we considered two experimentally relevant cases: 12% and 25%. The cellulose nanofibril was modelled as a crystalline layer based on crystallographic data of cellulose Iβ with unit cell parameters a = 0.82 nm, b = 0.78 nm, c = 1.038 nm, β = 90°, α = 90° and γ = 96.6°, and atomic coordinates of cellulose molecules in a crystal cell.28 Layer thickness was 2.3 nm (Fig. 1a). The layer consisted of 64 cellulose oligomers, each containing 16 glucose residues. To model phosphorylation, part of the primary hydroxyl groups (Fig. S1, ESI) on both surfaces were replaced by ionized phosphate groups (with charge −1) such that the charged groups were evenly distributed on each surface (Fig. 1b and c). The modified layer was placed into the simulation cell which was filled by water molecules and ions. These structures were used to obtain the initial configurations for simulations as follows: they were energy minimized using the conjugate gradient method in order to relax possible steric clashes created at the system construction. The MD simulations were performed as describe in section “Simulation procedures” below. The cell dimensions were 8.9 nm along the cellulose molecules and 8.7 nm and 10 nm in the two transverse directions. The thickness of the water layer was 7.7 nm. This distance is large compared to the length scales characterizing electrostatic interactions in the simulated systems (such as the lateral separation for counterions and Gouy–Chapman length; details are provided in ESI).
image file: c8cp05767b-f1.tif
Fig. 1 Side view of the model cellulose layer in the simulation box (a). Spatial distribution of charged groups at 12% (b) and 25% (c) phosphorylation. Cellulose molecules are orientated along the z-axis. Red: phosphate group oxygens, blue: oxygens of the primary hydroxyl groups. Coordinate axes are shown in the upper left corner for clarity.

Due to the length scales accessible by atomistic simulations, they cannot differentiate between nanofibers and nanocrystals since the main differences between them are morphological;29 at larger scales, nanofibers are long and flexible whereas nanocrystals are short and more rod-like. At the length scales studied here, there is no difference between the two. Since nanofibers and nanocrystals have the same crystal structure and, correspondingly, their surface structures are the same, ion adsorption can be assumed to be similar at the nanometer scale.

The following systems were studied: (1) phosphorylated cellulose at 12% in water with counterions (the total number of atoms in the system is 79[thin space (1/6-em)]855), (2) phosphorylated cellulose at 25% in water with counterions (79[thin space (1/6-em)]766 atoms), (3) 12% phosphorylated cellulose in CaCl2 solution (78[thin space (1/6-em)]835 atoms) (4) 25% phosphorylated cellulose in CaCl2 solution (78[thin space (1/6-em)]746 atoms). The solution concentration was 0.5 mol kg−1 (170 Ca2+, 340 Cl and 18[thin space (1/6-em)]887 H2O). The surface charge densities were 0.033 C m−2 (0.2 e nm−2) and 0.066 C m−2 (0.4 e nm−2) at 12% and 25% phosphorylation, respectively. 16 (12%) and 32 (25%) Na+ ions were added to ensure charge neutrality of the simulation box.

All the simulations were performed with the Gromacs simulation suite version For ions, the CHARMM2731 force field was used and CSFF32 (Carbohydrate Solution Force Field) was employed for cellulose. This combination of force fields has been successfully used in previous studies of native cellulose in CaCl2 solution.33

To model phosphorylated cellulose, additional parameters are needed for the deprotonated phosphate groups. The valence bond parameters as well as the valence and dihedral angle potentials for the phosphate groups were taken from CHARMM27.31 To obtain the charge distribution of the deprotonated phosphate group, quantum chemical calculations using the Hartree–Fock method were done with different basis sets (6-31G(d), 6-31+G(d), 6-31G(d,p), 6-31+G(d,p)). Mulliken charge distribution was used for partial charge assignment since Mulliken charges are generally used in CHARMM force fields.31 The charge distribution for the phosphorylated cellulose molecule was obtained from calculations for two cellulose oligomers differing in their mutual arrangements of the glucopyranose residues with substituted and unsubstituted primary hydroxyl groups (Fig. S2, ESI). The criterion for choosing the basis set (6-31+G(d,p)) and the corresponding partial charges was the requirement to preserve intrafibrillar hydrogen bonds. Namely, the charges on the atoms forming these bonds should remain the same as for native cellulose in CSFF. The details are provided in ESI (see Tables SI, SII, and Fig. S3). In addition, we also determined charges based on the Löwdin population analysis.34 Details are provided in ESI and discussion is included in Summary and Conclusions. Importantly, the qualitative observations are independent of the charge assignment although there are some quantitative differences.

Simulation procedures

(1) Each system was simulated for 200 ns with 1 fs time step. All simulations were performed using a fully periodic simulation box in the NPT ensemble at P = 1 bar and T = 300 K.

(2) Constant temperature was maintained by the Nosé–Hoover thermostat35,36 and pressure coupling was done semi-isotropically using the Parrinello–Rahman barostat37 with coupling times of τT = 0.2 ps and τP = 0.5 ps, respectively.

(3) Long-distance Coulomb interactions were calculated using the particle-mesh Ewald method.38

(4) TIPs3P39 water molecules were kept rigid throughout the process using SETTLE;40 TIPs3P is the CHARMM compatible version of the TIP3P model.

(5) The first 100 ns were considered as equilibration, and data collection and averaging were done over the interval 100–200 ns.

(6) The main criterion for equilibration was the analysis of the time-evolution of the radial distribution functions for cations at phosphate groups (see Fig. S4, ESI). Other quantities were also monitored to verify that equilibrium had indeed been reached. Typical error bars in our simulations were 5–7% for the calculated numbers of water molecules and of the adsorbed and free ions, 3–10% for mean square displacements of ions.

Results and discussion


Phosphorylation and the presence of salt can affect the conformations of cellulose molecules. To study the effects of phosphorylation and CaCl2 on cellulose conformations in the surface layer, the distributions of the dihedral angles, which define the main chain and phosphate group conformations, were measured (details in ESI, see Fig. S5–S8). The results show that phosphorylation affects the orientation of the glucopyranose cycle with the phosphate group only slightly and the differences between 12% and 25% phosphorylation are not significant. Moreover, the main chain conformation is virtually independent of the interactions of the phosphate groups with the salt ions (Fig. S6, ESI).

Density profiles

Some of the differences in the conformational changes between the systems with water and CaCl2 solution (Fig. S7 and S8, ESI) can be explained by the differences in the interactions of the two types of cations (Na+ and Ca2+) with phosphorylated surfaces. Fig. 2a shows the number density profiles of the different components. The number density profiles were calculated by dividing the space of the simulation cell into equally-sized slices of thickness 0.01 nm, and then inside each slice counting the number of atoms belonging to specific molecules or ions divided by the volume of the slice.
image file: c8cp05767b-f2.tif
Fig. 2 The number density profiles for cellulose (red line) and water (black and green lines), (a). The positions of the phosphorylated regions, (b). The arrow shows the direction normal to the surface. The position y = 2.0 nm chosen in (a) for ease of the presentation as the reference position, is clearly seen in (b). This position was measured from the position of the oxygen atoms of the primary hydroxyl groups of the first layer of the cellulose crystal.

As the figure shows, the degree of phosphorylation appeared to have no influence on the penetration depth. However, when CaCl2 is added water molecules attracted by counterions (interacting with different oxygen atoms of phosphate groups) penetrate deeper.

The Ca2+ and Na+ ions can interact attractively with four sites of the phosphate groups (OP1, OP2, OP3 and O6; Fig. 2b). To estimate the effect of the degree of substitution on the structure of the surface/ion interfaces and on the interactions of the Ca2+ and Na+ ions with the phosphorylated surfaces, the number density profiles were compared (Fig. 3). The changes are more pronounced for the surfaces in the CaCl2 solution (Fig. 3c and d).

image file: c8cp05767b-f3.tif
Fig. 3 The number density profiles for water molecules, Na+, Ca2+, Cl ions and oxygen atoms of the phosphate groups: for the surfaces with 12% and 25% substitution in water (a and b) and in CaCl2 solution (c and d). Illustration of cation's positions near the phosphate groups (e). The oxygen atoms of the phosphate group are shown by red, brown, olive and orange lines, water by dark blue lines, Na+ ions by magenta lines, Ca2+ ions by green lines and Cl ions by blue lines.

The bound surface charges of the phosphate groups are spread out over several Ångströms normal to the interface. The Na+ and Ca2+ ions are concentrated in two distinct layers. In contrast to the Na+ ions which penetrate deeper past the phosphate groups, most of the Ca2+ ions prefer to be in or near the layers of the oxygen atoms OP2 and OP3 (Fig. 3c and d). The total number of Na+ and Ca2+ ions attracted to the surface increases with increasing degree of phosphorylation. The attracted cations are immobilized (Fig. S9, ESI) due to their adsorption by the negative surface charges. The data in Table 1 shows that practically all cations that enter the interfacial layer (obtained from the number density profiles) are adsorbed in the surface region. This is determined by the number of cations in close contact with the phosphate groups (obtained from data shown in Fig. 4).

Table 1 The number of ions in the interfacial layer: N0 stands for the number of ions in the layer and NP for the number of ions per phosphate group. There are 16 and 32 phosphate groups for 12% and 25% phosphorylation, respectively
12% CaCl2 25% CaCl2 12% water 25% water
N 0 N P N 0 N P N 0 N P N 0 N P
Ca2+ 20.9 1.3 32.2 1.0
Na+ 2.8 0.2 13.2 0.4 9.0 0.6 22.8 0.7

image file: c8cp05767b-f4.tif
Fig. 4 The numbers of Na+, Ca2+ and Cl ions in the first coordination shell of P atoms in the system with pure water (a) and CaCl2 solution (b).

Before the discussion of the effect of the degree of phosphorylation on the number of the adsorbed cations demonstrated on Fig. 4, it is necessary to understand how the counterion layers are formed in our systems. To do that, let us compare the counter-ion density profiles with the predictions from the classical Stern theory.41

Comparison with Stern theory

In the Stern model, a layer of ions distributed from the charged surface to the bulk is divided into an inner layer (the layer of adsorbed ions – the Stern layer) and an outer layer (the diffuse layer – the Gouy layer). The Stern layer is separated from the Gouy layer at the so-called Stern plane at a distance that corresponds roughly to the radius of an ion. In our case for Na+ ions it is located at y = 2.38 nm. In the Gouy layer the Gouy–Chapman can be applied. The Gouy–Chapman relation is given as
image file: c8cp05767b-t1.tif(1)
where ρPB is counterion density, y is the distance from the surface, lB is the Bjerrum length (∼0.7 nm in pure water), σS is the surface charge density and μ is the Gouy–Chapman length,
image file: c8cp05767b-t2.tif(2)
where q is the elementary charge.

The bare surface charge density for the 25% phosphorylated surface is 0.4 e nm−2. However, taking into account the adsorbed ions neutralizing part of the surface charge, the effective surface charge density (σS) is 0.28 e nm−2. To calculate the renormalized σS we subtracted the charge of the adsorbed ions from the bare surface charge and the resulting charge was divided by the surface area. Similar procedure has been followed in other studies, see e.g.ref. 24.

It is known42 that if the surface charge density does not exceed 20 μC cm−2 (1.25 e nm−2) the Gouy–Chapman theory generally yields reliable results. This, however, depends on the valence of the ions interacting with the surface.43 Calculation of the density profiles for Na+ counterions in the system without salt by using eqn (1) yields low counterion densities in the interfacial regions (Fig. 5a).

image file: c8cp05767b-f5.tif
Fig. 5 Na+ counterion density profiles: directly from MD simulation (black line), calculated using eqn (1) with parameters determined from system properties (solid red line), and fitted by eqn (1) at varied lB and σS (dashed red line) (a). The dielectric constant calculated using eqn (3), (blue line) and a fit based on a sigmoid function eqn (4) (blue dashed line). The density profiles of cellulose and water are shown by red and green lines, correspondingly (b). Inset: The artificial region of dielectric constant profile and water density profile in this region.

To determine if the surface charge density or/and Bjerrum length become renormalized, we performed fitting using the Gouy–Chapman expression (eqn (1)) using σS and lB as free parameters. The result is shown in Fig. 5a. The fit yields σS = 0.25 e nm−2 (close to the σS value in the simulated system when accounting the ions in the Stern layer) and lB = 10 nm. Using lB and σS determined directly from the system without fitting yields very low ion densities in the interfacial regions and is in poor agreement with the MD data. The main reason for the above is the significant change in the dielectric constant as a function of distance from the surface.

Dielectric constant in the interfacial region

These changes can be characterized by using the Clausius–Mossotti relation and are illustrated in Fig. 5b for the system with 25% phosphorylation. The Clausius–Mossotti relation is given as
image file: c8cp05767b-t3.tif(3)
where ρ1(y) and ρ2(y) are the density profiles for cellulose and water, and ρ1,0 and ρ2,0 are the densities of pure cellulose and pure water in the simulated system. ε1 is the dielectric constant for pure cellulose (ε1 = 1.944), and ε2 is the dielectric constant for water (for TIP3P in CHARMM ε2 = 9239).

As Fig. 5b shows, in the region y = 2.9–3.2 nm the value of the calculated dielectric constant is higher than in bulk water (y > 3.2 nm) by about 10%. However, this increase is artificial and related only to the method of the calculation. As can be seen from eqn (3), the Clausius–Mossotti relation takes into account only changes in the densities. Therefore, the only reason for the increase in the value of the dielectric constant is the increase in water density in this region (see the inset in Fig. 5b) due to the high solvation of the phosphate groups. It is thus clear that interfacial regions cannot be treated entirely in terms of theories based on uniform density distributions.

The strong electric field due to the charged phosphate groups leads also to decreased mobility of water molecules in the region of the low dielectric constant.45 The obtained dependence (excluding the artificial region) fits very well to a sigmoid function, which is commonly used for the description of dielectric constant in the interfacial region46 (Fig. 5b).

image file: c8cp05767b-t4.tif(4)
In eqn (4), y0 (the position of the inflection point) and λ (the dielectric correlation length) were adjustable parameters, and ε1 and ε2 were defined by eqn (3). Fitting gives y0 = 2.72 nm and λ = 0.11 nm. These parameters define the range of the inhomogeneity in dielectric permittivity and the distance for recovering the bulk value of dielectric constant which are usually unknown.

Effect of phosphorylation degree on cation adsorption

The observed strong decrease in the dielectric constant in the interfacial layer allows us to explain the effect of the degree of phosphorylation on the number of adsorbed cations demonstrated in Fig. 4. The analysis shows that the number of the adsorbed Na+ ions per phosphate group increases whereas the number of Ca2+ ions decreases with increasing degree of phosphorylation. For Ca2+, this behaviour demonstrates that the level of the adsorbed ions is determined primarily by the surface area per binding site; increase in the number of binding sites results in shorter distances between them. Taking into account the much lower value of the dielectric constant in the interface region than in bulk solution (Fig. 5b), one can suggest that the charge–charge repulsion between the Ca2+ ions in the interfacial layer is the limiting factor for their adsorption. This leads to the observed decrease.

Why Na+ ions bind closest to the surface?

The question remains as why Na+ ions prefer to be in such a close vicinity to the surface (Fig. 3e)? As mentioned above, this is similar to what has been observed in lipid bilayer simulations of charged mitochondrial membranes25 albeit no detailed explanation was provided. Both systems have a complex surface structure and the dielectric constant changes rapidly from about 3 to about 80 (bulk water). Approaches using simpler Poisson–Boltzmann based theoretical models and implicit solvent Monte Carlo simulations22,23 are not able to capture such behavior: in contrast to the current observations and those from charged membranes,25 such approaches predict that the ions with a higher valence-to-volume ratio, Ca2+ ions, should be localized deeper in the interfacial layer (the valence-to-volume ratio for Ca2+ is higher as can be easily estimated using the CHARMM31 and OPLS47 parameters for ions).

The qualitative discrepancy between the detailed atomistic models and the simplified models shows that inclusion of surface structure and explicit solvation are crucial for accurate modelling and predictions; in simplified models charged surfaces are typically considered as homogeneous or localized charge distributions, or in the best cases as decorated surfaces (Fig. S11a, ESI).

During the simulations, the phosphate groups adopt conformations (Fig. S6 and S7, (ESI) show the preferred dihedral angles) which allow for the Na+ ions to form close contacts with three oxygen atoms: O6 and OP2 of the phosphate group, and O5 of the glucopyranose cycle (the peak for Na+ in Fig. 3; see also Fig. S10, ESI). This region is practically devoid of water molecules (see Fig. S10c and f, ESI) and thus, the main free energy losses are due to: (1) the removal of a part of the hydration shell of the Na+ ions and (2) the entropic losses associated with Na+ localization. These losses should be compensated by gains in entropy due to the removal of water molecules from the hydration shells and in electrostatic energy. The transfer of a water molecule from the immobile state to bulk gives an entropy gain of 8.4 kJ mol−1.48 However, the loss in energy per water molecule (taking into account coordination number for Na+ (4–8)49) is about 50 kJ mol−1.50 Na+ falls into the region with very low dielectric constant (Fig. 5b) and obtains maximum gain in electrostatic energy allowing for the energy and entropy losses to overlap. In contrast, in the region where Ca2+ ions are localized (Fig. 3), there must be additional losses due to dehydration of the phosphate groups. However, to the best of our knowledge, there is no data on the hydration energy of the phosphate groups in phosphorylated cellulose.

The data for a singly deprotonated phosphate ion can be used as an estimate. The data differs widely between studies51–53 but using even the smallest value52 gives the energy of ∼24 kJ mol−1 per water molecule. Thus, although dehydration provides a gain in entropy due to the release of water molecules from the phosphate groups, the energy loss is much greater. Moreover, in the region around a Ca2+ ion (Fig. 3), the dielectric constant is much larger and the gain in electrostatic energy is less than in the region around a Na+ ion. Thus, it is more favorable for the Na+ ions to be closer to the cellulose surface.

As for the Ca2+ ions, although the gain in electrostatic energy close to the surface (Fig. 3) would be greater, it appears that the loss in the hydration energy is too large – the hydration shell of Ca2+ is much larger than that of Na+. The large difference in the hydration energies of the Ca2+ and Na+ ions (−1579 kJ mol−1 and −406 kJ mol−1, respectively50) explains why Ca2+ ions prefer to be in the layer more distant from the surface.

Orientation of water molecules, electrostatic potentials and charge density profiles

Charged surfaces with attracted ions strongly influence the distribution of water molecules at the interface. Fig. 6a and b show the orientation of water molecules given as 〈cos(α)〉, where α is the angle between the water dipole and the surface normal as a function of the distance from the surface. In addition, charge density due to the explicit free charges (Na+, Ca2+, Cl) and the bound charges (partial atomic charges) on the phosphate groups and water molecules are shown. Fig. 6c and d show the electrostatic potentials obtained by integrating the corresponding charge distributions.
image file: c8cp05767b-f6.tif
Fig. 6 Orientation of water molecules given as 〈cos(α)〉, where α is the angle between the opposite to the water dipole and the surface normal, as a function of distance from the surfaces with 12% (dashed lines) and 25% (solid lines) phosphorylation in water (a) and in CaCl2 (b) (blue lines – charge density profiles). The electrostatic potentials due to Na+, Ca2+ and Cl, and the bound charges on the phosphate groups and on the water molecules for the surfaces with 12% (dashed lines) and 25% (solid lines) of phosphorylation in water (c) and in CaCl2 (d).

At short distances and without added salt, water molecules are mainly oriented such that the hydrogens point away from the surface (negative 〈cos(α)〉). Orientation of water dipoles is affected by the Na+ ions and the oxygen atoms of the phosphate groups. At larger distances, the water molecules are oriented with hydrogen atoms facing the surface as reflected by positive 〈cos(α)〉 values. This is due to the oxygen atoms of the phosphate groups. The considerable difference in the absolute values of 〈cos(α)〉 between these two layers is due to the differences in the amount of water encapsulated in them. In the more distant layer, the amount of water is much greater and hence the relative number of the oriented water molecules is smaller.

In the system with added salt, 〈cos(α)〉 is negative at short distances. The corresponding absolute value of the peak is higher than in water and the distribution becomes wider. The reason for the larger absolute value is that although the number of adsorbed Na+ ions in this region is less than in water (Fig. 3b and d), appearing here calcium ions increase the number of water molecules with similar (as for Na+) orientation relative to the surface. The reason for peak widening is the slight shift of the Na+ peak position to shorter distances (from 2.3–2.5 nm to 2.25–2.4 nm) and the appearance of a peak in Ca2+ ion adsorption at 2.45–2.6 nm.

The data in Fig. 6c and d demonstrates the electrostatic potentials in water and in salt solution. In water (Fig. 6c) the electrostatic potential of the 25% substituted surface is less negative (in comparison with 12% phosphorylation) at short distances (y = 2.5–2.6 nm) due to the larger amount of Na+ ions concentrated in this layer. At larger distances (y = 2.7–2.9 nm) this is due to the increased number of diffuse Na+ ions screening the negatively charged phosphate groups. In CaCl2 solution (Fig. 6d) the potential of the 25% substituted surface at short distances (y = 2.5–2.6 nm) is practically the same as for 12% due to the additional effect of the adsorbed Ca2+ ions. At distances 2.6–2.8 nm the electrostatic potential becomes more negative due to a lower (in comparison to 12% phosphorylation) number of Ca2+ ions per phosphate group.


We have performed explicit solvent all-atom MD simulations to investigate cellulose nanofibril phosphorylation and the effects of the resulting surface charge both with and without added CaCl2. Native cellulose was used as the reference system. In CaCl2 solution, the conformations of the grafted phosphate groups show strong dependence on the degree of phosphorylation, but in water (and in the presence of counterions only) the dependence is negligible. Importantly, the observed behaviour of mono- and divalent cations is qualitatively different from the predictions of both simplified Poisson–Boltzmann type analytical theories and implicit solvent simulations which predict that the ions with higher valency (Ca2+ here) should bind closest to the surface. Our analysis shows that hydration free energy and the surface structure of cellulose are responsible for the observed behaviour. Mean field theories and many simplified simulation setups neglect such effects. It was also observed that Na+ ions are localized close to the glucoside cycles both in water and in CaCl2 solution. In contrast, Ca2+ ions always remained concentrated further away.

Our analysis shows that this behaviour is due to competition between electrostatic interactions and hydration free energies of the ions. The results also have the important consequence that in mixtures of other monovalent cations, ions as large as Na+ ions will be adsorbed in the inner layer but larger ions, such as potassium, may not be integrated into this layer (1) due to steric hindrances if their sizes are too large or/and (2) when the energy gain from electrostatic interactions does not compensate the loss of hydration free energy. To obtain accurate quantitative predictions, two factors, namely, steric hindrances and electrostatics must be taken into account correctly. Not only the atomic level structural details of the surface, but also the ionic sizes are important. Besides the fact that the ions have different sizes in different force field, the ionic sizes are often rescaled to obtain better agreement with experimental data. Such adjustments have been made for protein, nucleic acid, and lipid systems (in particular, in CHARMM27).54 There is, however, no universal way for such adjustments.

To correctly account for electrostatics, not only ionic sizes but the partial charges are important. We used Mulliken population analysis for partial charge assignment, because Mulliken charges are generally used in CHARMM force fields.31 Different basis sets were tested and the only one that satisfied the requirement of preserving the crystal structure of the cellulose fibril after phosphorylation was the 6-31+G(d,p) basis set (see Table SII, ESI). It is, however, known that the Mulliken population analysis overestimates electronegativity of oxygen atoms and polarization effects.55,56 To check our results we used the Löwdin population analysis34 which, in contrast to Mulliken analysis, tends to underestimate atomic electronegativity suggesting smaller polarization effects.56 The Löwdin atomic partial charges satisfying the requirement of preserving the intrafibrillar hydrogen bonds were obtained with the basis set 6-31+G(d) (see Table SIII, ESI). Similarly to Mulliken charges, the simulations with Löwdin charges (Fig. S12, ESI, with density profiles) demonstrated that the mono- and divalent cations (Na+ and Ca2+) concentrate in separate interfacial layers with Na+ ions closer to the surface. Moreover, for both sets of the partial atomic charges, a decrease in the adsorbed number of Ca2+ ions and an increase in the number of adsorbed Na+ ions per phosphate group with increasing degree of phosphorylation were observed (Fig. S13a and b, ESI).

As a consequence of the differences in the Mulliken and Löwdin charges there is a difference in the absolute number of adsorbed Na+ ions and Ca2+ ions. Stronger influence on the adsorption of Na+ ions was observed and it is due to the lower electronegativity of the Löwdin charge on oxygen atoms of the inner part of the interfacial layer (Table SIII, ESI). The lesser observed effect for Ca2+ ions is due to a nearly proportional decrease of the electronegativity of the oxygen atoms and electropositivity of the phosphorus atom of phosphate groups (Table SIII, ESI). To obtain correct quantitative predictions especially when dealing with a system with divalent (Ca2+) cations, it might be beneficial to include polarization effects and charge transfer. By now the developed versions of the polarizable force fields (in particular, CHARMM Drude polarizable force field57) are applicable for proteins, DNA, lipids and for unmodified carbohydrates only.

Finally, our results demonstrate the existence of an optimum degree of phosphorylation in terms of adsorbed Ca2+ ions as an important criterion for materials design, for example, when choosing the cellulose matrix, e.g., for bone implants by biomimetic mineralization with a preliminary soaking of the matrix in a CaCl2 solution.4 Finally, we would like to point out that the current simulations were performed for only one salt concentration. Studies over a wide range of concentrations with varying degrees of phosphorylation are necessary to define optimal parameters and conditions for biomimetic mineralization. The conclusions section should come in this section at the end of the article, before the acknowledgements.

Conflicts of interest

There are no conflicts to declare.


This study was supported by the Russian Ministry of Education and Science within State Contract 14.W03.31.0014 (MegaGrant). Calculations were performed with the use of the computational resources of ComputeCanada/SharcNet and resources of the Institute of Macromolecular Compounds and the resources of Center of collective use “Complex of modeling and data research mega-class facilities” NRC “Kurchatov Institute”. Unique identifier RFMEFI62114X0006.

Notes and references

  1. N. Lin and A. Dufresne, Nanocellulose in biomedicine: Current status and future prospect, Eur. Polym. J., 2014, 59, 302–325 CrossRef CAS.
  2. T. Oshima, K. Kondo, K. Ohto, K. Inoue and Y. Baba, Preparation of phosphorylated bacterial cellulose as an adsorbent for metal ions, React. Funct. Polym., 2008, 68, 376–383 CrossRef CAS.
  3. T. Oshima, S. Taguchi, K. Ohe and Y. Baba, Phosphorylated bacterial cellulose for adsorption of proteins, Carbohydr. Polym., 2011, 83, 953–958 CrossRef CAS.
  4. Y. Z. Wan, Y. Huang, C. D. Yuan, S. Raman, Y. Zhu, H. J. Jiang, F. He and C. Gao, Biomimetic synthesis of hydroxyapatite/bacterial cellulose nanocomposites for biomedical applications, Mater. Sci. Eng., C, 2007, 27, 855–864 CrossRef CAS.
  5. H. G. de Oliveira Barud, R. R. da Silva, H. da Silva Barud, A. Tercjak, J. Gutierrez, W. R. Lustri, O. B. de Oliveira and S. J. L. Ribeiro, A multipurpose natural and renewable polymer in medical applications: bacterial cellulose, Carbohydr. Polym., 2016, 153, 406–420 CrossRef PubMed.
  6. W. Hu, S. Chen, J. Yang, Z. Li and H. Wang, Functionalized bacterial cellulose derivatives and nanocomposites, Carbohydr. Polym., 2014, 101, 1043–1060 CrossRef CAS PubMed.
  7. I. Sulaeva, U. Henniges, T. Rosenau and A. Potthast, Bacterial cellulose as a material for wound treatment: properties and modifications. A review, Biotechnol. Adv., 2015, 33, 1547–1571 CrossRef CAS PubMed.
  8. S. Saska, L. N. Teixeira, P. Tambasco de Oliveira, A. M. Minarelli Gaspar, S. J. Lima Ribeiro, Y. Messaddeq and R. Marchetto, Bacterial cellulose–collagen nanocomposite for bone tissue engineering, J. Mater. Chem., 2012, 22, 22102–22110 RSC.
  9. H. J. Jiang, Y. L. Wang, S. R. Jia, Y. Huang, F. He and Y. Z. Wan, Preparation and characterization of hydroxyapatite/bacterial cellulose nanocomposite scaffolds for bone tissue engineering, Key Eng. Mater., 2007, 330–332, 923–926 CAS.
  10. C. Gao, G. Y. Xiong, H. L. Luo, K. J. Ren, Y. Huang and Y. Z. Wan, Dynamic interaction between the growing Ca–P minerals and bacterial cellulose nanofibers during early biomineralization process, Cellulose, 2010, 17, 365–373 CrossRef CAS.
  11. M. Ghanadpour, F. Carosio, P. T. Larsson and L. Wågberg, Phosphorylated Cellulose Nanofibrils: A Renewable Nanomaterial for the Preparation of Intrinsically Flame-Retardant Materials, Biomacromolecules, 2015, 16, 3399–3410 CrossRef CAS PubMed.
  12. A. Naderi, T. Lindström, C. F. Weise, G. Flodberg, J. Sundström, K. Junel, J. Erlandsson and A. Runebjörk, Phosphorylated nanofibrillated cellulose: production and properties, Nord. Pulp Pap. Res. J., 2016, 31, 20–29 CAS.
  13. V. Kokol, M. Božič, R. Vogrinčič and A. P. Mathew, Characterisation and properties of homo- and heterogenously phosphorylated nanocellulose, Carbohydr. Polym., 2015, 125, 301–313 CrossRef CAS PubMed.
  14. S. R. Qiu and C. a Orme, Dynamics of biomineral formation at the near molecular level, Chem. Rev., 2008, 108, 4784–4822 CrossRef CAS PubMed.
  15. F. C. Meldrum and H. Cölfen, Controlling mineral morphologies and structures in biological and synthetic systems, Chem. Rev., 2008, 108, 4332–4432 CrossRef CAS.
  16. Y. Levin, Electrostatic correlations: from plasma to biology, Rep. Prog. Phys., 2002, 65, 1577–1632 CrossRef CAS.
  17. J. H. Harding, D. M. Duffy, M. L. Sushko, P. M. Rodger, D. Quigley and J. A. Elliott, Computational techniques at the organic-inorganic interface in biomineralization, Chem. Rev., 2008, 108, 4823–4854 CrossRef CAS PubMed.
  18. B. Grohe, S. Hug, A. Langdon, J. Jalkanen, K. A. Rogers, H. A. Goldberg, M. Karttunen and G. K. Hunter, Mimicking the biomolecular control of calcium oxalate monohydrate crystal growth: effect of contiguous glutamic acids, Langmuir, 2012, 28, 12182–12190 CrossRef CAS.
  19. A. D. Glova, S. G. Falkovich, S. V. Larin, D. A. Mezhenskaia, N. V. Lukasheva, V. M. Nazarychev, D. A. Tolmachev, A. A. Mercurieva, J. M. Kenny and S. V. Lyulin, Poly(lactic acid)-based nanocomposites filled with cellulose nanocrystals with modified surface: all-atom molecular dynamics simulations, Polym. Int., 2016, 65, 892–898 CrossRef CAS.
  20. O. Lyubimova, S. R. Stoyanov, S. Gusarov and A. Kovalenko, Electric Interfacial Layer of Modified Cellulose Nanocrystals in Aqueous Electrolyte Solution: Predictions by the Molecular Theory of Solvation, Langmuir, 2015, 31, 7106–7116 CrossRef CAS.
  21. A. Paajanen, Y. Sonavane, D. Ignasiak, J. A. Ketoja, T. Maloney and S. Paavilainen, Atomistic molecular dynamics simulations on the interaction of TEMPO-oxidized cellulose nanofibrils in water, Cellulose, 2016, 23, 3449–3462 CrossRef CAS.
  22. G. Tresset, Generalized Poisson-Fermi formalism for investigating size correlation effects with multiple ions, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 061506 CrossRef.
  23. M. Valiskó, D. Boda and D. Gillespie, Selective Adsorption of Ions with Different Diameter and Valence at Highly Charged Interfaces, J. Phys. Chem. C, 2007, 111, 15575–15585 CrossRef.
  24. A. A. Gurtovenko, M. Miettinen, M. Karttunen and I. Vattulainen, Effect of Monovalent Salt on Cationic Lipid Membranes As Revealed by Molecular Dynamics Simulations, J. Phys. Chem. B, 2005, 109, 21126–21134 CrossRef CAS.
  25. S. Pöyry, T. Róg, M. Karttunen and I. Vattulainen, Mitochondrial Membranes with Mono- and Divalent Salt: Changes Induced by Salt Ions on Structure and Dynamics, J. Phys. Chem. B, 2009, 113, 15513–15521 CrossRef PubMed.
  26. J. J. De Yoreo, P. U. P. A. Gilbert, N. A. J. M. Sommerdijk, R. L. Penn, S. Whitelam, D. Joester, H. Zhang, J. D. Rimer, A. Navrotsky, J. F. Banfield, A. F. Wallace, F. M. Michel, F. C. Meldrum, H. Cölfen and P. M. Dove, CRYSTAL GROWTH. Crystallization by particle attachment in synthetic, biogenic, and geologic environments, Science, 2015, 349, aaa6760 CrossRef PubMed.
  27. R. A. Böckmann, A. Hac, T. Heimburg and H. Grubmüller, Effect of sodium chloride on a lipid bilayer, Biophys. J., 2003, 85, 1647–1655 CrossRef.
  28. Y. Nishiyama, P. Langan and H. Chanzy, 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 CrossRef CAS.
  29. X. Xu, H. Wang, L. Jiang, X. Wang, S. A. Payne and J. Y. Zhu, Comparison between Cellulose Nanocrystal and Cellulose Nano fibril Reinforced Poly(ethylene oxide) Nano fibers and Their Novel Shish-Kebab-Like Crystalline Structures, Macromolecules, 2014, 47, 3409–3416 CrossRef CAS.
  30. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindah, Gromacs: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25 CrossRef.
  31. A. D. MacKerell, N. Banavali and N. Foloppe, Development and current status of the CHARMM force field for nucleic acids, Biopolymers, 2001, 56, 257–265 CrossRef.
  32. M. Kuttel, J. W. Brady and K. J. Naidoo, Carbohydrate solution simulations: producing a force field with experimentally consistent primary alcohol rotational frequencies and populations, J. Comput. Chem., 2002, 23, 1236–1243 CrossRef CAS PubMed.
  33. N. V. Lukasheva and D. A. Tolmachev, Cellulose Nanofibrils and Mechanism of their Mineralization in Biomimetic Synthesis of Hydroxyapatite/Native Bacterial Cellulose Nanocomposites: Molecular Dynamics Simulations, Langmuir, 2016, 32, 125–134 CrossRef CAS.
  34. P. O. Löwdin, On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals, J. Chem. Phys., 1950, 18, 365–375 CrossRef.
  35. S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys., 1984, 52, 255–268 CrossRef.
  36. W. G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697 CrossRef.
  37. M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys., 1981, 52, 7182–7190 CrossRef CAS.
  38. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys., 1995, 103, 8577–8593 CrossRef CAS.
  39. A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin and M. Karplus, All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins, J. Phys. Chem. B, 1998, 102, 3586–3616 CrossRef CAS.
  40. S. Miyamoto and P. A. Kollman, Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models, J. Comput. Chem., 1992, 13, 952–962 CrossRef CAS.
  41. O. Stern, Zur Theorie Der Elektrolytischen Doppelschicht, Ber. Bunsen-Ges., 1924, 30, 508–516 CAS.
  42. G. Bolt, Analysis of the validity of the Gouy-Chapman theory of the electric double layer, J. Colloid Sci., 1955, 10, 206–218 CrossRef CAS.
  43. A. G. Moreira and R. R. Netz, Simulations of counterions at charged plates, Eur. Phys. J. E: Soft Matter Biol. Phys., 2002, 8, 33–58 CrossRef CAS.
  44. S. A. M. Tofail, D. Haverty, K. T. Stanton and J. B. McMonagle, Structural order and dielectric behaviour of hydroxyapatite, Ferroelectrics, 2005, 319, 117–123 CrossRef.
  45. I. Danielewicz-Ferchmin, E. M. Banachowicz and A. R. Ferchmin, Dielectric saturation in water as quantitative measure of formation of well-defined hydration shells of ions at various temperatures and pressures. Vapor-liquid equilibrium case, J. Mol. Liq., 2013, 187, 157–164 CrossRef CAS.
  46. P. J. Lenart, A. Jusufi and A. Z. Panagiotopoulos, Effective potentials for 1[thin space (1/6-em)]:[thin space (1/6-em)]1 electrolyte solutions incorporating dielectric saturation and repulsive hydration, J. Chem. Phys., 2007, 126, 044509 CrossRef PubMed.
  47. W. L. Jorgensen and J. Tirado-Rives, The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin, J. Am. Chem. Soc., 1988, 110, 1657–1666 CrossRef CAS.
  48. J. D. Dunitz, The Entropic Cost of Bound Water in Crystals and Biomolecules, Science, 1994, 264, 670 CrossRef CAS.
  49. M. Patra and M. Karttunen, Systematic comparison of force fields for microscopic simulations of NaCl in aqueous solutions: diffusion, free energy of hydration, and structural properties, J. Comput. Chem., 2004, 25, 678–689 CrossRef CAS PubMed.
  50. D. W. Smith, Ionic hydration enthalpies, J. Chem. Educ., 1977, 54, 540 CrossRef CAS.
  51. P. George, R. J. Witonsky, M. Trachtman, C. Wu, W. Dorwart, L. Richman, W. Richman, F. Shurayh and B. Lentz, ‘Squiggle-H2O’. An enquiry into the importance of solvation effects in phosphate ester and anhydride reactions, Biochim. Biophys. Acta, Bioenerg., 1970, 223, 1–15 CrossRef CAS.
  52. J. Li, T. Zhu, G. D. Hawkins, P. Winget, D. A. Liotard, C. J. Cramer, D. G. Truhlar, Â. De Bordeaux, C. De Liberation, F. T. Cedex, H. Fock and B. Perdew, Regular article Extension of the platform of applicability of the SM5. 42R universal solvation model, Theor. Chem. Acc., 1999, 103, 9–63 Search PubMed.
  53. Y. Marcus, Thermodynamics of Solvation of Ions, J. Chem. Soc., Faraday Trans., 1993, 89, 713–718 RSC.
  54. J. Yoo and A. Aksimentiev, New tricks for old dogs: improving the accuracy of biomolecular force fields by pair-specific corrections to non-bonded interactions, Phys. Chem. Chem. Phys., 2018, 20, 8432–8449 RSC.
  55. J. Head, Partial Optimization of Adsorbates on Clusters: Oxygen on Al(111) Models for Surface and Bulk Phenomena, in Cluster Models for Surface and Bulk Phenomena, ed. G. Pacchioni, P. S. Bagus and F. Parmigiani, Springer, Boston, MA, 1992, pp. 415–422 Search PubMed.
  56. N. J. Mayhall, K. Raghavachari and H. P. Hratchian, ONIOM-based QM:QM electronic embedding method using Löwdin atomic charges: Energies and analytic gradients, J. Chem. Phys., 2010, 132, 114107 CrossRef PubMed.
  57. K. S. Vanommeslaeghe and A. D. MacKerell, CHARMM additive and polarizable force fields for biophysics and computer-aided drug design, Biochim. Biophys. Acta, 2015, 1850, 861–871 CrossRef CAS PubMed.


Electronic supplementary information (ESI) available. See DOI: 10.1039/c8cp05767b

This journal is © the Owner Societies 2019