Free Energy Change during the Formation of Crystalline Contact between Lysozyme Monomers under Different Physical and Chemical Conditions

Abstract

We use the MM/GBSA method to calculate the free energies of dimer formation by binding two monomers with different combinations of precipitant ions, both embedded in the structure of monomers and in the crystallization solution. It shows that the largest difference in free energy values corresponds to the most accurate dimer model, which considers all precipitant ions in their structure. In addition, it shows that in the absence of precipitant ions in the solution of lysozyme molecules, a monomer is a more energetically favorable state.

1. Introduction

Despite significant progress made in recent years, protein crystallization remains a highly laborious, more empirical, and the least predicted stage of structural research, guided by the results of previous attempts.

Due to its availability and ease of its crystallization, the largest number of studies on the growth of protein crystals and the effect on the conditions of the crystal structures are devoted to hen egg white lysozyme (HEWL). It makes it possible to obtain sufficiently large (up to 2 mm [1]) lysozyme crystals in a wide range of conditions, and use them as model objects for studying the mechanism of protein crystallization [2]. According to [3], lysozyme crystals are almost all syngonies: tetragonal, monoclinic, orthorhombic, hexagonal, and triclinic. This study considers lysozyme dimers that are part of a tetragonal syngony crystal.

For many years, there have been various data on the content of the pre-crystallization solution of lysozyme. Some researchers assumed that it consists only of monomers [4], but most agreed that oligomers are formed in solution. In recent works [5,6,7], it was found that, before its crystallization in solution, lysozyme’s oligomers are formed, which are involved in the subsequent growth of the crystal. The authors of this study proposed to use the molecular dynamics (MD) method to assess the stability of oligomers and to study the behavior of atoms and bonds in their structure [8,9]. This technique consists of the MD modeling of isolated from the protein crystal lattice oligomers, and the subsequent analysis of the mobility of their atoms. Using the proposed method, the type of octamer formed in the pre-crystallization solution of lysozyme [8] has already been identified, and the influence of the precipitant ions built into the structure of the protein on the stability of monomer, dimer, and octamer has been determined [9]. However, it has not yet been clarified whether the change in free energy when binding two monomers correlates with the stability of dimers estimated by root mean square fluctuations (RMSF) graphs.

To calculate the change in free energy during the formation of a dimer, the molecular mechanics/generalized Born surface area (MM/GBSA) method was chosen. It was first used in [10], and was initially based on the solution of the Poisson–Boltzmann equation, but to speed up the calculations, it was optimized using the generalized Born equation [11,12]. The purpose of this approach is to determine the difference between the free energy values of solvated molecules in bound and free states. The change in free energy upon the binding of the receptor to the ligand (ΔG) can be calculated based on the following equation:

ΔG = ΔG_vac + ΔG_complex − (ΔG_lig + ΔG_rec),

(1)

where ΔG_vac is the change in free energy when a ligand binds to a receptor in a vacuum, ΔG_complex, ΔG_lig, and ΔG_rec represent the difference in free energy values when in a vacuum and in a solution for the complex, ligand, and receptor, respectively.

The solution of the generalized Born equation [12,13,14] for each of the three states (ΔG_complex, ΔG_lig, and ΔG_rec) determines the electrostatic component (ΔG_el) and, when adding an empirical term for hydrophobic interactions (ΔG_hydr), allows us to calculate the energy of solvation:

ΔG_solv = ΔG_el + ΔG_hydr,

(2)

where ΔG_vac is found using molecular mechanics by averaging the energies of the interaction of the ligand with the receptor obtained from an ensemble of MD trajectory snapshots (ΔG_mm) and considering the change in entropy (−TΔS):

ΔG_vac = ΔG_mm − TΔS.

(3)

However, if only calculating the difference in energies between states with almost the same entropy, then this contribution is ignored.

By this method, we calculated the free energies of dimer formation by binding two monomers with different combinations of precipitant ions, both embedded in their structure and in the crystallization solution. We found that the most noticeable decrease in free energy occurs at the highest considered concentration of the precipitant in the solution.

2. Materials and Methods

2.1. Construction of the Initial Models of Oligomers

We investigated the same 3 types of dimer models as in [9]: with Na and Cl ions associated with the protein molecule (3 Na ions and 4 Cl ions per monomer), with lysozyme-bound Na ions (one ion per monomer) and without precipitant ions incorporated into the dimer structure. The protonation states of amino acid residues at pH 4.5 (in accordance with the experimental conditions [15]) were assigned by the PROPKA server [16].

2.2. Molecular Dynamics Simulations

All calculations and preparations of structures were performed in the GROMACS software package version 5.0.4 [17] and were carried out in the same way for different dimer models, if it is not indicated. Amber ff99SB-ILDN [18] was chosen as the force field, since the new potential parameters for some torsion angles were added to it.

During the system equilibration and productive MD, three-dimensional periodic boundary conditions were applied, and long-range electrostatic forces were treated with the smooth particle mesh Ewald summation method (PME [19]) with cubic interpolation and a Fourier grid spacing of 0.16 nm. The short-range interaction cutoff was set at 1 nm. All bonds were constrained by the LINCS algorithm [20].

The prepared dimer structures were placed in the center of the cubic simulation box, with the minimum distance between its edge and the protein molecule set at 1 nm. The remaining box space was filled with TIP4P-Ew water [21] to perform MD calculations using the Ewald techniques. When simulating a dimer in a solution with a precipitant, water molecules were replaced by Na and Cl ions so that the concentration of NaCl in the box was 0.4 M. To neutralize the net charge of each system, a small number of Cl ions was added.

Energy minimization and equilibration of the systems were performed according to the following protocol. First, the energy was minimized by the steepest descent algorithm (50,000 steps), using a force constant of 1000 kJ M⁻¹nm⁻² for position restraints. Then, the boxes were sequentially equilibrated in NVT- and NPT-ensembles by the modified Berendsen thermostat (V-rescale) [22] and the Parrinello–Rahman method [23], respectively (for 100 ps each). The integration step was set at 2 fs, the temperature at 283 K, and the pressure at 1 atm.

The productive MD was calculated in the NPT ensemble using the modified Berendsen thermostat and the Parrinello–Rahman barostat. The equations of motion were integrated using the standard leap-frog algorithm [24]. The duration of each simulated trajectory of lysozyme dimers was 100 ns.

Before analyzing the trajectories, artifacts arising from quasi-infinite periodic boundary conditions were eliminated using the gmx trjconv command with the -pbc nojump flag, returning the dimer back to the simulation box.

2.3. Calculations of Binding Free Energy

The free energy of binding two monomers was computed by the molecular mechanics/generalized Born (Poisson–Boltzmann) surface area (MM/GB(PB)SA) method [9], which combines molecular mechanics and the solution of the Poisson–Boltzmann equation. The applied module gmx_MMPBSA [25] version 1.4.0, in combination with the MMPBSA.py script [26] and the AmberTools20 package [27], made it possible to use the MM/GB(PB)SA method with the generalized Born model [28] when processing the MD simulation results carried out using the GROMACS program. The calculation of the free energy for each trajectory was performed on 1000 frames.

3. Results

Based on the results of the MD simulation, the graphs of the mean square fluctuations of C_α atoms (RMSF, Figure 1) were plotted, from the comparison of which it is noticeable that curve 1, corresponding to the dimer model considering all precipitant ions, lies below the rest. This indicates the relative stability of this model. In Figure 1, the biggest difference in RMSF values (in the range from 0.15 to 0.27 nm) between different models is observed in the region of amino acid residues GLY126, CYS127, ARG128, LEU129, and GLY102 of the chain A. We note that these residues are located on the surface of the protein molecule, but do not participate in the formation of covalent or hydrogen bonds between monomers; the residue LEU129 is the terminal in the polypeptide chain, and GLY126, CYS127, and ARG128 are placed next to it. The RMSF value of the dimer averaged over all C_α atoms, modeled without Na and Cl ions both in the solution and in the protein structure, is 0.16 nm. When a precipitant is presented only in the solution (0.4 M), it is 0.15 nm; when the precipitant is in the solution and the built-in ions of Na or built-in ions of both Na and Cl are considered, it is 0.14 and 0.12 nm, respectively. Consequently, the more precipitant ions are accounted for in the simulated system, the more stable the dynamics of the dimer in the solution and less flexible its structure, since the RMSF values of its atoms decrease.

The maximum RMSF value in Figure 1 is 0.44 nm (for the CYS127 of chain A of the dimer with embedded Na ions), and this indicates that the studied trajectories of dimers are sufficiently stable to apply the MM/GBSA method to them to determine the energy change upon binding of monomers. Figure 2 shows the root mean square deviations (RMSDs) for the C_α atoms along the 100 ns trajectories of all systems under study. The RMSD values are in the range 0.09–0.43 nm, which indicates the relative stability of the dimers during the simulation. In Figure 3, the radii of gyration for all systems practically remain the same during the simulation process, in the range from 1.95 to 2.15 nm. It demonstrates that the compactness of the dimers did not change significantly during the MD simulation.

The MM/GBSA method was used to calculate the change in free energy upon the binding of monomers in solutions with different NaCl concentrations near the value corresponding to optimal lysozyme crystallization conditions (0.4 M) and near the point of 0 M as no dimers are observed in this case (

Table 1. Comparison of changes in free energy upon the binding of lysozyme monomers with different combinations of precipitant ions bound to the protein (Na and Cl, only Na, without Na and Cl) in solutions with the addition of NaCl at various concentrations. The standard deviation for ΔG is provided (the standard deviation of the mean value is a factor of n^1/2 lower, where n is the number of snapshots and equals 1000).

The Concentration of Precipitant in the Solution, M	Na and Cl, kcal/M	Na, kcal/M	Without Na and Cl, kcal/M
0.6	−7.85 +/− 6.1	−7.97 +/− 4.7	−7.61 +/− 4.9
0.4	−7.78 +/− 6.1	−7.87 +/− 4.7	−7.49 +/− 4.9
0.2	−7.43 +/− 6	−7.43 +/− 4.7	−7.08 +/− 4.8
0.1	−6.71 +/− 5.9	−6.56 +/− 4.7	−6.29 +/− 4.6
0.05	−5.58 +/− 5.7	−5.22 +/− 4.8	−5.1 +/− 4.4
0.01	−1.97 +/− 5.4	−1.11 +/− 4.8	−1.36 +/− 4
0	5.82 +/− 5.3	7.11 +/− 4.8	5.55 +/− 5

).

There were two types of precipitant ions studied:

1. Na and Cl ions in the aqueous solution surrounding dimer (their concentration is given in column 1 of Table 1). They were added to the simulation box by replacing water molecules;

2. Na and Cl ions which were found to be associated with the protein in the lysozyme crystal. We investigated three combinations of such ions: when both of them are bound with the monomers, when only Na ions are embedded in the protein structure, and without both Na and Cl ions incorporated into the lysozyme structure (columns 2, 3, and 4 of Table 1, respectively). These precipitant conditions were chosen because, as in our previous work [9], it was found that MD results on lysozyme oligomers stability investigations are consistent with the SAXS experiments [5,6] only when accounting for all precipitant ions associated with the protein.

Table 1 shows that the formation of a dimer is the most energetically favorable in a solution with the highest concentration of the precipitant (0.6 M); furthermore, with its increase from 0 to 0.6 M as well as considering precipitant ions embedded in the protein (all or only cations—not essential), the energy of the dimer formation gradually decreases. It should be noted that the association of monomers with only Na ions embedded in their structure is less probable than in a case of considering both Na and Cl ions at precipitant concentrations from 0 to 0.1 M, the same at 0.2 M, and became more favorable at concentrations of 0.4 and 0.6 M, which correspond the lysozyme crystallization conditions. This fact may indicate that precipitant cations and anions embedded in the protein crystal should be studied separately as the difference of their behavior during the protein crystallization was observed.

In the absence of a precipitant in solution, the monomeric state for lysozyme molecules is more energetically favorable. Moreover, according to Table 1, even a small amount of NaCl (0.01 M) initiates dimers formation and dramatically affects the interaction between monomers while the increase in NaCl concentration from 0.05 to 0.6 M leads to a gradual energy decrease.

4. Discussion

Crystallization is well known to occur with precipitants. The addition of precipitants was considered [29] to change the interaction between monomers from repulsion to attraction. The result of our calculations substantiates these assumptions, showing that the combination of lysozyme molecules into dimers in the presence of precipitants is energetically more favorable: according to the results of the analysis of the RMSF graphs and calculations of the change in free energy during the binding of monomers, the presence of a precipitant NaCl in a lysozyme solution leads to the formation of dimers from lysozyme molecules, and the presence of precipitant ions (at least cations) embedded in the protein shows that dimers remain stable and do not decompose into monomers.

The inequality in the dimers’ formation energy values calculated for models differing in the presence of Na and Cl ions bound to the protein confirms the assumption made in [9] that precipitant ions both associated with the protein structure and located in the solution surrounding the protein affect crystallization.