The constructed systems were minimised in the presence of a harmonic restraint of 5 kcal mol^-1 Å^-2 on the heavy atoms of the ligand and the backbones of the protein for 10 000 steps. Next, the system was equilibrated in a NVT assembly at a temperature of 300 K for 5 ns, then the system was equilibrated in a second phase in a NpT assembly maintaining a pressure of 1 bar for 10 ns with the same restraint. Finally, production in NpT at 310 K and 1 bar of pressure was performed for 300 ns, for a total of three independent replicates (3 x 300 ns) accumulating a total of 0.9 μs.

In the equilibration and production steps the Lavengin thermostat and the Nosé-Hoover Langevin piston were used to control temperature and pressure. Atoms of rigid hydrogens were maintained with SETTLE algorithm [34] in waters, while to solute the RATTLE/SHAKE algorithm was used [35]. Short-range, non-bonded and long-range non-bonded interactions were treated with a cut-off scheme at 12 Å and PME [36], respectively. The r-RESPA multiple time step scheme [37] was used in all cases with 2 fs time integration steps.

For this, the system was prepared and equilibrated with the standard protocol of QwikMD [38], with the only exception that a 40 Å extra padding was added into the Z axis to enable pulling. For pulling a potential bias of the form was added for the original potential through the colvars module [39], getting a modified potential (H ^±). For this, the centre of mass of the ligand ξ _L , was pulling with a speed vof 10 Å ns^-1 and к of 7 kcal mol¹ Å ² during 4 ns, enough to reach the unbound state. During the pulling the C _α atoms were restrained with a force of 5 kcal mol^-1 Å ^-2. A total of 20 independent replicates were performed (total 80 ns). The work was computed for numerical integration using W= v∫₀ ^t1 dt'F(t').

Based on these simulations, the potential of mean force (PMF(ξ),Φ(ξ) was reconstructed by means of Jarnzynski identity as show in the following equations [40]:

If rearranging the equation, it shields the Ф(ξ):

and expanding this equation for the 2^nd order cumulant:

Where Ф(ξ) is the potential of mean force, i.e the free energy profile over the collective variable (ξ), Wis the work done along ξ and β is the inverse of thermal energy (1/k _B T),where k _B is Boltzmann's constant.

On the other hand, to calculate the ΔG ^o _bind standard, we used the following expression:

then:

Where ΔG ^o _bind in kcal mol^-1, k _B in kcal mol^-1 К¹,С^сошр=1 ligand Box^-1 is the concentration used in the simulation in ų and C°=1/1661 Å^-3 is the standard concentration yielding 1 ligand x 1661 Å ³. These last two terms are added to pass from the computational system to the standard system. Here Ф(ξ) ^bound and Ф(ξ) ^unbound correspond to Ф(ξ) the associated and dissociated thermodynamics state (i.e macrostates) respectively.

In order to evaluate changes in the rigidity/flexibility of the bromodomain structure, root-mean-square fluctuations (RMSF) analyses of the backbone were performed. Protein flexibility and dynamics are essential for the molecular recognition process and play a fundamental role in virtually all biochemical processes in living organisms [45]. The results in Figure 2A show that all the replicates present three flexible regions in the protein and present RMSF values from 1 to 4 Å. Furthermore, Figures 2A, 2B and 2C show the structure of the receptor and the location of the flexible regions. R'₁ ranges from PRO-36 to THR-50, R'₂ from ARG-58 to VAL-72 and finally R'₃ from ASP-85 to ASP-95. The aforementioned flexible regions are structurally composed of loops that play a key role in molecular recognition.

Figure 2

The stability and structural changes of the bromodomain and bro-mosporine were evaluated through root-mean-square deviations (RMSD). Figure 3A and 3B indicated different stability values for both; bromodomain presented high values with an average around 5 Å, and bromosporin presented low values 2 Å. The results also indicated that bromosporine presents more than one conformational change. According to the RMSD density plots we noticed three forms of recognition that bromosporin may have in order to interact favourably with the receptor. The obtained RMSD values of the ligand and complex present a correlation because the conformational changes of the complex, specifically the high plasticity of the pocket, cause freedom of movement of the ligand and could generate several rearrangements during the simulation.

Figures 3D, 3E and 3F using a 2D free energy landscape approach revealed that the ligand scans through three different conformational states had similar results obtained based on the RMSD of the ligand. In addition, based on the RMSD values the 3D structure of the receptor and ligand in the initial structure were obtained, as well as the one showing the maximum RMSD. Both structures were aligned for comparison purposes and to observe the regions showing the highest variability. Figure 3C indicates that the greatest variations in the structure, with respect to the initial one, are in the N-terminal and C-terminal regions.

Figure 3

Figure 4 is a set of 3D images showing the bromodomain binding cavity and its associated ligand at three different times. In each image a different conformation of the ligand within the binding cavity is observed. The first image shows the conformation of the ligand at 28 ns with an area of interaction with the receptor determined by the area bounded by the red line. The second image shows a modified conformation of the ligand and the interaction area increases to 193 ns. Finally, the third image shows a new rearrangement of the ligand within the cavity and an increase in the contact area with the receptor. In addition, the images indicate conformational trends of the pocket during the simulation in all replicates, indicating the importance of pocket plasticity.

The relative freedom that can be observed of the ligand during the course of the simulation is due to the constant movement of the loops surrounding the enzyme pocket. The movements allow the ligand to interact with key residues in the complex and start rearranging its atoms until it finds the thermodynamically stable location. The bromosporine at the end of the simulation can fit better because it interacts with the ZA loop which is deeper with in the enzyme pocket than the BC loop [46].

Figure 4

In molecular dynamics simulations, non-bonded energy interactions such as Van der Waals and Coulomb were calculated. Figure 5A and 5B shows the average of the replicas of both energies. Van der Waals interactions contribute more to molecular recognition because they present a strong average energy value of -40 kcal mol^-1compa-red to Coulomb interactions with an average weak value of -27 kcal mo^l-1. Figure 5C shows a Kruskal-Wallis statistical study performed on both interactions to observe the significant differences between each replicate with a value of n=5000. The analyses indicated that there is a significant difference between the replicates in both Van der Waals and Coulomb force.

Figure 5

In Figure 6A, we present the average strength profiles as a function of distance obtained from 20 SMD simulations for the ligand bromosporine bound to the bromodomain complex. At the beginning of each simulation (t = 0 ns), the ligand is in a bound state, interacting with the protein residues at the binding site. We observe the separation force of the bromosporine bound to the bromodomain complex. In the 20 replicates performed, it is observed that there is a peak of 175 pN near 5 Å. indicating that there was bond breaking or interactions within the complex. After the peak, near 10 Å of separation of the complex, the force values drop and do not present any additional energy barrier. In addition, the work of the system was studied. Figure 6B showed a trend in all replicates, where at the beginning of the simulation, the work has an exponential trend up to 25 Å. After that, the work continues to increase gradually up to 30 Å and then the system stabilises, causing the work to not increase and generating an asymptote in the replicas. Finally, we calculate the potential of mean force (PMF) using the Jarzynski equation; Figure 6C shows the PMF profiles as a function of stretching distance between the 20 replicates performed along the coordinate set. The profile consists of two peaks, the first one is near 10 Å which is the highest value with 60 kcal mol^-1 and the second one starts after 20 Å and reaches a value of 40 kcal mol^-1 at a distance of 25 Å. The obtained energy (AG° _bind ) value was 43.79 kcal mol^-1.

Figure 6

Figure 5A shows 3D images of hydrogen bond interactions and hydrophobic interactions with their respective distances. Key amino acids, especially ASN-87 and TRP-93 interacting with bromosporine and surrounding the hydrophobic bromodomain pocket are observed. In addition, Figures 5B, 5C and 5D show the percentage occurrences of hydrophobic, hydrogen bonding and ϖ-stacking interactions respectively of the bromodomain protein amino acids in the three replicates. Hydrophobic interactions indicated that the amino acid TRP-93 is the one that presented the highest frequency, reaching 20 % followed by two Valines at position 36 and 41 that have a similar frequency of 7.5 %. On the other hand, in the hydrogen bond the amino acid ASN-87 was the one that clearly presented the highest and most significant values with respect to the remainder of the amino acids. Finally, ϖ-stacking interactions presented the lowest percentage of occurrence with values less than 1 %. Both of the aforementioned residues provide the major contribution in the binding of the bromosporine inhibitor to the complex. Therefore, these residues may be effective targets for the development of selective bromodomain inhibitors in L. donovani.

Figure 7