Coarse-grained Simulation of Anionic Surfactant Sodium Dodecyl Sulfate

Through analyzing the deficiency of the current coarse-grained (CG) model, a new CG model for the ionic surfactant was proposed based on the Martini force field and iterative Boltzmann inversion method. In this model, the electrostatic interaction can be tackled by using a self-defined piecewise function to avoid the disadvantage of using coarse-grained solvents, and the VDW interaction parameters were derived by iterative methods. Using the improved model, the radial distribution function of NaCl and SDS solution in all-atom OPLS can be completely reproduced. The successful setup of the new coarse-grained model provides a good example of the construction of a high-precision coarse-grained force field.


Introduction
Molecular dynamics simulation is a powerful approach to studying the properties of lipid aggregates, self-assembly on a gas-liquid interface, surfactant interaction with new carbon materials and the interaction between surfactants and polymers. [1,2] However, many processes that occur in liquids, soft materials and biomolecular systems occur over length and time scales that are well beyond the current capabilities of atomic-level simulation, such as using molecular dynamics model (all-atom, AA and united-atom, UA).
As such, new and novel approaches continue to be developed that can access longer time and length scale phenomena. One such approach is coarse-grained (CG) simulation. [3] The post representative CG method is the Martini force field developed by the group of Marrink and Tieleman for many different lipids, [4][5][6][7] protein, [8][9][10] polymers [11][12][13] and sugar [14][15] and other more. [16][17][18] The Martini model is based on a four-to-one mapping, i.e. on average four heavy atoms plus associated hydrogens are represented by a single interaction center. This model aims for a broader range of applications without the need to reparameterize the model each time. A top-down approach is used by the extensive calibration of the nonbonded interactions of the chemical building blocks against experimental data, in particular, thermodynamic data such as oil/water partitioning coefficients. And this greatly improves the universality of the force field. Currently, the Martini force field is ideally suited to study the formation of micelles in solutions. [19][20][21][22][23][24][25] In a study comparing some surfactants, the Martini model reproduced experimental Critical micelle concentrations (CMCs) and aggregation numbers for non-ionic surfactants reasonably well. [26][27] Also, the group of Klein (W. Shinoda, R. Devane, and M. L Klein) established the new coarse-grained intermolecular force field (in short SDK) for a series of zwitterionic lipids. The model is designed to reproduce experimental surface/interfacial properties using simple function forms. Although there are many topologies of phospholipid molecules and surfactant molecules in the existing Martini force field and SDK force field, [28][29] they are not universality for the traditional typical surfactant systems, especially the ionic surfactant systems. For example, the Martini force field belongs to the semi-quantitative coarse-grained field, it is necessary to reiterate the parameters repeatedly to obtain a more reasonable simulation value. The process is more complicated.
In this paper, we will use take the most common anionic surfactant SDS (sodium dodecyl sulfate) as an example, based on the Martini force field to improve and establish a coarse-grained force field with higher accuracy. The effectiveness of the force field is tested by reproducing the basic physicochemical properties of SDS in the aqueous solution. It is expected that a set of coarsening force fields belonging to the surfactant system will be modified according to this method.

Component Van Der Waals Particle System
In the single-component Van der Waals particle systems, we used the Martini coarsegrained force field to perform molecular dynamics simulations on three simple particle systems, to analyze the relationship between the LJ potential energy function and the radial distribution function. The simulation was conducted with a pressure P of 1 bar and a temperature of 303K, and the particle mass (m) was set to 72 with the number of 400. The simulation time was 2 ns, and the constant σ was set 0.47 nm and ε were set 0.1 KJ/mol, 2 KJ/mol, and 4 KJ/mol, respectively. All the simulations were performed with Gromacs 5.0 [30] package. It can be seen from Figure.1 that the radial distribution function among the particles is single-peak in all systems because there are only one species of the particle in them. The height of the peak gradually decreases and the shape of the peak also changes from a sharp peak to a broad peak, when the interaction strength between the particles decreases. It should be noted that no matter how the intensity of the peak changes, the left side of the single peak always changes around the σ value (0.47nm). Specifically, the highest of the peak is closer the σ value when the intensity of the peak is low, whereas it is farther to the σ value. Therefore, we can adjust continuously the σ value and ε value between particles in the coarse-grained system, based on the radial distribution function between CG units in all-atom or united-atom force field simulation results. Finally, the fitting CG force field can be obtained in the single-component system.

LJ Potential Energy Function and Radial Distribution Function in Two-Component Van Der Waals Particle System
In the experiment, the single-component system was too simple and more multicomponent systems were used. Next, the relationship between the radial distribution function and potential energy function will be analyzed in the van der Waals system containing two kinds of particles. Assume that there are solute particles (S) and solvent particles (W) in the system, and the interaction parameters of solvent molecules are As generally seen in Figure 2a. and Figure 2c., the change of peaks is not obvious in appearance. It is worth noting that the main peak of RDF will be lower than the height of the latter peak, in the low polarity and high solubility system (seen in Figure 2d.). In summary, the RDF curve between solute molecules has the following rules: The main peak height is determined by the polarity and solubility of solute molecules. The solute with strong polarity ( large εss value) and low solubility (small εsw) has the highest peak value. While the intensity of small peaks relative to the main peak is mainly determined by the solubility, larger solubility results in smaller peaks. On the contrary, it is weak. Although the above rules are all derived with the same conditions of σ, when σ changes, the above rules are suitable even if the complexity of the peaks increases.

Radial Distribution Function of Charged Particles in Aqueous Solution
There is electrostatic force expect Van der Walls forces in the solution of charged particles. Using sodium chloride as an example, we first study the characteristics of the radial distribution function in a solution of a simple charge particle by the optimized potentials for liquid simulation all-atom (OPLS-AA) force field. The action parameters of chloride ion and sodium ion can be seen in Table 1. The water molecules is represented by SPC/E [31] (extended simple point charge) water model.
The simulation conditions were set as follows: The number of sodium ions, bromide ions and water molecules were 50, 50 and 2000, respectively. The simulation step was 5fs and the simulation time was 5ns. The temperature and pressure were the same as those of the previous two systems.
Addition Na + and Clin aqueous systems, Figure 3. shows that the relationship between the main peak of the RDF and the σ value no longer follows the rules of the first two systems. This is due to the fact that the electrostatic attraction between Na + and Clparticles is much greater than the Van der Waals forces. So it appears on RDF curve that Table 1. The σ and ε values of Na + and Cl -inOPLS-AA force field [32,33] particle σ(nm) ε(kJ/mol) the equilibrium distance between Na + and Clis much smaller than the σ value between them. Due to the large repulsive force between Na + and Na + , Cland Cl -, the intensity of the peaks on the RDF curve is very small, which is almost close to 1. Next, we only focus on the RDF between Na + an Clappearance. There are two peaks in Figure 3c. The first peak is very high and located exactly to the left of σNaCl, while the second is much weaker. This shows the electrostatic effect between Na + and Clis the strongest, when the distance between them is less than the minimum distance of Van der Waals action. However, when this distance is exceeded, the electrostatic interaction between Na + and Clrapidly decreases due to the hydration effect. What is worth noting is how the second small peak is formed? From the perspective of location, it is not formed through Van der Waals action.  stable, and this state corresponds exactly to the second peak in Figure 3c. So when the distance between Na + and Clcontinues to increase in NaCl aqueous, both Na + and Clare in a fully hydrated state, the electrostatic force between them is relatively weak.

Discussion
The Coarse-grained simulation of aqueous NaCl solution In the previous coarse-grained model, simple ions (such as Na + , Cl -,etc.) are often coarse-grained into one particle along with the hydrated water molecules. In the simulation process, the water molecules considered as the largest consumption of calculated resources, coarsely granulates multiple water molecules into one particle through the "four-in-one" or "three-in-one" rule. Because the coarse-grained water molecules have no polarity, they can not effectively shield the charge of charged particles like the water molecules in AA or UA force field. In these coarse-grained force fields, therefore, the relative permittivity is modified to weaken the electrostatic force between the charged particles. By analyzing the hydration state of Na + and Clin water, it can be concluded that the relative permittivity of aqueous solution will change with the distance between ions, which can be expressed by a piecewise function: In the above equation, r0 represents the radius of hydration of the ions. When the distance r between ions is smaller than r0, the relative permittivity εr can be expressed as an exponential function. when r is greater than r0, εr is expressed as a constant 78 (the exponential function of water). Since most of ions have a hydration radius of 0.3-0.4nm, the r value is uniformly taken as 0.4nm in our model. Table 2 lists the parameter in our coarse-grained force field. The non-bond parameters σ and ε of Na + and Clalso come from the OPLS field. The model of water directly uses the "four-in-one" coarse-grained model in the Martini force field. Therefore, it is only necessary to adjust the interaction about Na + (Cl -) with W (the water's CG particle), respectively. It is worth noting that we use the LJ(12-3) function, but not LJ  function, about the non-bonding interaction between Na + (Cl -) and water coarse-grained particles W. This is because electrostatic interactions exist between charged ions and water molecules, and the LJ(12-3) function can more accurately reflect the interaction between the two. It can be seen from Figure 5., our simulation results can reproduce RDF curve obtained using the OPLS force field. In Figure 5c, it should be pointed out that the second peak of the RDF curve can not be reproduced using the coarse-grained force field. It is due to the absence of the polarity in the coarse-grained water molecule model and the inability to truly reflect the hydration conditions about anions and cations.

The Coarse-grained of Sodium Dodecyl Sulfate Aqueous Solution
A set of simple rules can be summarized to adjust the interaction parameters of each coarse-grained particle, by analyzing the interaction between the RDF of particles in simple particle system and LJ potential energy and electrostatic potential energy. The most common anionic surfactant, sodium dodecyl sulfate (SDS), is used as an example to illustrate the adjustment process of parameters of the coarse-grained field of ionic surfactant. shown in Figure 8. The head of the SDS molecule consists of one S atom and four O atom, and they are divided into a coarse-grained particle with a unit of negative charge(SO4 shown in Figure 7.). For the carbon atoms in the hydrophobic chain, the principle of "three-in-one" is adopted. That is three carbon atoms together with hydrogen atoms are used as a coarse-grained unit, and the hydrophobic tail is divided into four coarse-grained particles(C1, C2, C3and C4 shown in Figure 7.). The four CG particles are classified into two types. The CG particle of carbon attached to the head group(C1) are CE type ,the other three CG particle(C2, C3 and C4) are CT types. In other CG models, inorganic salt ions (such as Na + and Cletc), together with hydrated water molecules, are often used as a CG unit. In our CG model, Na + is seen as an independent CG particle(NA shown in Figure 8.) with a positive unit charge. The advantage of this is that the simulation results can truly reflect the role Na + in the self-assembly of SDS molecules and increase the accuracy of the simulation.

The Bonded parameters in SDS Coarse-grained Model
In the SDS molecule, the potential energy between CG particles is in the form of a harmonic potential function： Equation (2) and (3)  respectively. The above parameters can be obtained by comparing with the results of using a fine force field to simulate SDS molecule, see Table 3. From the parameters of the intramolecular action, it can be seen that both the stretching vibration and the bending vibration of the bond angle are much large than the standard Martini force field. This is because the chain length of surfactants is shorter than that of large biomolecules, and therefore their softness is much lower.

Non-Bonded Parameters in the SDS Coarse-grained model
In the first half, we have summarized the relationship between the potential energy function and the RDF in a system where only van der Waals forces exist and a charged particle system. On this basis, we will try to estimate the potential energy parameters based on the RDF of partiles in in the CG system containing chemical bonds, thus, complete the CG process of the complex system. Unlike simple particle systems, the existence of chemical bonds leads to a stronger correlation between particles, and the RDF between them is more complicated. The SDS aqueous solution system is used as an example to describe the process of obtaining the model parameter in the ionic surfactant CG force field.
First, we obtain the RDF curve of SDS molecules in aqueous solution using the fine force field (OPLS-UA force field). The steps are as follows: In a cubic box with side length of 5.5nm, 10 SDS molecules were randomly placed, and then 4502 water molecules were added. Since the force field of the SDS molecules has been reported many times, there is no further verification of the reliability of the of force field. The water molecules adopts the SPC/E model. The simulations were performed in the NpT ensemble as taken in previous paper [34,35] , the T (temperature) was kept at 303K by applying a Berendsen thermostat with the time constant to 0.1ps, and the P (pressure) was held constant at 1 bar using a Parrinello-Rahman barostat,with the time constant to 0.5ps.
After the above setting, the simulation step can be set to 5fs. The simulation was performed in the molecular dynamics package Gromacs 5.0 with a simulation time of 20 ns. The last 10 ns of the trajectory was used to analyze the RDF between CG cells, as shown in Figure 8.
In the CG simulation system, the number of molecules such as SDS and water, and the simulations conditions such as temperature and pressure were all set exactly the same as those in the above system. The treatment for Van der Walls forces is the shift method.
The values of rvdw_switch and rvdw are set as 0.9nm and 1.2nm, respectively. The electrostatic force still use a segmented approach, as shown in equation (1). Although the radius of SO4 is large than Na + or Cl -, its hydration radius is not much different, so the parameters in the piecewise function of SO4 are exactly the same as those of Na + or Cl -. Compared with the molecular model of the atomic level, the simulation step size can be modified to 20fs due to the greatly increased distance between CG units. In dealing with the intermolecular potential energy, we dealt with Van der Waals forces and electrostatic forces separately. Therefore, the RDF curves of CG system and the atomic (OPLS-UA force field) system can be realized the coincidence by repeatedly adjusting the Van der Walls potential energy parameters between the CG particles and comparing the RDF images of the CG system and the atomic system, the results shown in FIgure 8. Figure 9 shows the snapshots of SDS in atomistic (a) and the coarse grained (g) simulation, respectively. In essence, our CG method is consist with Boltzmann inverse iterative method. In the process of adjusting coarse granulation parameters, the following matters should be paid more attention: First, the σ value of the Van der Walls interaction between the CG units needs to be determined, because the σ value has the greatest influence on the RDF image between the two CG units. As mentioned earlier, the σ value between most CG particles is on the left side of the main RDF peak, except for the σ value between the CG particles with positive and negative charges. So the σ value can e directly estimated from the RDF curve.
The second, the intensity of CG RDF curve can generally be adjusted by the value of ε, the higher the peak intensity, the larger ε value.
The third, when fitting the potential energy parameters between charged CG units such as SO4-SO4, and SO4-NA, we must consider the adjustment of the RDF among other particles. This process is actually a matter of balancing electrostatic potential energy and Van der Waals potential energy. After adjusting, potential energy relationships between CG units in the SDS molecule can be obtained, as shown in Table 4.

Conclusion
In this article, we propose a new CG model for ionic surfactants by improving the electrostatic treatment and the fitting method of Van der Walls in the original CG model.
The model can well reproduce the RDF curves of NaCl and SDS aqueous solution system.
This initially realizes the establishment of a high-precision CG force field, and lays a good foundation to accurately simulate the large scale surfactants.

Methods
Like the traditional molecular force field, the common CG force field also uses the same potential energy expression form: That is, the energy of the entire molecular system ( mol E ) is divided into the intramolecular potential energy( intra E ) and the intermolecular potential energy ( inter E ).
The parameters of intra E can be obtained by comparing the corresponding CG units with the united-atom or the all-atom model, and repeatedly adjusting the stretching and the bending vibration constants in the CG model. The treatment of inter E will directly affect the precision and universality of the CG force field, compared with intra E . Therefore, the key task in developing the CG force field is to fit the non-bonded potential energy For the van der Waals forces between different coarse-grained units, the most commonly Lennard-Jones12-6 [30] potential energy functions form was used (see in equation 6).
Wherre ij  represents the minimum distance between interaction particles and ij  represents the strength of the interaction between them as shown in Figure 10.