DNA interaction with a benzimidazole-containing ligands

^{The work done on the Abalone program.}

There are a large number of substances that interact with DNA. In general, these are proteins carry out a variety of biological functions. Most proteins interact with DNA in a wide groove of the double helix. That is not surprising, since a wide groove is more accessible and provides a greater variety of contacts, which can be used to build as a general interaction as well as recognition of unique sites.

In contrast low molecular weight substances such as toxins, antibiotics, dyes are often bound with a narrow groove of DNA. The narrow groove has a width which allows it to accept the small molecules with the thickness of a one - two atoms. In addition to the small thickness, a substance that binds with a narrow groove, should have an elongated banana-like shape so that it could if do not form a spiral, similar to the DNA double helix, but at least to form a short bend of the right size. In order to enter like a key into the lock the molecule must be either a rigid or semi-rigid. For example such properties are oligomers of benzimidazole, based on which was synthesized the number of DNA-binding ligands. The figure shows the dye Hoechst 33258, which contains two benzimidazole residues.

In this example we will demonstrate the main stages of the development of molecular-mechanical model, but will not to do this in detail. There are developed residues in the chains constructor, and we'll use them as needed. Nevertheless, we will describe how these models have been constructed.

First of all, we need benzimidazole residues. Draw a graph of the benzimidazole dimer and add the methyl residues to the ends. They will help us to create versatile building blocks for polymers. We contract the hand-drawn graph in three-dimensional structure by the optimizer. Then, a more accurate optimization can be done with the help of quantum chemistry. For example, optimize the structure by semi-empirical method RM1; then ab initio DFT PBEO in the MINI basis set and, finally, in the 6-31+G(d) basis set. In principle, it would be nice to have a final calculation in a larger basis, such as 6-311+G(2d,p), taking into account electron correlation at MP2, but our model is too high already. Therefore, we restrict ourselves to 6-31+G(d) with the DFT. This is enough to get fairly accurate values of the valence parameters. Multi-stage optimization is twofold. First, the calculations in large basis sets are very durable, so before we go to them we should do everything that can be done by the less expensive methods. Second, if we draw a model with the very far from reality geometry, quantum-chemical calculation may not converge. The larger (and more expensive) basis set we use the more accurate initial model we must take in order to reduce the risk of procedure divergence. When we moved from MINI basis to the 6-31+G(d) we just had a problem with the divergence, and to overcome it we had to make several attempts of optimization from the multiple starting points. For the same two reasons it makes sense not only gradually expand the basis set, but to do optimization in several steps even in each of these basis sets. Preliminary optimization is best done using coarse convergence criteria, and only after a dozen or so iterations it makes sense moving to more accurate criteria.

So, we have a model with a fairly accurate geometry. But this model is still not suitable for molecular-mechanical calculations. We need to develop parameters of valence interactions, the parameters of non-valence interactions, such as partial charges on atoms and Lennard-Jones potentials, as well as to set up conditions in which we will simulate. For example, we have to decide whether we'll put it in explicit water or use an implicit solvent model.

Our system has a very rigid conformation. The only confirmation parameter, which can vary widely, is the torsion angle between the benzimidazole residues. In order to get the torsion parameters we had to carry out quantum chemical calculations so large model. The exact values of other parameters of the valence interactions are not very important in such rigid system and we can take them the same as available parameters for aromatic systems. Focusing on molecular-mechanical types of atoms in the AMBER99SB force field, we adopt the following classification of atoms:

With this classification, most of the valence parameters are already defined in the AMBER99SB force field. The rest we are adding similar, except for the torsion angles between residues. These parameters we choose so that the difference in energy that occurs in torsional rotation was close to obtained in quantum-chemical calculations.

The figure shows the torsional barrier, calculated by DFT in the 6-31+G(d) and in the MINI basis sets. The third graph shows effect of our parameters added to the AMBER99SB force field. We have tried to approximate the results obtained in a larger basis, but does not reproduce the exact value of the barrier. Note the different behavior of the graphs near the planar conformation (zero angle). The graph looks like a sine wave in the MINI basis. While in the 6-31+G(d) basis it shows the distortion that is connected with the overlap of the hydrogen atoms in a planar conformation. It is known that the 6-31+G(d) basis is sufficient to reproduce the torsional energy, but too rough to calculate the energy of Van der Waals forces. Therefore, we tried to reproduce the energy in the 90 degrees region. In the of 0 degrees region we have to rely on the molecular mechanical parameters of the Lennard-Jones potential.

Since we have used in the existing AMBER atom types, Lennard-Jones parameters are already defined. We did not modify them.

It remains to determine the partial charges on atoms. Unfortunately, this is the most uncertain part of the procedure. In the AMBER charges are determined by fitting the electrostatic field generated by the partial charges to the field calculated by quantum-chemical method. This is a very sensible approach, since we need partial charges to approximate the electrostatic field. But there are some problems:

• There is a large arbitrariness in the choice of points for the electrostatic potential calculation. The results strongly depend on this choice.
• This fit works well only for the boundary atoms. The charges on the atoms are in the depths of the molecule, away from the points at which we calculate the electrostatic potential, and especially the atoms in a symmetrical environment, badly calculated.
• The charges calculated for the molecule in different conformations can be very different.
• If we determine not only the charges for a particular molecule, but for the remainder which will be used for construction of large molecules, there is a problem of portability, as well as the problem of joining with other residues. For this reason, the remainders are often made out of charge, even for this we have to significantly deviate from the quantum chemical results.

Nevertheless, potentially derived charges remain the best choice.

Having the residues for construction of the DNA and DNA-binding ligands, we can leave quantum chemistry and to do molecular mechanical simulations. We also need the model of the solvent. If we are going to model in explicitly represented water, we need to build a model with counterions and put it in the water cube. Such models are the most accurate, but they require large computational resources. If we omit water, it is necessary to use some compensating potential. DNA models clearly illustrate this situation since the DNA properties are largely determined by electrostatic forces. If we simply remove the water (leaving the counterions) and run molecular dynamic calculation the molecule rapidly gather into a lump (Fig. a). If we remove the counterions then the like charges untwist DNA helix and stretch the threads in opposite directions (Fig. b).

The simplest solution is to remove counterions and use distance-dependent dielectric constant, that is to accept ε = r (Fig. c). Another way is to reduce the charges on the phosphate groups from -1 to -0.3 - 0.5 electrons (Fig. d). In both cases we obtain the model, about preserving the structure of DNA. However, they are very rough. The better models are obtained if not simply assume the dielectric constant equal to the distance but to use a specially selected function. These functions have sigmoidal shape. Several their types are described in the literature. Another method is the Generalized Born model (Fig. e). However, calculations in the Generalized Born model costing 3-5 times more expensive than in the vacuum. If the effective Born radii of atoms are assumed to be constant the calculations can be significantly accelerated. This is the Sheffield model¹⁾ (рис. f). Calculations are slower 2-fold only, while a B-DNA form remains reasonable (for proteins “Sheffield model” gives bad results). That's what we use.

So, we have the ligand and DNA models. We chose an implicit water model. Now we can assemble complex. Because the shape of the DNA narrow groove is similar to our ligands shape it’s easy to do. We can either put the ligand in the groove with the mouse (Fig.) or place it at a certain distance and then drawn into the groove with the aid of constrains. In any case on the first phase these manipulations should be done with the “frozen” DNA. Otherwise, we can greatly spoil its shape.

At this stage (with frozen atoms) should not be carefully optimized. It's sufficient to eliminate the strong overlap of atoms.

Then we can make more thorough optimization unfreezing DNA.

We have optimized model. It can be used as a start for the molecular dynamics. But before we go into dynamics under normal conditions, it can be further optimized with the help of dynamics at low temperature. For example, in the range 50-100 K. The calculation should be long enough for at least a hundred picoseconds (nanosecond better). In our case, 1-ns dynamics at 50 K did not lead to significant changes.

We perform molecular dynamic simulation of the resulting model. Exchange of the 6th replicas at temperatures of -25 - +25 ° C, 2.5 ns for each replica (total time of 15 ns). We see that the complex is stable, the DNA retains the B-form, the ligand retains its position in the narrow groove. The ligand position is fixed by two hydrogen bonds formed by NH groups of benzimidazole rings and the oxygen and nitrogen atoms of DNA base pairs.

In the case of AT pairs in the narrow groove of DNA exhibited two potential hydrogen bond acceptor, in the case of the GC pair the donor and acceptor hydrogen bond. Guanine hydrogen prevents ligands entering deeply into the narrow groove. Therefore, most of the ligands are donors of hydrogen bonds and binds to AT-rich DNA sites.

Numerous attempts to create similar ligands which recognize GC-containing sites are mainly associated with the replacement of a ligand NH hydrogen bond donor to the oxygen, nitrogen or sulfur, which can form hydrogen bonds in the opposite direction, accepting the NH hydrogen group of guanine.

So, we have a model of the ligand in the narrow groove of DNA. Let us make the next step. It is known that narrow groove ligands can bind to DNA, not only individually, but also in pairs. There is evidence of this type of binding and for benzimidazole ligands. Make a model. First of all prepare a “sandwich” of two ligand molecules.

To convenient manipulations of this pair make a fictitious bond.

We extend a narrow groove in our previous model. This can be done with the help of restrains. Freeze the dye, one strand of DNA and nucleotides at the ends of the second strand. We introduce repulsive restrains for the second strand and ligand, and spend a short dynamics at low temperature (50 K, followed by 10 K). The figure shows the resulting model with an extended groove. Restrains were applied to selected atoms.

Add a sandwich of two ligands overlapping one of the dies with the dye in the model. Erase the auxiliary bond and excess dye. Eliminate the strong atomic overlap by the short optimization.

Remove the freezing, and optimize the model using the conjugate gradients and then dynamics at 10-50K. Then half a nanosecond replica exchange molecular dynamics (250-300K).

The complex is stable. The shape of DNA is not much distorted, despite the two dye molecules in the narrow groove. It is interesting that the relative position of the dye molecules is not conserved. They easily slide against each other in the depth of the grooves. At the same time, their position along the DNA strands is quite stable; the network of hydrogen bonds is preserved.

Now try to make a GC-binding ligand by replacing donor of the hydrogen bond to the acceptor, namely the benzimidazole residue to benzoxazole.

Such a substance is really capable to bind with the GC-site of DNA. But what happens on the AT-site? Unfortunately, the lack of a hydrogen bond does not lead to sufficient destabilization of the complex and the ligand binds both to a GC- and to AT-sites. The figure shows the complex of benzoxazole-containing ligand with the same dodeca-nucleotide CGCGAATTCGCG. Distortions in the geometry of the complex are not observed, the interaction energy is large.

At present there is progress in the creation of GC-ligand, but a really good solution has not been found.