Free energy decomposition of protein-protein interactions.

A free energy decomposition scheme has been developed and tested on antibody-antigen and protease-inhibitor binding for which accurate experimental structures were available for both free and bound proteins. Using the x-ray coordinates of the free and bound proteins, the absolute binding free energy was computed assuming additivity of three well-defined, physical processes: desolvation of the x-ray structures, isomerization of the x-ray conformation to a nearby local minimum in the gas-phase, and subsequent noncovalent complex formation in the gas phase. This free energy scheme, together with the Generalized Born model for computing the electrostatic solvation free energy, yielded binding free energies in remarkable agreement with experimental data. Two assumptions commonly used in theoretical treatments; viz., the rigid-binding approximation (which assumes no conformational change upon complexation) and the neglect of vdW interactions, were found to yield large errors in the binding free energy. Protein-protein vdW and electrostatic interactions between complementary surfaces over a relatively large area (1400--1700 A(2)) were found to drive antibody-antigen and protease-inhibitor binding.     

Free energy calculations for the relative binding affinity between DNA and lambda-repressor

The change in the binding free energy between DNA and lambda-repressor resulting from a base substitution, thymine (T)-->deoxyuracil (abbreviated as U), was evaluated by the free energy perturbation method on the basis of molecular dynamics simulations for the DNA-lambda-repressor complex in water with all degrees of freedom and including long-range Coulomb interactions. The binding free energy change that we calculated (1.47 +/- 0.40 kcal/mol) was in good agreement with an experimental value (1.8 kcal/mol). We clarified why the small difference between T and U (CH(3) in T is replaced with H in U) caused such a significant change in the binding free energy: The substitution of CH(3) in T with H in U lowered the dissociated-state free energy level due to the gain of the hydration free energy. Furthermore, the T-->U substitution raised the free energy level in the associated state due to the loss of the favored van der Waals (vdW) interactions with the lambda-repressor amino acid residues. In other words, the amino acid residues of lambda-repressor can recognize the CH(3) in T through the vdW interactions with the CH(3). This recognition is enhanced in a water environment, because the hydrophobic CH(3) prefers the amino acid residues of lambda-repressor to water molecules.     

Do electrostatic interactions destabilize protein–nucleic acid binding?

The negatively charged phosphates of nucleic acids are often paired with positively charged residues upon binding proteins. It was thus counter-intuitive when previous Poisson-Boltzmann (PB) calculations gave positive energies from electrostatic interactions, meaning that they destabilize protein-nucleic acid binding. Our own PB calculations on protein-protein binding have shown that the sign and the magnitude of the electrostatic component are sensitive to the specification of the dielectric boundary in PB calculations. A popular choice for the boundary between the solute low dielectric and the solvent high dielectric is the molecular surface; an alternative is the van der Waals (vdW) surface. In line with results for protein-protein binding, in this article, we found that PB calculations with the molecular surface gave positive electrostatic interaction energies for two protein-RNA complexes, but the signs are reversed when the vdW surface was used. Therefore, whether destabilizing or stabilizing effects are predicted depends on the choice of the dielectric boundary. The two calculation protocols, however, yielded similar salt effects on the binding affinity. Effects of charge mutations differentiated the two calculation protocols; PB calculations with the vdW surface had smaller deviations overall from experimental data.     

Prediction of protein-protein binding free energies

We present an energy function for predicting binding free energies of protein-protein complexes, using the three-dimensional structures of the complex and unbound proteins as input. Our function is a linear combination of nine terms and achieves a correlation coefficient of 0.63 with experimental measurements when tested on a benchmark of 144 complexes using leave-one-out cross validation. Although we systematically tested both atomic and residue-based scoring functions, the selected function is dominated by residue-based terms. Our function is stable for subsets of the benchmark stratified by experimental pH and extent of conformational change upon complex formation, with correlation coefficients ranging from 0.61 to 0.66.     

A statistical rescoring scheme for protein-ligand docking: Consideration of entropic effect

Computational prediction of protein-ligand binding modes provides useful information on the relationship between structure and activity needed for drug design. A statistical rescoring method that incorporates entropic effect is proposed to improve the accuracy of binding mode prediction. A probability function for two sampled conformations to belong to the same broad basin in the potential energy surface is introduced to estimate the contribution of the state represented by a sampled conformation to the configurational integral. The rescoring function is reduced to the colony energy introduced by Xiang et al. (Proc Natl Acad Sci USA 2002;99:7432-7437) when a particular functional form for the probability function is used. The scheme is applied to rescore protein-ligand complex conformations generated by AutoDock. It is demonstrated that this simple rescoring improves prediction accuracy substantially when tested on 163 protein-ligand complexes with known experimental structures. For example, the percentage of complexes for which predicted ligand conformations are within 1 A root-mean-square deviation from the native conformations is doubled from about 20% to more than 40%. Rescoring with 11 different scoring functions including AutoDock scoring functions were also tested using the ensemble of conformations generated by Wang et al. (J Med Chem 2003;46:2287-2303). Comparison with other methods that use clustering and estimation of conformational entropy is provided. Examination of the docked poses reveals that the rescoring corrects the predictions in which ligands are tightly fit into the binding pockets and have low energies, but have too little room for conformational freedom and thus have low entropy.     

Molecular basis of inhibitory peptide maurotoxin recognizing Kv1.2 channel explored by ZDOCK and molecular dynamic simulations

Inhibitory peptide-channel interactions have been utilized to characterize both channels and peptides; however, the fundamental basis for these interactions remains elusive. Here, combined computation methods were employed to study the specific binding of maurotoxin (MTX) peptide to Kv1.2 channel. In the first stage, numerous predicted complexes were generated by docking an ensemble of all 35 NMR conformations of MTX to Kv1.2 channel with ZDOCK program. Then the resulted complexes were clustered and classified into four main binding modes, based on experimental information and interaction energy analysis after the energy minimization and molecular dynamics (MD) simulations. By examining the stability of the plausible candidates through unrestrained MD simulations and calculation of the binding free energies, a final reasonable MTX-Kv1.2 complex was identified, with an overall high degree of correlation between the calculation and experiment on mutational effects. In the obtained complex structure model, MTX mainly used its beta-sheet domains to associate the channel mouth instead of the well-recognized functionally important S5P linkers of Kv1.2 channel. Structure analysis characterized that the most essential Tyr(32) residue of MTX was surrounded by a "pocket" formed by many nonpolar and polar residues of Kv1.2 channel, and revealed a pore-blocking Lys(23) and an important Lys(7) stabilized by strong electrostatic interactions with Asp(379) of Kv1.2. Furthermore, a stepwise structural arrangement for both ligand and receptor was found to accompany the tighter interaction of MTX into the target channel. The starting conformation of MTX, the side-chain conformation of the most important residue Tyr(32), and proper introduction of flexibility for candidate complexes were demonstrated to be considerably important factors for obtaining the final reasonable complex structure model. All these findings should not only be helpful for identifying more plausible K(+) channel-inhibitory peptide complex structures, but also provide intrinsically valuable structural biology information to interpret binding affinities, specificities, and diversity of K(+) channel-nature toxin interactions.     

Calculation of HyHel10-lysozyme binding free energy changes: Effect of ten point mutations

The change in free energy of binding of hen egg white lysozyme (HEL) to the antibody HyHel-10 arising from ten point mutations in HEL (D101K, D101G, K96M, K97D, K97G, K97G, R21E, R21K, W62Y, and W63Y) was calculated using a combination of the finite difference Poisson-Boltzmann method for the electrostatic contribution, a solvent accessible surface area term for the non-polar contribution, and rotamer counting for the sidechain entropy contribution. Comparison of experimental and calculated results indicate that because of pKa shifts in some of the mutated residues, primarily those involving Aspartate or Glutamate, proton uptake or release occurs in binding. When this effect was incorporated into the binding free energy calculations, the agreement with experiment improved significantly, and resulted in a mean error of about 1.9 kcal/mole. Thus these calculations predict that there should be a significant pH dependence to the change in binding caused by these mutations. The other major contributions to binding energy changes comes from solvation and charge charge interactions, which tend to oppose each other. Smaller contributions come from nonpolar interactions and sidechain entropy changes. The structures of the HyHel-10-HEL complexes with mutant HEL were obtained by modeling, and the effect of the modeled structure on the calculations was also examined. “Knowledge based” modeling and automatic generation of models using molecular mechanics produced comparable results. Proteins 33:39–48, 1998. © 1998 Wiley-Liss, Inc.     

Contribution of protein conformation to stereochemistry and reactivity of the active center of heme proteins and enzymes. The existence of horseradish peroxidase conformations and their possible role in the catalysis mechanism

In the spectral region 350-800 nm at 4.2 K we measured magnetic circular dichroism (MCD) spectra of the pentacoordinated complex of protcheme with 2-methylimidazole, deoxyleghemoglobin, neutral and alkaline forms of reduced horseradish peroxidase in the equilibrium states, as well as in non-equilibrium states produced by low-temperature photolysis of their carbon monoxide derivatives. Earlier the corresponding results have been obtained for myoglobin, hemoglobin and cytochromes P-450 and P-420. The energies of Fe-N (proximal His) and Fe-N(pyrroles) bonds and their changes upon ligand binding in heme proteins and enzymes were compared with those in the model heme complex thus providing conformational contribution into stereochemistry of the active site. The examples of weak and strong conformational "pressure" on stereochemistry were analysed and observed. If conformational energy contribution into stereochemistry prevails the electronic one the heme stereochemistry remains unchanged on ligand binding as it was observed for leghemoglobin and alkaline horseradish peroxidase. The change of bond energies in myoglobin and hemoglobin on ligand binding are comparable with those in protein free pentacoordinated protoheme, giving an example of weak conformational contribution to heme stereochemistry. The role of protein conformation energy in the modulation of ligand binding properties of heme in leghemoglobin relative to those in myoglobins is discussed. The most striking result were obtained in the study of reduced horseradish peroxidase in the pH region of 6.0-10.2. It was found that such different perturbations as ligand binding and heme-linked ionization of the distal amino acid residue induce identical changes in heme stereochemistry. Neither heme-linked ionization in the carbon monoxide complex nor the geometry of Fe-Co bond affect the heme local structure of photoproducts. These and other findings suggest a very low conformation mobility of horseradish peroxidase whose protein constraints appear to allow only two preferable geometries of specific amino acid residues that form the heme pocket. The role of the two tertiary structure constraints on the heme in the mechanism of horseradish peroxidase function is discussed. It is supposed that one conformation produces a heme environment suitable for two-electron oxidation of the native enzyme to compound I by hydrogen peroxide while another conformation changes the heme stereochemistry in the direction favourable for back reduction of compound I by the substrate to the resting enzyme through two one-electron steps. The switch from one tertiary structure to another is expected to be induced by substrate bind     

Probing electrostatic interactions and ligand binding in aspartyl-tRNA synthetase through site-directed mutagenesis and computer simulations

Faithful genetic code translation requires that each aminoacyl-tRNA synthetase recognise its cognate amino acid ligand specifically. Aspartyl-tRNA synthetase (AspRS) distinguishes between its negatively-charged Asp substrate and two competitors, neutral Asn and di-negative succinate, using a complex network of electrostatic interactions. Here, we used molecular dynamics simulations and site-directed mutagenesis experiments to probe these interactions further. We attempt to decrease the Asp/Asn binding free energy difference via single, double and triple mutations that reduce the net positive charge in the active site of Escherichia coli AspRS. Earlier, Glutamine 199 was changed to a negatively-charged glutamate, giving a computed reduction in Asp affinity in good agreement with experiment. Here, Lysine 198 was changed to a neutral leucine; then, Lys198 and Gln199 were mutated simultaneously. Both mutants are predicted to have reduced Asp binding and improved Asn binding, but the changes are insufficient to overcome the initial, high specificity of the native enzyme, which retains a preference for Asp. Probing the aminoacyl-adenylation reaction through pyrophosphate exchange experiments, we found no detectable activity for the mutant enzymes, indicating weaker Asp binding and/or poorer transition state stabilization. The simulations show that the mutations' effect is partly offset by proton uptake by a nearby histidine. Therefore, we performed additional simulations where the nearby Histidines 448 and 449 were mutated to neutral or negative residues: (Lys198Leu, His448Gln, His449Gln), and (Lys198Leu, His448Glu, His449Gln). This led to unexpected conformational changes and loss of active site preorganization, suggesting that the AspRS active site has a limited structural tolerance for electrostatic modifications. The data give insights into the complex electrostatic network in the AspRS active site and illustrate the difficulty in engineering charged-to-neutral changes of the preferred ligand.     

Focusing of the electrostatic potential at EF-hands of calbindin D(9k): titration of acidic residues

Biological functions for a large class of calmodulin-related proteins, such as target protein activation and Ca(2+) buffering, are based on fine-tuned binding and release of Ca(2+) ions by pairs of coupled EF-hand metal binding sites. These are abundantly filled with acidic residues of so far unknown ionization characteristics, but assumed to be essential for protein function in their ionized forms. Here we describe the measurement and modeling of pK(a) values for all aspartic and glutamic acid residues in apo calbindin D(9k), a representative of calmodulin-related proteins. We point out that while all the acidic residues are ionized predominantly at neutral pH, the onset of proton uptake by Ca(2+) ligands with high pK(a) under these conditions may have functional implications. We also show that the negative electrostatic potential is focused at the bidental Ca(2+) ligand of each site, and that the potential is significantly more negative at the N-terminal binding site.     

Protein-protein binding is often associated with changes in protonation state

pK(a) values of ionizable residues have been calculated using the PROPKA method and structures of 75 protein-protein complexes and their corresponding free forms. These pK(a) values were used to compute changes in protonation state of individual residues, net changes in protonation state of the complex relative to the uncomplexed proteins, and the correction to a binding energy calculated assuming standard protonation states at pH 7. For each complex, two different structures for the uncomplexed form of the proteins were used: the X-ray structures determined for the proteins in the absence of the other protein and the individual protein structures taken from the structure of the complex (referred to as unbound and bound structures, respectively). In 28 and 77% of the cases considered here, protein-protein binding is accompanied by a complete (>95%) or significant (>50%) change in protonation state of at least one residue using unbound structures. Furthermore, in 36 and 61% of the cases, protein-protein binding is accompanied by a complete or significant net change in protonation state of the complex relative to the separated monomers. Using bound structures, the corresponding values are 12, 51, 20, and 48%. Comparison to experimental data suggest that using unbound and bound structures lead to over- and underestimation of binding-induced protonation state changes, respectively. Thus, we conclude that protein-protein binding is often associated with changes in protonation state of amino acid residues and with changes in the net protonation state of the proteins. The pH-dependent correction to the binding energy contributes at least one order of magnitude to the binding constant in 45 and 23%, using unbound and bound structures, respectively.     

Integrin alphaIIbbeta3:ligand interactions are linked to binding-site remodeling

This study tested the hypothesis that high-affinity binding of macromolecular ligands to the alphaIIbbeta3 integrin is tightly coupled to binding-site remodeling, an induced-fit process that shifts a conformational equilibrium from a resting toward an open receptor. Interactions between alphaIIbbeta3 and two model ligands-echistatin, a 6-kDa recombinant protein with an RGD integrin-targeting sequence, and fibrinogen's gamma-module, a 30-kDa recombinant protein with a KQAGDV integrin binding site-were measured by sedimentation velocity, fluorescence anisotropy, and a solid-phase binding assay, and modeled by molecular graphics. Studying echistatin variants (R24A, R24K, D26A, D26E, D27W, D27F), we found that electrostatic contacts with charged residues at the alphaIIb/beta3 interface, rather than nonpolar contacts, perturb the conformation of the resting integrin. Aspartate 26, which interacts with the nearby MIDAS cation, was essential for binding, as D26A and D26E were inactive. In contrast, R24K was fully and R24A partly active, indicating that the positively charged arginine 24 contributes to, but is not required for, integrin recognition. Moreover, we demonstrated that priming--i.e., ectodomain conformational changes and oligomerization induced by incubation at 35 degrees C with the ligand-mimetic peptide cHarGD--promotes complex formation with fibrinogen's gamma-module. We also observed that the gamma-module's flexible carboxy terminus was not required for alphaIIbbeta3 integrin binding. Our studies differentiate priming ligands, which bind to the resting receptor and perturb its conformation, from regulated ligands, where binding-site remodeling must first occur. Echistatin's binding energy is sufficient to rearrange the subunit interface, but regulated ligands like fibrinogen must rely on priming to overcome conformational barriers.     

Probing metal ion binding and conformational properties of the colicin E9 endonuclease by electrospray ionization time-of-flight mass spectrometry

Nano-electrospray ionization time-of-flight mass spectrometry (ESI-MS) was used to study the conformational consequences of metal ion binding to the colicin E9 endonuclease (E9 DNase) by taking advantage of the unique capability of ESI-MS to allow simultaneous assessment of conformational heterogeneity and metal ion binding. Alterations of charge state distributions on metal ion binding/release were correlated with spectral changes observed in far- and near-UV circular dichroism (CD) and intrinsic tryptophan fluorescence. In addition, hydrogen/deuterium (H/D) exchange experiments were used to probe structural integrity. The present study shows that ESI-MS is sensitive to changes of the thermodynamic stability of E9 DNase as a result of metal ion binding/release in a manner consistent with that deduced from proteolysis and calorimetric experiments. Interestingly, acid-induced release of the metal ion from the E9 DNase causes dramatic conformational instability associated with a loss of fixed tertiary structure, but secondary structure is retained. Furthermore, ESI-MS enabled the direct observation of the noncovalent protein complex of E9 DNase bound to its cognate immunity protein Im9 in the presence and absence of Zn(2+). Gas-phase dissociation experiments of the deuterium-labeled binary and ternary complexes revealed that metal ion binding, not Im9, results in a dramatic exchange protection of E9 DNase in the complex. In addition, our metal ion binding studies and gas-phase dissociation experiments of the ternary E9 DNase-Zn(2+)-Im9 complex have provided further evidence that electrostatic interactions govern the gas phase ion stability.     

The role of flexibility and hydration on the sequence-specific DNA recognition by the Tn916 integrase protein: a molecular dynamics analysis

The N-terminal domain of the Tn916 integrase protein (INT-DBD) is responsible for DNA binding in the process of strand cleavage and joining reactions required for transposition of the Tn916 conjugative transposon. Site-specific association is facilitated by numerous protein-DNA contacts from the face of a three-stranded beta-sheet inserted into the major groove. The protein undergoes a subtle conformational transition and is slightly unfolded in the protein-DNA complex. The conformation of many charged residues is poorly defined by NMR data but mutational studies have indicated that removal of polar side chains decreases binding affinity, while non-polar contacts are malleable. Based on analysis of the binding enthalpy and binding heat capacity, we have reasoned that dehydration of the protein-DNA interface is incomplete. This study presents results from a molecular dynamics investigation of the INT-DBD-DNA complex aimed at a more detailed understanding of the role of conformational dynamics and hydration in site-specific binding. Comparison of simulations (total of 13 ns) of the free protein and of the bound protein conformation (in isolation or DNA-bound) reveals intrinsic flexibility in certain parts of the molecule. Conformational adaptation linked to partial unfolding appears to be induced by protein-DNA contacts. The protein-DNA hydrogen-bonding network is highly dynamic. The simulation identifies protein-DNA interactions that are poorly resolved or only surmised from the NMR ensemble. Single water molecules and water clusters dynamically optimize the complementarity of polar interactions at the 'wet' protein-DNA interface. The simulation results are useful to establish a qualitative link between experimental data on individual residue's contribution to binding affinity and thermodynamic properties of INT-DBD alone and in complex with DNA.     

An insight into the inhibitory selectivity of 4-(Pyrazol- 4-yl)-pyrimidines to CDK4 over CDK2

Revealing selectivity mechanism of cyclin-dependent kinases (CDK) and their inhibitors is an important issue to develop potential anticancer drugs. The substituted 4-(Pyrazol-4-yl)-pyrimidines are potent inhibitors of CDK4 but not of the highly homologous CDK2. In order to reveal the inhibitory selectivity of these inhibitors to CDK4 over CDK2, we select one of substituted 4-(Pyrazol-4-yl)-pyrimidines as a representative (marked as A1 hereunder) and perform molecular docking, molecular dynamics simulations and binding free energy analysis for CDK4/A1 and CDK2/A1, respectively. The electrostatic and van der Waals (vdW) interactions of the A1 inhibitor with CDK4/CDK2 are discussed. The computed binding free energies based on the MM-PBSA method are consistent with experimental bioactivity ranking of A1 inhibitor to CDK4/CDK2. On the other hand, the conformational characteristics of CDK2 and CDK4 induced by A1 inhibitor are analysed and revealed. Results demonstrate that the vdW interactions considerably contribute to binding of CDK4/CDK2 with A1 inhibitor and are similar in size. The hydrogen bonding between A1 inhibitor and CDK4/CDK2 is considerably favourable to the binding, in which the hydrogen bond between the NH group of the pyrazole group of A1 and the residue Asp158 of CDK4 plays a crucial role in inhibitory selectivity of A1 inhibitor to CDK4 over CDK2. The electrostatic interaction energy differences between the corresponding residues of CDK4/A1 and CDK2/A1 confirm the above inference. The conformational changes of CDK2 and CDK4 induced by A1 inhibitor influence the selectivity of A1 inhibitor to CDK4/CDK2.     

Empirical calculation of the relative free energies of peptide binding to the molecular chaperone DnaK

We describe a methodology to calculate the relative free energies of protein-peptide complex formation. The interaction energy was decomposed into nonpolar, electrostatic and entropic contributions. A free energy-surface area relationship served to calculate the nonpolar free energy term. The electrostatic free energy was calculated with the finite difference Poisson-Boltzmann method and the entropic contribution was estimated from the loss in the conformational entropy of the peptide side chains. We applied this methodology to a series of DnaK*peptide complexes. On the basis of the single known crystal structure of the peptide-binding domain of DnaK with a bound heptapeptide, we modeled ten other DnaK*heptapeptide complexes with experimentally measured K(d) values from 0.06 microM to 11 microM, using molecular dynamics to refine the structures of the complexes. Molecular dynamic trajectories, after equilibration, were used for calculating the energies with greater accuracy. The calculated relative binding free energies were compared with the experimentally determined free energies. Linear scaling of the calculated terms was applied to fit them to the experimental values. The calculated binding free energies were between -7.1 kcal/mol and - 9.4 kcal/mol with a correlation coefficient of 0.86. The calculated nonpolar contributions are mainly due to the central hydrophobic binding pocket of DnaK for three amino acid residues. Negative electrostatic fields generated by the protein increase the binding affinity for basic residues flanking the hydrophobic core of the peptide ligand. Analysis of the individual energy contributions indicated that the nonpolar contributions are predominant compared to the other energy terms even for peptides with low affinity and that inclusion of the change in conformational entropy of the peptide side chains does not improve the discriminative power of the calculation. The method seems to be useful for predicting relative binding energies of peptide ligands of DnaK and might be applicable to other protein-peptide systems, particularly if only the structure of one protein-ligand complex is available.     

Insights into activation and RNA binding of trp RNA-binding attenuation protein (TRAP) through all-atom simulations

Tryptophan biosynthesis in Bacillus stearothermophilus is regulated by a trp RNA binding attenuation protein (TRAP). It is a ring-shaped 11-mer of identical 74 residue subunits. Tryptophan binding pockets are located between adjacent subunits, and tryptophan binding activates TRAP to bind RNA. Here, we report results from all-atom molecular dynamics simulations of the system, complementing existing extensive experimental studies. We focus on two questions. First, we look at the activation mechanism, of which relatively little is known experimentally. We find that the absence of tryptophan allows larger motions close to the tryptophan binding site, and we see indication of a conformational change in the BC loop. However, complete deactivation seems to occur on much longer time scales than the 40 ns studied here. Second, we study the TRAP-RNA interactions. We look at the relative flexibilities of the different bases in the complex and analyze the hydrogen bonds between the protein and RNA. We also study the role of Lys37, Lys56, and Arg58, which have been experimentally identified as essential for RNA binding. Hydrophobic stacking of Lys37 with the nearby RNA base is confirmed, but we do not see direct hydrogen bonding between RNA and the other two residues, in contrast to the crystal structure. Rather, these residues seem to stabilize the RNA-binding surface, and their positive charge may also play a role in RNA binding. Simulations also indicate that TRAP is able to attract RNA nonspecifically, and the interactions are quantified in more detail using binding energy calculations. The formation of the final binding complex is a very slow process: within the simulation time scale of 40 ns, only two guanine bases become bound (and no others), indicating that the binding initiates at these positions. In general, our results are in good agreement with experimental studies, and provide atomic-scale insights into the processes.     

Binding of 1-substituted carbazolyl-3,4-dihydro-β-carbolines with DNA: Molecular dynamics simulation and MM-GBSA analysis

Molecular Mechanics-Generalized Born-Solvent Accessibility free energy calculations were used to analyse DNA binding affinity of 1-substituted carbazolyl-3,4-dihydro-β-carboline molecules. In this study, DNA structure with sequence of d(CGATCG)2 was used for simulations. 15 ns molecular dynamics simulations of the studied complexes were performed. The calculated free energy was compared with experimental antitumor activity (IC₅₀). The predicted free energies decreased with the increase of IC₅₀ values. It was shown that molecules 1–6 bind to DNA via intercalation mode, while molecules 7–9 bind through groove binding mode. Also, it was found that the vdW energy term (ΔE_vdW) and the non-polar desolvation energy (ΔG_SA) are the favorable terms for binding energy, whereas net electrostatic energies (ΔE_ele + ΔG_GB) and conformational entropy energy (TΔS) are unfavorable ones.