Effect of the Simulation Box Size and Precipitant Concentration on the Behavior of Tetragonal Lysozyme Dimer

The 10-nanosecond simulation of a lysozyme dimer, which is a fragment of the tetragonal lysozyme crystal structure, has been carried out by the molecular dynamics method at different simulation box sizes and precipitant concentrations in a solution. The dimer stability has been estimated by calculating the root-mean-square fluctuations of protein atoms. It is shown that the box size does not significantly affect the mobility of protein atoms on a relatively short trajectory, while the effect of the precipitant concentration on this trajectory is noticeable.


INTRODUCTION
Study of protein crystallization and the properties of protein crystals is undoubtedly interesting both for the development of approaches to search for protein crystallization conditions and investigate the difference between the structures of a protein molecule in the crystalline state and in a solution. It was found for several proteins that the pre-crystallization solution contains oligomers built from protein molecules, the structure of which is well-described by the structures of oligomers found in crystals of the examined proteins [1][2][3].
The molecular dynamics (MD) method is widely used to explore the crystallization of proteins and the properties of protein crystals. The results of the MD simulation of protein crystal properties were reported in [4][5][6][7][8][9]. In addition, the crystal formation is studied applying various computational techniques. In particular, the dependences of the protein crystal stability on the ion concentration in the crystallization solution and on the charges of amino-acid residues on the protein surface were investigated in [10]. In [11], an attempt was made to simulate the formation of a protein crystal. The MD method was used also to study the properties of protein oligomers present in the precrystallization solution [12,13].
It is noteworthy that numerous studies on the protein crystallization using the computational biology methods are also based on the analysis of the mobility of atoms of protein molecule under different conditions [10,12,13]. It was shown that the stability of a biological tetrameric molecule of human hemoglobin in the MD simulation in a solution is affected, along with other factors, by the simulation box size [14]. It is unknown, however, to what extent this fact [14] (which is indicative of the effect of chosen simulation box parameters on the protein oligomer stability) can be taken into account in an MD simulation of oligomers formed in a crystallization solution of proteins. To answer this question, we used the MD method to study a lysozyme dimer, which is a fragment of the tetragonal lysozyme crystal structure [15], in two different simulation boxes and at three different precipitant concentrations in a solution in order to establish to what extent these parameters affect the stability of protein dimers.

EXPERIMENTAL
Preparation of dimer model. The molecular model of a possible crystal growth unit was built using the data on the tetragonal lysozyme crystal structure (PDB_ID: 6QWY) solved with a resolution of 1.35 Å [15]. The crystal belongs to the sp. gr. P4 3 2 1 2 with the unit-cell parameters a = b = 77.147 Å, c = 37.139 Å, and α = β = γ = 90.00°. The independent part of the unit cell contains one lysozyme monomer.
The simulation was performed using periodic boundary conditions. The protein dimer was placed at the center of a cubic box. The simulation box size was set such that the distance between any protein atom and the box edge was at least 1 nm in one experiment and at least 2.5 nm in the other; the box edge lengths were 8.6 and 11.6 nm, respectively. The interaction potentials were taken into account only for the atoms localized within a radius of 1 nm.
To fill the simulation box with an explicit solvent, a 4-site TIP4P-Ew water model was chosen, in which the Ewald summation methods can be applied [20]. To reproduce more accurately the experimental conditions of protein crystallization, lysozyme dimers were simulated in an aqueous solution with a precipitant. To this end, some of water molecules were replaced with sodium and chlorine ions so that the salt concentration was 0.4 M (25 mg/ml) in the cubic box with an edge of 11.6 nm and 0.2, 0.4, and 0.6 M in the box with an edge of 8.6 nm. The total charge of the system was neutralized by adding a small number of chlorine ions.
Before the beginning of each MD computation, the energy of the system was minimized by the steepest descent method (50 000 steps) until the force acting on any atom became weaker than 1000 kJ/(M nm). Then, the system was thermostated for 100 ps by the advanced Berendsen method [21] and barostated by the Parrinello-Rahman method [22] (100 ps); in this case, the limiting forces (k pr = 1000 kJ/(M nm 2 )) were applied to all nonhydrogen protein atoms and all bonds were constrained using the LINCS algorithms [23].
The MD simulation was performed in an isothermal-isobaric ensemble using a Berendsen thermostat and a Parrinello-Rahman barostat (at P = 1 atm for all the systems). The integration was made using a standard leap-frog integrator [24], the integration time step was 2 fs, and the atomic coordinates were saved every 1000 ps. The electrostatic interactions were calculated by the particle mesh Ewald (PME) summation [25] with the fourth-order cubic interpolation and a Fourier grid spacing of 0.16 nm. The duration of each calculated trajectory of protein dimers was 10 ns.

RESULTS AND DISCUSSION
Basing on the MD simulation data, after structural alignment of trajectories to the initial position, the root-mean-square fluctuations (RMSFs) of atoms were calculated and plotted, which are the fluctuations of atomic coordinates averaged over the entire simulation time (10 ns). It is convenient to use these values to estimate the mobility of protein atoms under different conditions. Figure 1 shows the RMSFs for С α atoms of lysozyme dimer in the solution at a precipitant concentration of 0.4 M for different simulation box sizes (the box edges are 2.5 and 1.0 nm). It can be seen that, on a relatively short trajectory of 10 ns, the box size does not significantly affect the mobility of protein atoms, although there is some difference between the mobilities of atoms in certain protein parts. Note that the residues with a difference of more than 0.0275 nm between the RMSF values are SER81, ALA82, GLY102, TRP111, ARG112, ASN113, GLY117, and GLY126 in monomer 1 and HIP15 (HIS in the protonated form), PRO70, GLY71, THR89, ALA90, ARG114, LYS116, GLY117, THR118, ASP119, and ALA122 in monomer 2. Figure 2 shows the RMSFs for С α atoms of lysozyme dimer in a solution at precipitant concentrations of 0.2, 0.4, and 0.6 M. It can be seen that the effect of precipitant concentration on the mobility of protein atoms along the same trajectory is most pronounced. In particular, for ASN37, ALA42, ASN44, ARG45, THR47, ARG68, THR69, PRO70, ALA107, ALA122, TRP123, TRP124, GLY126, and LEU129 in monomer 1 and for GLY67, ARG68, and THR69 in monomer 2, the difference between the RMSF values for the same residues but at different precipitant concentrations is more than 0.05 nm. For PRO70 in monomer 1, the difference between the RMSF values of C α atoms is more than 0.1 nm for the systems with NaCl concentrations of 0.2 and 0.6 M.
The results obtained show that, in the MD estimation of the stability of lysozyme dimers, the precipitant concentration in the solution is more important than the simulation box size; even on a relatively short (10 ns) trajectory, the mobility of the С α protein atoms changes by more than 30%, depending on the precipitant concentration, while the effect described in [14] is hardly present in the investigated systems. Note that, with an increase in the number of sodium chloride ions in the range from 0.2 to 0.6 M, the mobility of С α atoms of monomer 1 gradually decreases, while the dependence of the mobility of С α atoms of mono-mer 2 on the precipitant concentration is nonmonotonic.

ACKNOWLEDGMENTS
This work has been carried out using computing resources of the federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at NRC "Kurchatov Institute," http://ckp.nrcki.ru/.

FUNDING
This study was supported in part by the Ministry of Science and Higher Education of the Russian Federation within a state assignment for the Federal Scientific Research Center "Crystallography and Photonics" of the Russian Academy of Sciences; the Russian Foundation for Basic Research, project no. 19-29-12042 MK; and the National Research Centre "Kurchatov Institute," order no. 1360.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.