Team:Thessaloniki/Measurement
Overview
Our team has invested a lot of effort into developing models that simulate all parts of our project. The need to verify the experimental results regarding the calculation of the binding affinity of the DNA-Protein reaction was covered implementing Molecular Dynamics techniques in our model.
To assess the experimental results, we focused on calculating in silico the binding free energy of the reaction. Molecular Dynamics simulations were conducted on the DNA-Protein complex using the MM/GBSA and the MM/PBSA methods. Comparing the two and observing the controversy in the two different methods, we evaluated the experimental data and provided valuable feedback regarding the credibility of the results.
Molecular Dynamics were conducted in order to simulate the binding reaction between the complex of NF-KB P65 HOMODIMER and the DNA binding site. The input data for the simulation were derived from the 1RAM structure provided by the Protein Data Bank (PDB). From a variety of MD software we concluded in using two of the most accurate suits. Οverall, we used the AMBER suit for both MM/GBSA and MM/PBSA calculations and the Schrodinger-Desmond software, for on MM/GBSA calculation.
Binding Afinity
Binding affinity is the strength of the binding interaction between a single biomolecule to the ligand/binding partner. This strength is absolutely proportional to the binding free energy ΔG. ΔG of binding is the calculated difference between the energy of the reactants and the products.The lower the binding energy, the stronger the bond between receptor and ligand.
Figure 1: The gibbs energy of a system during a reaction. The reactants create an intermediate complex with high energy and then, while the final products are forming, the energy drops.
In order to experimentally calculate ΔG we have to measure the concentrations of protein, DNA and Protein-DNA complex at the equilibrium state. Ιn more detail, using the concentration data we can calculate the equilibrium constant and therefore the Binding Free Energy of the reaction.
Molecular Dynamics Theory
The process of calculating the Gibbs Energy through a Molecular Dynamic simulation is quite different from the theoretical approach.
As we describe in more detail below, we first produce a trajectory of the reaction occuring for 100ns. Then we calculate the Binding Free Energy of each frame of the trajectory for the complex, the receptor and the ligand, starting from the frame that the system reached equilibrium. The total Gibbs Binding Free Energy is calculated according to the following equation:
We minimized the idle time of our experiments by using a single trajectory approach, which means that all ensembles are extracted from a single molecular dynamic trajectory. That being said, we are making the assumption that the receptor and ligand ensembles are comparable in the bound and the unbound state.
At this point of our experiments we came up on the obstacle of not knowing which method to use to calculate the Binding Free Energy. Molecular Dynamics simulations were not used from professors at our home University. Therefore, we contacted Professor Kolocouris A. at National and Kapodistrian University of Athens. After a skype session with Mr Kolocouris, we were introduced to one of his Phd students, Dimitris Stamatis. Dimitris provided assistance on many levels of our work, among others he helped us chose the best method to calculate the Binding Free Energy.
Among the multiple binding free energy calculation methods, we chose the most used end-state methods, MM/GBSA and MM/PBSA, as they are less time consuming and less theoretically rigorous than any other corresponding pathway-based methods.
According to MM-GB/PB-SA methods, in order to calculate the binding free energy, it is being broken down to its constituents. Every aspect of the equation (2) is calculated separately according to the following equations:
For each term on the right side of (Eq. 2):
(Eq. 3)
(Eq. 4)
(Eq. 5)
The table below shows the method used for each component:
Poisson-Boltzmann & Generalized Born
The Poisson-Boltzmann equation (PB) describes the theoretical electrostatic interactions of a solute in a solvent which contains ions. Usually this equation is used with approximations as it is computationally expensive.
An approximation to the linearized Poisson-Boltzmann equation is the Generalized Born model. Here, the solute is considered to be composed of spheres whose internal dielectric constant differs from the external solvent.
Input Preparation
We derived our input files (1RAM) from Protein Data Bank. The input structure of the DNA - NF-KB P65 HOMODIMER, did not contain the hydrogens nor the ions for the neutralization. Therefore, prior to conducting the simulation we had to add the hydrogens. Then we applied the appropriate force fields. A force field describes the forces that are applied on each atom. Finally ions and water molecules were added to the solvent.
It is important to note that at the AMBER suit we used the leaprc.protein.ff14SB and the leaprc.DNA.OL1 force fields, while in Desmond, we used the OPLS_2005 force field.
Figure 2: Before input preparation the the structure consists of the protein and the DNA without hydrogens ions or water (left). After adding all the components, the system is a box that contains the complex and the ions as solutes and water as solvent (right).
Minimization
The structure obtained after the preparation contains a geometry which may not correspond to the actual minima for the forces we are using, because the addition of components. This will result in conflicts and overlaps between the atoms of the system. So prior to any other simulations we are minimizing the structure. Minimization results in rearranging the positions of the atoms in order to satisfy the implemented force fields. The minimization time was fairly short in comparison with the total simulation.
Temperature and Pressure
Τill this point our whole system has a temperature of 0 K. Our wet lab experiment was conducted at 310 K, so we have to heat the system up at the desired temperature, with a method that will equally distribute the heat. It is required that our system will be heated from the initial temperature on a slow rate to avoid disrupting the solvation box limits.
Figure 3: The temperature of the system in Kelvin for the first 120 ps of simulation. The graph was obtained from the simulation with AMBER.
Initially in our solvated system the pressure is zero. The experiment is being conducted under constant, atmospheric, pressure. Heating the system under constant pressure results in changes in the volume of the solvated box. To avoid this, during the heating we keep the volume constant and the pressure can change. Once the desired temperature is achieved, the restraints are reversed. The desired outcome is to preserve a constant average pressure, and not a constant pressure in each frame.
In the figure (4) below we observe that for the first 20ps of the simulation, the pressure remains close to zero. This is not irregular as the simulation occured for constant volume. After 20ps, although we observe the pressure rate fluctuating wildly. As the simulation time passes by, the pressure stabilizes at 1 atm rate. We observe the stabilization of the mean value at an estimated time of 60ps.
Figure 4: The pressure value is calculated in atm (left) and the volume value in A3 (right) of the system between 20-120 ps of simulation. The graphs are obtained from simulation analysis with AMBER suite.
Equilibrium
After successfully preparing our system, we proceed with the generation of the trajectory files. To complete an accurate Molecular Dynamics simulation it is crucial to apply a sufficient time period for the simulation to occur till the equilibrium state. Reaching the equilibrium state is equivalent to reaching the minimum energy value. We conducted multiple simulations in order to evaluate the time our system needs to equilibrate. We came to the conclusion that 100 ns were more than enough for our simulation. It is important to state that the data needed to calculate the free energy of binding are derived the time period that our system is at the equilibrium state.
We are aware that our system has reached equilibrium when the average distance between the heavy atoms of the system (C,N,O) is stabilizing or fluctuates around an average value. This evaluation measure Root-mean-square deviation of atomic positions (RMSD) was used to reach our conclusions. The RMSD measure is reliable, according to bibliography, in order evaluate that a system has reached equilibrium.
Figure 5 : RMSD of all atoms of the complex (up-left), RMSD of DNA atoms (up-right) and the RMSD of the protein backbone (down middle). The protein Backbone contains the A carbons of the protein.The graphs are obtained from the simulation conducted with Schrodinger suite.
MM-GB/PB-SA
As described in theory above, we used MM/GBSA and MM/PBSA methods to analyze the frames of the simulation that has reached equilibrium.Using the mean values of the results we conduct the required calculations as shown in Eq.2 above. It is important to be mentioned that these methods presents further capacities ruther than calculating the binding free energy. Within the needs of our project it was not required to use any other capacity of the methods.
Software
We used two software packages for the Molecular Dynamics simulations and the analysis of the results. On our first approach to Molecular Dynamics we used Schrodinger suite. Schrodinger provides a user friendly environment and can be accessed on a free non commercial license.
The second software we used is AMBER, a suite of biomolecular simulation programs. At the beginning of our experiments we used AmberTools19, the free module that AMBER community provides.
In detail, we used Schrodinger free version, distributed from DE Shaw to prepare the system and conduct the Molecular Dynamics simulations using a GPU. After the simulation procedure was completed we converted all files derived from Schrodinger to be compatible with AmberTools19 file format. We quickly abandoned the above methodology as it was too difficult, time consuming and most importantly unable to provide us with trustworthy results.
We reached out to AMBER community and Schrodinger for assistance. We were delighted to hear that both companies could provide us with a free full licence of the programs for the competition.
Meanwhile, we came in conduct with Professor Kolocouris Antonios, Associate Professor at the department of Pharmaceutical Chemistry and Organic Chemistry at the National and Kapodistrian University of Athens. Dr Kolocouris referred us to his M.Sc. student, Stamatis Dimitris, member of the Pharmaceutical chemistry laboratory at the National and Kapodistrian University of Athens. Dimitris from that moment and after was our dry lab mentor and provided assistance at every problem we occured.
Having received the full licenses for both Schrodinger and AMBER we redesigned our simulations.
We first approached Schrodinger software as it provides a user friendly graphic environment and minimizes the time needed to prepare the simulations. Nonetheless, the results obtained from the RMSD of the ligand (Figure 5) depict high abnormal values. Additionally, the calculated Binding Free Energy from the MMGBSA method is also surprisingly high as the results fluctuate around the value of 200 kcal/mol. We came to the conclusion that Schrodinger software could not support efficiently our simulations. The results derived were not compatible with the expected values. After conducting the support team of Schrodinger we were informed that primal modules of Schrodinger lack support for DNA-Protein binding.
Consequently, we conducted our simulations using AMBER software. The calculations derived from AMBER are presented in detail on the result section.
The code used for the calculations can be found in our GitHub account.
Schrodinger - Desmond
As mentioned above, calculations for Binding Free Energy derived from Schrodinger was abnormally high and we discarded the results. More specifically, from the frames at equilibrium state of a 100 ns simulation, the calculated free binding energy with a MMGBSA method was -283.27 kcal/mol. This results to abnormal values of equilibrium constants, K_bind = 3.734*10^(-210) and K_detachment = 2.678 * 10^(209). The K_detachment value according to bibliography is fluctuated between 464 to 341 nM. The error appearing in the difference between the values indicates that either the software or the method (MMGBSA) is appropriate for these calculations.
Using AMBER we also performed a simulation of 100 ns. From this simulation we used some of the frames that were at equilibrium to perform the MMGBSA and the MMPBSA methods.