2018 年 17 巻 1 号 p. 80-84
Solvation change induced by a chemical reaction can move a macromolecule in a specific direction. The hydrodynamic effect on a solvation motor was studied by using molecular dynamics simulations with a simple model. We examined the dependence of motor displacement induced by the chemical reaction on the size of the basic cell. The size dependence given by the simulation sat on a straight line as a function of the system size, and the displacement in an infinite system was estimated from the plot.
In a chemical reaction, solute molecules (reactants) transform into different molecules (products). In a solution, the solute–solvent molecule interaction changes in the reaction. For example, in an aqueous solution the solute–water interaction becomes larger after a charge separation reaction. The change is estimated to be 102–103 times the thermal energy at the standard temperature. This excess potential energy dissipates and a part of the energy transforms to kinetic energy of the solute molecule. If the charged part is not located in the center of the solute molecule, there is a possibility that the solute molecule moves in a specific direction.
Chemical-reaction-induced specific motion of a group in a protein is expected to be induced based on a similar mechanism. When a protein has an acceptor group for an electron, the group is negatively charged after the reaction. The electron-transfer reaction causes an increase in potential energy between the group and water immediately after the reaction. Then, the group should move in a specific direction.
In our previous study, molecular dynamics (MD) simulations were carried out to study the induced specific motion [1]. Here, a solvation motor is defined as a macromolecule whose motion is driven by solvation change. A simple system was adopted as a motor molecule and solvent molecules (Figure 1). In the simple system, the motor macromolecule had a reaction site on the surface. The induced specific motion appeared when the interaction change was larger than 103 times due to the change in solvation structure around the reaction site. The specific motion was confirmed under the NVE ensemble.
A solvation motor in a solvent. The diameters of the motor part (M), reaction site (R), and solvent (S) are 8.4, 2.8, and 2.8 Å, respectively. The distance between the centers of M and R is 3.4 Å. The molecular axis is defined as the R → M direction at the beginning of the reaction (t = 0) and is fixed thereafter. The surface of R sticks out from that of M by 0.6Å.
Examination of the hydrodynamic effect was inadequate in our previous study. The motion of the macromolecule induces a flow of solvent molecules and the range of the hydrodynamic effect is long [2,3,4,5]. In the previous study, a periodic boundary condition was adopted in the molecular simulation. It seems that the basic cell was smaller than the range of the hydrodynamic effect. A similar problem has been pointed out in the MD simulation to obtain the diffusion coefficient of a macromolecule in a solution [2,3,4,5]. There are several theoretical studies to obtain the diffusion coefficient of a macromolecule in explicit solvent molecules without simulation [6,7,8].
However, some studies have focused on the size effect with the MD simulation. The above difficulty in the calculation of the diffusion coefficient was pointed out and discussed by Fushiki [3]. In the study, a method for the calculation of the diffusion coefficient in an infinite system was proposed. In the method, some sizes of the basic cell were examined in the MD simulation. The diffusion coefficient of a macromolecule was obtained based on an extrapolation. Fushiki's method reveals that the hydrodynamic effect cannot be ignored in the MD simulation for a solution containing macromolecules. In our present study, we examined dependences of the induced specific motion on the size of the basic cell.
We examine a simple system composed of one motor molecule (MR) and solvent (S) particles. The motor molecule has a motor part (M) and a reaction site (R). The R − S and S − S interactions are described by the 6–12 Lennard-Jones (LJ) potential:
| (1) |
| (2) |
The size of S(σSS) was determined to be 2.8 Å
[10] to give a packing fraction of
ϕ = 0.383, which is almost the same as the value for water when
ρ = 1.0 g/cm3. The size of
M(σMM) is three times larger than that of
S, σMM = 3 σSS =
8.4 Å, and the size of R(σRR) is set as 2.8 Å.
We set the masses of these particles as mR =
mS = 3 × 10−26 kg and
mM = 26 mR. The summation of
mM and mR is equal to
33mS because we assume that the mass is
proportional to the volume, and the total volume of M and
R is almost equal to 33 times the solvent volume. The
interaction parameters of M, R, and S are
set as εMM = εRR =
εSS = 0.6517 kJ/mol =
1.082
| (3) |
In this study, only the R − S interaction is changed from εRS to ε΄RS = 103εRS. During the "chemical" reaction time, theR − S interaction parameter is ε΄RS . After the reaction time, ε΄RS returns to εRS. The S − S and S − M interactions do not change before, during, or after the reaction. The R − M intramolecular interaction was ignored in this model; thus, the parameter εRM was set as 0. The chemical reaction time is set to 100 fs, which leads to the largest displacement in our previous work [1]. To discuss the cell size dependence, the number of particles (M+R+S) is changed as N = 864, 2048, 4000, and 8788. The cell size was set so that the number density (weight density) of the solvent was ρ = 33.33 nm−3(1.0 g/cm3). Periodic boundary conditions were applied. The MD simulation period after the beginning of the reaction was 35 ps.
The MD simulations were performed using the velocity Verlet integration algorithm [12] with a time step of 0.5 fs. The M − R length was constrained using the RATTLE method [13]. The S − S interactions over 2.5 σSS = 7.0 Å were neglected. The simulations before the reaction were carried out with the Nosé–Hoover thermostat [14,15] at a constant temperature, T = 300 K. The initial conditions for the following processes were prepared in these NVT calculations. After the beginning of the reaction period, however, the simulations were carried out in NVE conditions. Error bars in Figures were obtained from the standard deviation, the number of trials, and multiplier (1.96) for the 95% confidence interval.
Fushiki gave a method to estimate the self-diffusion coefficient of a tagged molecule in a solution using linear-response theory and linearized hydrodynamics [3]. In the method, the diffusion coefficient in the infinite system is estimated by using the system size dependence of the diffusion coefficient. Here, a system size r0 is defined as follows:
| (4) |
Here, the solvation motor is almost spherical. The following discussion does not depend on whether it is a slip condition or a stick condition. Here, we assume the stick condition. The diffusion coefficient for a spherical macromolecule is:
| (5) |
Here, the average trajectory of the solvation motor MR as a function of time can be obtained by MD simulation. The average trajectory asymptotically approaches a value x∞ and it is defined as the displacement, i.e., the average trajectory x(t) at t = ∞. Because the average external force is zero after the reaction time, the equation of motion of the solvation motor MR is:
| (6) |
| (7) |
Substituting γ = v0m/x∞ into (5) yields the next relation:
| (8) |
This relation shows that the diffusion coefficient D is proportional to x∞. If we assume that the size dependences of T and v0 are smaller, we can expect a linear relationship between displacement x∞ and D. When there is a linear relationship, the displacement as a function of a'/r0 sits on a straight line.
MD simulations for the solvation motor MR were carried out and one of the average trajectories is shown in Figure 2. The average trajectory is the position in the direction along the molecular axis. The shape of the average trajectory as a function of time is similar to the solution of the equation of motion equation (7). The displacement is defined as the position at t = ∞. The displacements are estimated at 35 ps because the average motion of solvation motor MR is stopped by the solvent molecules.
Motion of solvation motor MR averaged over 5000 trajectories with N = 8788. Each motion is the position in the direction along the molecular axis. The average position is plotted as a function of time. The displacement x∞ is defined as the position at t = ∞.
Figure 3 shows the basic cell size dependence of displacement of a solvation motor MR as a function of system size a'/r0. The MD results seem to sit in a straight line. Thus, it seems that there is a linear relationship between the displacement and the diffusion coefficient of solvation motor MR. The displacement in the infinite system is estimated as 8.7 Å from the straight line at a'/r0 = 0. The calculated results in Figure 3 are not trivial because the order of displacement after the reaction is much larger than that in the diffusion behavior.
Basic cell size dependence of displacement of a solvation motor MR as a function of system size a'/r0. Circles are 5000 MD trials averages. The straight line is the least-squares fit.
Four system sizes, i.e., N = 864, 2048, 4000, and 8788 were examined. It seems reasonable because the diffusion coefficients calculated by MD as a function of a'/r0 sit on a straight line for the system sizes. The plot for the solvation motor MR is shown in Figure 4 and the diffusion coefficients calculated from the velocity autocorrelation function actually sit on a straight line. Therefore, we thought the system sizes are reasonable. Figure 4 also shows the importance of this cell size analysis in the calculation of the diffusion coefficient. The coefficient for the smallest system (864 solvent molecules) is about 70% of that for an infinite system. This difference in the diffusion coefficient is very large when we discuss the solvation effect based on equation (5).
The diffusion coefficient of the solvation motor MR along the molecular axis as a function of system size a'/r0. The R − S interactions do not change in this calculation. Circles are average MD simulation results of 5000 trials and the straight line is the least-squares fit of the results.
In the present calculation, the size dependence of displacement as a function of system size a'/r0 was a straight line. However, several dependencies of the displacement, such as packing fraction dependence, are not clear. In particular, the effects on the interaction change at the reaction part should be checked again because the interaction change increases the kinetic energy of the system. When the interaction change is kept, the average kinetic energy of a solvent molecule, which is almost proportional to temperature, becomes smaller as the system size becomes larger. Moreover, the diffusion coefficient is proportional to the system temperature. Therefore, there is a possibility that the reaction affects the diffusion and the displacement behaviors.
MD simulations for a solvation motor in an LJ solvent were carried out and the displacement along the molecular axis obtained to study the basic cell size dependence of the displacement. The displacement is plotted as a function of a'/r0. The size dependence sits on a straight line, which means that the displacement in the infinite system can be estimated from the straight line at a'/r0 = 0.
The authors dedicate this note to Professor Katsuyuki Nobusada who passed away on January 15th, 2018. This note was submitted on that day. The authors discussed this study in his room occasionally. He encouraged the authors at that time.The authors thank Prof. A. Yoshimori, Dr. Y. Nakamura, and Dr. Y. Uematsu for their comments. This work was supported by JSPS KAKENHI Grant Numbers JP17K14550, JP16H00774, JP16K05512, and JP15K05249. The computation was mainly carried out using the computer facilities at Reserch Institute for Information Techology, Kyushu University.
