The Octahydroindene Carboxyl Substructure from Dihydrobetulinic Acid is Essential to Inhibit Topoisomerase IB from Leishmania donovani

The dihydrobetulinic acid is a known competitive inhibitor of topoisomerase IB from Leishmania donovani, a validated target for developing new antileishmanial drugs. However, its binding mode and interaction pocket have not been established yet. We combined docking and molecular dynamics simulations to identify the most probable binding pocket. Our best model strongly suggests a cavity involving the residues arginine 314 and arginine 410 from chain A, and the catalytic tyrosine 222 from chain B as the interaction site and a substructure of this terpene inhibitor as essential for the process of molecular recognition. Then, a new class of inhibitors with increased affinity could be designed by structure-based approaches.


Introduction
Visceral leishmaniasis disease (VL) or Kala-azar is the second largest parasitic killer in the world (only malaria is more fatal) and infects an estimated 200-400 thousand people each year. 1,2VL is caused by the Leishmania donovani (L.donovani) protozoan and transmitted by phlebotomine sand flies.After biting, the parasite migrates to internal organs such as liver, spleen and bone marrow causing hepatosplenomegaly, severe anemia and damages to the host's immune system. 3This disease is almost always fatal if untreated, and its pharmacological therapy is limited. 3entavalent antimonials and meglumine antimoniate are the first-line treatment in most parts of the world despite the low efficacy and adverse reactions.Amphotericin B is the second choice used in different formulations, such as the liposome, and pentamidine is another antileishmanial drug used in endemic areas where antimonials are inefficient.Oralmiltefosine and paromomycin have been developed and introduced in the place of antimonials. 4Although miltefosine shows moderate toxicity, its terotogenicity imposes a problem due to the time required for the complete treatment, while the toxicity of paromomycin and the rise of resistance 4 contribute to enhance the limitations for the antileishmanial chemotherapy.Then, limited therapeutic options, emergence of resistant strains, toxicity and teratogenicity associated with current antileishmanial drugs demand continuous development of new chemotherapeutic agents 5 and the identification of more suitable molecular targets as well. 6NA (deoxyribonucleic acid) topoisomerases have emerged as striking candidate targets for the therapy of leishmaniasis. 7These ubiquitous enzymes are involved in many vital cellular process catalyzing changes in the topological state of duplex DNA during replication, transcription, recombination and DNA repair, by introducing transient single strain breaks in the nucleic acid backbone. 8The DNA topoisomerases are classified into type I and II, both of equal importance as therapeutic targets against leishmaniasis. 9,10Leishmania topoisomerase IB possess heterodimeric structure consisting of a large subunit of 635 amino acids (LdTOP1L) and a small subunit of 262 amino acids (LdTOP1S). 11opoisomerase inhibitors for antileishmanial therapy are of immense interest.These inhibitors are classified into two classes: the class I inhibitors stabilize the formation of topoisomerase-DNA covalent complex (cleavable complex) and they are also known as "topoisomerase poisons", while inhibitors from class II avoid the formation of the covalent complex and inhibit the catalytic activity of the enzyme. 8Camptothecin (CPT), the most studied topoisomerase I inhibitor, is an antitumor agent and promotes the formation of the protein-DNA cleavable complex leading to apoptosis-like cell death in L. donovani. 12On the other hand, dihydrobetulinic acid (DHBA), a derivative from the natural product betulinic acid, has been reported as a highly potent antileishmanial agent and a promising compound to generate a lead.DHBA is dual inhibitor that prevents the formation of topoisomerase-DNA binary complex.Besides that, it is a competitive inhibitor that also induces the apoptosis. 13For this reason, the investigation of its mode of interaction, the identification of the binding site, and the residues that are involved in this molecular recognition are key information to aid the synthesis of new inhibitor derivatives with improved potency and selectivity for topoisomerase type I (topIB).Establishing a model to describe the molecular interaction can reveal which residues from topIB are more likely to effectively interact with DHBA.Also, this could provide the fragment pattern of DHBA that is more important for the molecular recognition and complex formation.Accordingly, such investigation can also provide a new road to rationally design novel inhibitors with the aim to overcome the drug resistance observed with poison inhibitors. 13Here, we describe computational studies carried out in an attempt to elucidate the mechanism of interaction between DHBA and topIB from L. donovani.

Experimental
TopI from Leishmania donovani deposited in protein databank (PDB ID 2B9S) was used as the model for the in silico studies.Vanadate ion covalently attached to catalytic tyrosine 222 (TYR222) residue from chain B was removed from the X-ray structure and this chain was submitted for reconstruction of missing residues by applying a protocol of stereochemistry correction and minimization. 14The original chain B was then removed from the original structure, and substituted accordingly by that one modeled.
The TYR222 of the reconstructed topIB was set as center for docking of DHBA.The program AutoDock version 4.2 with AutoDockTools v1.5.6 15 was used for this task.Non-polar hydrogens were merged onto the receptor, flexibility of the DHBA was based on the torsions automatically detected, and Gasteiger charges added to both.From the -OH group of the TYR222 was set a grid dimension of 60 Å in all dimensions keeping the other default options of the program.
The best pose predicted by AutoDock was used for building the complex formed by DHBA and topIB to perform a molecular dynamics simulation.The program Desmond v3.8 16 was used for this task.The complex was immersed in water (SPC model) in an orthorhombic box with 10 Å of distance from the complex to the wall in x, y and z dimensions.Neutralization of this system with appropriate counter ions was followed by 150 mmol L -1 addition of NaCl, with the optimized potentials for liquid simulations (OPLS-AA) force field set.The Smooth Particle Mesh Ewald was employed to treat long-range electrostatic interactions with a cutoff of 9 Å for short-range coulombic interactions, and the geometrical constrains of the molecules and bond lengths with hydrogens was treated by the shake algorithm.Minimization of the complex was performed using the steepest descent method until it converged to 1.0 kcal mol -1 Å -1 with a consecutive equilibration phase following a basic multistep protocol.Briefly, this approach included: (i) relaxation of the system by Brownian dynamics using an ensemble NVT(number of particles in the system, volume and temperature) at 10 K with the protein and the ligand restrained with exception for the hydrogens; (ii) a simulation with NVT ensemble and Berendsen thermostat in the same temperature; (iii) simulation with NPT (number of particles in the system, pressure and temperature) ensemble and the Berendsen thermostat and barostat at 10 K and 1 atm of pressure, keeping restraints with exception for hydrogens; (iv) addition of water molecules; (v) NPT ensemble simulation with Berendsen thermostat and barostat, keeping temperature at 300 K and the pressure of 1 atm with restraints for the solution, followed by a rapid constant for the temperature relaxation and slow relaxing of pressure; (vi) an ensemble NPT simulation, with the Berendsen thermostat and barostat at 300 K, 1 atm of pressure, a fast relaxing constant for temperature and a normal relaxing constant for the pressure.After this protocol was established, a molecular dynamics (MD) run of 20 ns at 300 K was performed.Analyses of the entire simulation run were accomplished and the intermolecular interactions identified accordingly.

Docking
Many positively charged residues are located on the surface of the internal channel of topIB, which are essential for anchoring the negatively charged DNA in a right position to allow its catalytic activity (Figure 1).This agrees with the formation of the binary enzyme-DNA complex as observed in 2B9S 17 and described in details by Bosedasgupta et al. 18 The two pockets where DHBA could interact with topIB (Figure 1) showed similar estimative of binding energy as revealed by docking, and both corroborates for the competitive mechanism of inhibition as determined before by Chowdhury et al. 20 The two cavities where DHBA can bind belong to chain A, in a region with high topological and electrostatic complementary (Figures 2a and 2c).

Molecular dynamics
Time-dependent interactions between topIB and the two poses of DHBA identified by docking were evaluated by running a MD simulation of 20 ns for each one in its respective pocket.For the first pocket (binding site 1), the H-bonds formed between the carboxyl from DHBA and the side chains of the arginine 314 (ARG314) from chain A and TYR222 residues were previously identified by docking.However, MD simulation uncovered additional H-bond with the guanidine group from arginine 410 (ARG410) from chain A (Figure 3).
From the root mean square fluctuation (RMSF) graph (Figure 4a) is possible to notice that the octahydroindene carboxyl substructure (Figure 4b) stays in the proximity of the protein in comparison to the initial trajectory of the complex while the rest of DHBA structure (Figure 4d) fluctuates in higher distances from the main chain of topIB along the trajectory.The carboxyl group showed a great change in its position to follow interactions with the side chains of ARG314, ARG410 and TYR222, thus accounting for the most hydrophilic forces keeping the DHBA-topIB     complex.This is corroborated by the analysis of the interaction fraction (Figure 4c), which depicts H-bonds involving almost only this group of DHBA and the residues mentioned before for more than 40% of the simulation run.Then, this part of the molecule can be a convenient fragment to generate new inhibitors of topIB.Also, interactions of the type residue-water-DHBA were present in less than 30% of the trajectory and involved mainly the residues GLU323, ARG410, CYS451, and HIS453.
Reducing the interaction cutoff (i.e., from 40% to 20%, Figure 5) allowed recognizing an additional H-bond, between the hydroxyl group of DHBA and the side chain of LYS319.
The second binding pocket recognized by docking depicts the DHBA molecule surrounded by a hydrophobic environment and H-bond interactions with the residues LYS41, LYS269 and LYS407 (Figures 2c and 2d).Analysis of the interactions that occurred during the MD simulation (Figure 6) indicated that H-bond between the hydroxyl group of DHBA and LYS41 is only briefly established, and its rupture is inevitable because of the high mobility of the loop where this residue locates.This suggests a lack of enough adaptability to attain the interactions with topIB due to a significant conformational changing of this region.
Similarly, the rigidity of the DHBA does not allow following the structural changes of topIB.This is also in agreement with the RMSF plot (Figure 7a), which depicts that the fluctuations of DHBA heavy atoms is practically  also possible to notice that water bridges may occur, and probably this is a result of interactions being lost with time due to the structural mobility of topIB.
The results from Figures 6 and 7 imply that the absent flexibility of DHBA that accounts for its low adaptability to a binding pocket is more critical for the second case.Combining with the significant variation of the topology in this cavity the result is an unlikely DHBA-topIB complex to achieve in comparison to the first one (pose 1, with lower binding affinity evidenced by docking).In particular, the type of molecular interactions observed in Figure 7 for LYS269 and LYS407 could be a result for an enhancing distance between DHBA and topIB as the simulation time increases, most probably as a result of this lack of DHBA adaptability due to the dynamic changes of its initial interaction pocket.As a consequence, we can expect a less stable complex in comparison with the other one.
The competitive inhibition mechanism of DHBA requires that this compound be located in a convenient binding pocket of the enzyme to prevent DNA ligation.Based on this biochemical fact we proposed the following mechanism: (i) the DNA needs to be positioned in the region of its binding recognition located in the channel of topIB, and this is a key step that allows the nucleophilic attack of the residue TYR222 located in the chain B; (ii) the binding of DHBA prevents such specific interaction with DNA, and the catalytic activity of topIB is then blocked.This region shows residues with positive charges, that are complementary to the negative charge of DNA, and it is also the feature for the molecular recognition of DHBA.The two possible binding sites suggested by our model are in agreement with these premises and they were chosen due to the computed values of binding energy, which are equivalent when considering the associated error of this computational procedure.The analysis of the first binding pocket shows  H-Bond interaction between the -OH group from DHBA and topIB is not maintained for a cutoff of 30%, similarly to what was observed for the first binding pocket.According to this simulation, the H-bonds were practically driven by the interaction between carboxyl group from DHBA and the residues LYS407 and LYS269.absent with respect to its initial structure.Although the carboxyl group can also show ionic interactions, it is that interaction with the catalytic residue TYR222 from chain B may be enabled when chains get closer and both the DNA-topIB interaction and accessibility of TYR222 for the catalytic attack are prevented.The second pocket option showed a surface rather less exposed, yet equally able to accommodate the DHBA in a convenient cavity also complementary in terms of electrostatic potential, but with the DHBA positioned far from TYR222, thus the interaction between them was not observed.By running 20 ns of MD simulation with the first pose (−8.92 kcal mol -1 ) it was identified a supplementary interaction between the carboxyl group and the ARG410 residue that contributes to enhance the strength of the forces initially considered for this region.Probably, this unrecognized interaction during the docking is the reason to obtain a lower value in the estimative energy of binding.Importantly, it was uncovered that the octahydroindene fragment practically accounts for the most interactions of DHBA with this region on topIB.Also, it revealed that the interaction between the -OH group from DHBA and the main chain of GLU318 is easily broken and replaced by interaction with LYS319, although during a short period of the trajectory (i.e., less than 3 ns in a discontinued manner).These residues are located in a loop with significant mobility, which contributes for the rupture of such interactions.Apparently, the absence of an additional fragment with higher flexibility in this part of DHBA together with the rigidity of the ring system do not allow for this molecule to follow the conformational changes of this topIB region.This dynamical changing of this pocket from topIB and the low flexibility of DHBA would contribute to further disrupt this complex.Then, modifications in the hydroxyl position of DHBA molecule by introducing convenient flexible and hydrophilic groups would allow exploring interactions more efficiently with other residues surrounding this pocket to render a more stable complex.On the other hand, this finding also emphasized that a simplification of the DHBA ligand could be another interesting approach to provide novel classes of inhibitors as this octahydroindene carboxyl fragment preserves the stronger interactions with topIB.This minimal structure could offer a promising fragment-based approach to generate new lead classes.
Although the second pose allowed establishing important interactions with LYS residues by the docking procedure, MD simulation uncovered the high influence that protein mobility exerts on this pocket, promoting the breaking of intermolecular H-bond interactions involving the hydroxyl group.For this reason, it would be expected a more variable binding pocket in terms of structure and mobility of charged groups, which is less compatible to a rigid structure of DHBA when comparing to the first case.
The more significant structural changes in this binding pocket and the absence of accumulated positive charges in a convenient subpocket like in the first pose would explain why the entire molecule of DHBA on this binding pocket does not accommodate properly.
A redocking procedure was carried out with the snapshot from the trajectory that showed the greatest number of intermolecular interactions between DHBA and topIB regarding the first binding pocket.The results showed that, in reality, a better energy of interaction is possible to achieve.Then, while the structure preparation for docking only consider a local minimum, even when considering the flexibility of the side chains from the residues of the binding pocket, the simulation can reveal conformations that can be adopted and interactions that are disrupted while others are formed, and this is a paramount molecular behavior that is essential for any ligand to complex with its counterpart target.The most important fact here, which is in agreement with previous reported biochemical investigation is the crucial role of the carboxyl group (Figure 8).

Conclusions
The computational models obtained by this study evidenced the most probable binding pocket for the molecular interaction of DHBA with topIB.Also, it revealed the fundamental role established by the octahydroindene carboxylic group of the DHBA to explore key interactions with topIB.The high hydrophobicity of the buried surface in both cavities correlates for the pentacyclic triterpenoid lupane structure to be accommodated, and the potential electrostatic surface depicted that the DHBA molecule fits well into both cavities of the internal channel from topIB identified by docking.These bindings are in agreement with the competitive mechanism of this ligand.However, a more precise study based on time-dependent interaction by applying a MD approach uncovered how interactions can behave along the trajectory of both complexes and their importance for the molecular recognition.While some H-bonds are hold very shortly or easily disrupted others are preserved during a great part of the entire trajectory.Careful analyses of these dynamic changes, and including a redocking procedure with the greatest interactions observed for the simulation trajectory with the first binding pocked allowed the indication of this one as the most probable for DHBA to interact with topIB.
This study contributed to extend our understanding on how the interactions between DHBA and topIB could behave in a time-dependent molecular basis.The structural motion of topIB is expected to have a great influence on the binding of the DHBA molecule.For the first binding pocket important fluctuations were identified, thus allowing a good yet partial molecular adaptability to fit complementary and strong interactions with key residues surrounding this pocket.Therefore, this complex would probably be favored instead of the second one, as the complementary forces of DHBA would be suffice to keep it, yet with a slight overall molecular flexibility observed for this ligand.The introduction of convenient flexible groups at the -OH portion of DHBA could be a strategy to provide new DHBA derivatives for exploring more effectively its interaction with topIB.Additionally, it was also highlighted that a simplification of the DHBA ligand is capable to provide novel classes of topIB inhibitors based on a molecular fragment of DHBA that explores more effectively the interaction with the catalytic subpocket of topIB.According to our simulation study, this octahydroindene carboxyl substructure is the essential part of DHBA that inhibits topIB from L. donovani, and also could be viewed as a suitable fragment from which new molecules can be synthesized.It is worthwhile to mention that terpenes like oleanolic and ursolic acid bear this similar substructure, 21 which also suggests their investigation with topIB from Leishmania species, as similar biological effects of those terpenes are also observed.Most important, this MD study revealed a very relevant chemical fragment that deserves a deep investigation, thus increasing chances towards the development of more promising antileishmanial leads through a fragment-based approach.

Figure 1 .
Figure 1.TopIB from L. donovani (PDB ID 2B9S) after reconstruction of its dimer and the docking result with DHBA: (a) the potential electrostatic surface of chain A depicts the internal part of the channel with many positively charged residues (blue regions) essential for anchoring DNA during the catalytic process; (b) the binding pockets identified by the docking procedure in a closed view of the internal channel surface from chain A and the catalytic TYR222 from chain B in cartoon.From left to right the poses with computed binding energies of −8.92 and −9.25 kcal mol -1 , respectively.Figures generated with Maestro.19 Figure 1.TopIB from L. donovani (PDB ID 2B9S) after reconstruction of its dimer and the docking result with DHBA: (a) the potential electrostatic surface of chain A depicts the internal part of the channel with many positively charged residues (blue regions) essential for anchoring DNA during the catalytic process; (b) the binding pockets identified by the docking procedure in a closed view of the internal channel surface from chain A and the catalytic TYR222 from chain B in cartoon.From left to right the poses with computed binding energies of −8.92 and −9.25 kcal mol -1 , respectively.Figures generated with Maestro.19 Figure 1.TopIB from L. donovani (PDB ID 2B9S) after reconstruction of its dimer and the docking result with DHBA: (a) the potential electrostatic surface of chain A depicts the internal part of the channel with many positively charged residues (blue regions) essential for anchoring DNA during the catalytic process; (b) the binding pockets identified by the docking procedure in a closed view of the internal channel surface from chain A and the catalytic TYR222 from chain B in cartoon.From left to right the poses with computed binding energies of −8.92 and −9.25 kcal mol -1 , respectively.Figures generated with Maestro.19

Figure 2 .
Figure 2. The two best poses for DHBA docking onto topIB and the molecular interactions observed.(a) Topology of the binding site 1 (second best pose in terms of energy) identified by docking showing the main residues and their electrostatic surface potential; (b) the 2D depiction of the binding pocket interactions; (c) topology of the binding site 2 (best pose) identified by docking showing the main residues and their electrostatic surface potential; (d) the 2D depiction of the binding pocket interactions.Figures generated with Maestro, 19 with a 4.0 Å cutoff for the representations in 2D view.H-bond interactions in the tridimensional binding sites are depicted in yellow dashed lines.
Figure 2. The two best poses for DHBA docking onto topIB and the molecular interactions observed.(a) Topology of the binding site 1 (second best pose in terms of energy) identified by docking showing the main residues and their electrostatic surface potential; (b) the 2D depiction of the binding pocket interactions; (c) topology of the binding site 2 (best pose) identified by docking showing the main residues and their electrostatic surface potential; (d) the 2D depiction of the binding pocket interactions.Figures generated with Maestro, 19 with a 4.0 Å cutoff for the representations in 2D view.H-bond interactions in the tridimensional binding sites are depicted in yellow dashed lines.

Figure 3 .
Figure3.The main H-bonds observed between DHBA and topIB during MD simulation carried out with the pose 1 (first pocket), considering 40% of the interaction cutoff.A percentage number higher than 100% is due to the guanidine group from side chain of ARG residues, which donates more than one atom in H-bond for the respective O of the carboxyl group.Then, multiple H-bonds are possible to occur with the same O atom from DHBA.

Figure 4 .
Figure 4. Analysis of the MD simulation: (a) RMSF from the complex topIB-DHBA; (b) ligand atom index of DHBA structure with carbons in black and oxygen in red; (c) a bar plot of interaction fraction (IF) from the simulation indicates that H-bonds between DHBA and the residues ARG314, ARG410 and TYR222 are the prevalent forces that keep this complex (more than 100% of IF means multiple interactions, i.e., with more than one specific H atom from the same residue); (d) the octahydroindene carboxyl substructure according to RMSF is the main part of DHBA closer to the topIB binding pocket (from the RMSF is possible to notice that the atoms 2, 4, 6, 7, 9, 13, 14, 20, 21, 27, 28, 29, 30, 31 of DHBA are in the closer proximity of topIB when considering all the simulation time).

Figure 5 .
Figure 5. Analysis of the MD simulation carried out with the DHBA located in the first binding pocket considering 20% of the interaction cutoff.

Figure 6 .
Figure 6.Percentage of interaction observed between DHBA and topIB for the MD simulation of the second pose (best energy score).H-Bond interaction between the -OH group from DHBA and topIB is not maintained for a cutoff of 30%, similarly to what was observed for the first binding pocket.According to this simulation, the H-bonds were practically driven by the interaction between carboxyl group from DHBA and the residues LYS407 and LYS269.

Figure 7 .
Figure 7. Analysis of the molecular interactions observed for DHBA and topIB considering the pose 2. (a) RMSF plot also suggests low conformational adaptability of DHBA and a distance of carboxyl group that is higher from the residues surrounding it in comparison to pose 1; (b) the interaction fraction of this complex evidences that overall interactions are less expressive during the simulation time for this pose.

Figure 8 .
Figure 8. Result from the procedure of a redocking considering the snapshot from the trajectory of the first binding pocket that showed greater interactions between DHBA and topIB.(a) 3D representation of the best pose identified, which showed an estimative of binding energy according to the docking score of −12.58 kcal mol -1 ; (b) 2D representation of the interactions identified in this redocking procedure.