Computational Studies of the Effect of Shock Waves on the Binding of Model Complexes.

We have simulated effects of a shock wave in water that would result from the collapse of a cavitation bubble on binding in model complexes. We have considered a benzene dimer, a pair of uracil molecules, a complex of fragments of the X-linked inhibitor of apoptosis and caspase-9, and a fragment of c-Myc oncoprotein in binding with its dimerization partner Max. The effect of the shock waves was simulated by adding a momentum to a slab of solvent water molecules and observing the system as the slab moved and caused changes. In the cases of the small molecular pairs, the passage of the shock waves lead to dissociation of the complexes. The behavior of the protein systems was more complex, yet significant disruption of the binding and geometry was also observed. In all the cases, the effects did not occur during the immediate impact of the high-momentum solvent molecules, but rather during the expansion of the compressed system that followed the passage of the waves. The rationale of the studies was in attempting to understand the strong effects that irradiation with a low-intensity ultrasound can have on biomolecular systems, because such ultrasound irradiation can cause cavitation bubbles to be produced and collapse, thus leading to local shock wave generation. The long-term objective is to contribute to future design of synergetic ultrasound and chemical drug strategy of protein inhibition.

Introduction

Low-intensity ultrasound is known to produce biomedical effects that cannot be explained by energy density of the sonic irradiation as such. Examples of these phenomena include synergy with chemical drugs[1] and destruction of tumor cells.[1b,2] It has been hypothesized that such effects proceed via cavitation and the subsequent collapse of the cavitation bubbles.[2,3] The low-pressure phase of the ultrasound irradiation leads to a drop in pressure, which leads to the emergence of microbubbles. The bubbles collapse during the higher-pressure phase, and local shock waves are produced as a result, with pressures reaching tens of kilobars.[4] These waves travel through the medium, interact with membranes and other parts of the tissue (including tumor cells) and cause damage that results in death of the cells. Targeted use of these effects lead to applications in minimally invasive microsurgery and nanosurgery.[5]

The use of shock waves, regardless of the origin, has also been suggested in increasing the permeability of cell membranes and delivery of molecular agents into cells.[6]

Studies of the shock waves and interactions of them with membranes have been undertaken. These projects included some computational modeling,[7] but a comprehensive picture of the microscopic damage produced by the shock waves resulting from low-intensity ultrasound-induced cavitation is lacking, especially at the all-atom level.

We have been studying effects of the shock waves on molecular dimers, including protein–ligand ones. The main questions to ask are whether such waves can be employed to destroy certain complexes and, thus, to promote alternative binding (such as binding of inhibitors introduced as chemical drugs). We published an article that described dissociation of the X-linked inhibitor of apoptosis (XIAP) with the caspase-9 protein, based on a crude semimacroscopic model of the process.[8] The conclusion of that work was that such complexes can dissociate as a result of irradiation with waves of low frequencies (as ultrasound frequencies are compared to characteristic frequencies of molecular vibrations), because of the high amplitude of these waves. In essence, the molecular complex is compressed in the high-pressure phase, thus excess mechanical energy is stored in it, and this energy can be channeled into dissociation of the complex in the following low-pressure phase, so that the frequency as such is not relevant, as long as there is sufficient time for the dissociation process. Thus, it should be possible to affect biomolecular processes, depending on protein binding, by applying low-intensity ultrasound and causing the resulting cavitation and microbubble collapse.

In the present work, we study effects of shock waves on a number of model bound systems at a more-detailed all-atom level with the aid of molecular dynamics (MD). It has been shown that the influence of shock waves on molecular systems can be simulated by supplying additional linear momentum to a slab of solvent water molecules.[7,9] This is the method we have applied to studying the disruption and dissociation of bound molecular complexes, as a result of the passage of the microbubble collapse-induced shock waves.

We have selected four model systems for this study. Benzene and uracil pairs in water serve as an example of nonpolar but polarizable π-electron system dimer and a dimer (benzene) with a more-pronounced electrostatic component (uracil). The fragment of the XIAP in complex with the caspase-9 represents a protein binding case relevant for apoptosis (which is known to be induced in some cavitation and the following shock wave cases). Finally, the c-Myc–Max zipper pair represents this different type of protein interactions, as well as stands for a very valid cancer-related target.[10] Thus, the reported research is viewed as a step in studying the cavitation-induced shock wave damage that can be imposed on various molecular complexes. In future work, such a damage can be applied to change the naturally occurring binding affinities and binding of inhibitors, and, thus, promote or restrict protein functions.

The rest of the article is organized as follows. The procedures are described in the second section (“Methods”). Results and discussion are presented in the third section (“Results and discussion”), followed by the fourth section (“Conclusions”), which summarizes our findings.

Methods

General Simulation Setup

The simulations were carried out with the molecular dynamics (MD) module of the TINKER software package.[11] OPLS-AA force field[12] was employed. Five systems were simulated: pure water (for reference), water with the benzene dimer, water with the uracil dimer, water with a fragment of X-linked inhibitor of apoptosis (XIAP)-bound complex with the 9-caspase, and water with solvated fragment of the c-Myc–Max zipper-type complex.

In each case, simulations were started by equilibrating a cell of ca. 8780 water molecules (for the benzene and uracil simulations) or ca. 17 200 water molecules (for the proteins) and the solutes (except for the case of pure water). The cells were elongated in the X-direction, along which the shock wave would propagate. Periodic boundary conditions were used. The initial equilibration was done with the NPT (isobaric–isothermal) ensemble at 25 °C and 1 atm. In all of the cases, 8.5 Å cutoffs were employed for nonbonded interactions with the interaction energies feathered over the last 0.5 Å of the distance. The time step was equal to 1 fs. The length of the equilibration was at least 500 ps for all of the systems.

Having equilibrated the system, we proceeded to simulate the effect of the shock wave. Approximately 20% of all the water molecules (the solvent slab at the left end of the simulation cell) were given additional linear momentum. The value of the momentum was set to I = 40 mPa s, which is a value similar to that suggested in the literature.[7] The additional velocity V in the positive X-direction for all the affected atoms was calculated in accordance with eq 1:[7]Here, m is the mass of one water molecule, and A and N represent the area of the solvent slab with the modified velocities of the atoms and the number of the water molecules in the slab, respectively. The resulting values of the additional velocities are ca. 5–8 km/s or 50–80 Å/ps, which has the same order of magnitude as the typical velocities of a shock wavefront, but it has been pointed out that it is the momentum and not the velocity that is important in this procedure.[7,9]

This initial configuration of the shock wave propagation in the equilibrated cell is shown in Figure 1a. The velocities have already been set, but the density is not yet distorted. The solute molecules (in this case, a fragment of the c-Myc in complex with Max) is visible near the center of the cell. The MD simulations proceed further in the NVE ensemble. The periodicity is preserved only in the Y- and X-directions, and the value of the period in along the X-axis is set to infinity. The time step is reduced to 0.2 fs at this stage of the simulations.

Figure 1

Three stages of the simulations: (a) a slab of water molecules (one-fifth of the system on the left) is given an additional liner momentum; (b) a region of increased density has been created and is approaching the solutes; (c) the front of the shock wave has passed the location of the solutes.

As shown in Figure 1b, an area of increased solvent density is built and moves toward the solute. The front then passes the dimer (Figure 1c) and continues toward the right-side end of the cell. Eventually, the solvent in the cell dissipates (as there is no periodic boundary conditions in the X-direction), and the model loses its physical relevance, such that the simulations must be stopped. We define the time when the simulations are completely relevant to be from the beginning to the moment when the density of solvent around the solutes drops to the level below the original average density. After this moment, the model still retains some of its qualitative relevance for a while, especially in the sense of further development of tendencies that emerged during the previous part of the simulations, but the exact correspondence to the physical picture is no longer present.

Molecules and Protein Fragments Used in the Simulations

The simulations of the benzene and uracil pairs employed systems shown in Figure 2. Typical configurations after equilibration of the solvated dimers are given in the figure.

Figure 2

Typical configurations of the benzene (a) and uracil (b) dimers. Water molecules are not shown, for the sake of clarity.

The protein fragments are shown in Figures 3 and 4. Figure 3 shows the simulated fragment of the XIAP bound to the caspase-9. We included one of the main binding parts of the complex. This part consisted of residues GLY335–HIS346 in chain A and LEU384–PHE404 in chain B of the Protein DataBank (PDB) file system 1NW9. Figure 4 shows the fragment of the c-Myc–Max complex. In this case, residues ALA946–LEU960 from chain A and ALA246–GLN261 of chain B from the PDB system 1NKP were used.

Figure 3

Simulated fragment of the X-linked inhibitor of apoptosis (XIAP) bound to the caspase-9.

Figure 4

Simulated fragment of the c-Myc–Max complex.

The OPLS-AA Force Field

A very brief description of the OPLS-AA force field is given here. The total energy (Etot) was calculated by adding the electrostatic interactions (Eelectrostatic), the van der Waals energy (EvdW), harmonic bond stretching and angle bending (Estretch and Ebend, respectively), and the torsional energy (Etorsion) terms:The electrostatic term is calculated as a sum over all atom–atom pairs i ≠ j:where the factor f is set to zero for 1,2- and 1,3-pairs (atoms which belong to the same valence bond or angle), to 0.5 for 1,4-interactions (atoms in the same dihedral angle) and to 1.0 for all the other pairs.

The standard Lennard-Jones formalism is employed for the van der Waals energy:

Geometric combining rules are applied as (ε = (εε)1/2,σ = (σσ)1/2). Harmonic bond stretches and angle bending were used and the torsional term is obtained as a Fourier series:

Results and Discussion

Simulations of Pure Water

Pure water was simulated to have benchmark data on propagation of the shock wave in the absence of any solutes. After the initial equilibration, the size of the simulation cell was 345.5 Å × 26.1 Å × 26.1 Å, the first side being along the X-axis. After we have increased the momentum of the water slab in the ca. 69.1 Å leftmost part of the system, we immediately divided the entire observed space into 120 bins in the X-direction. The middle ∼40 bins or cells (∼40–80 cells) were the ones in which water molecules were located at that time. The additional cells on the right side and left side were added to monitor any molecules that would escape the boundaries of the original system. We then followed the density of the system (number of the water molecules in those cells) and average velocity of molecules in the cells to obtain the distribution of density and velocity in the system. The same method was used to monitor the system density and velocity distribution in the simulations that included solvated molecular dimers.

Figure 5 shows the evolution of the density distribution in the pure water system after the increase of the linear momentum of the solvent slab at the leftmost part of the total simulation cell. At the beginning, the density is roughly uniform (the black line). The moving water molecules then create a spike in the density in the left part of the system (as illustrated for the 1 ps time by the red line). The spike moves to the right and reaches the middle of the system between 2 and 3 ps of the simulation time. Once the front of the shock wave reaches the rightmost edge of the system, density bins beyond those filled in the beginning become populated, as can be seen for the density curves at 5, 6, and 7 ps. Finally, the density inside the simulation cell falls to its initial equilibration value at ca. 7 ps. Thus, we can say that results for the first 7 ps after the generation of the shock wave are physically more relevant than those produced for the times after 7 ps, and the latter should only be considered in qualitative and approximate descriptions of the process.

Figure 5

Evolution of density distribution of the system in the pure water simulation. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

Figure 6 shows distributions of the water molecule velocities at various moments after generation of the wave. The initial additional velocity of ∼50 Å/ps is added to the leftmost fifth of the simulated system (the black line). It can be observed that the height of the maximum peak decreases as the wave propagates to the right and its momentum dissipates, but only until the last three curves (t ≥ 5 ps) are reached, as at that point the front reaches the edge of the system and the dissipation at the front does not occur any more.

Figure 6

Evolution of velocity distribution of the system in the pure water simulation. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

Simulations of Solvated Benzene Dimer

Results of simulations of the benzene dimer are given on Figures 7–10. The size of the equilibrated cell was 343.6 Å × 25.8 Å × 25.8 Å, not very different from the size observed for the pure water system. Figures 7 and 8 demonstrate the evolution of the density and velocity distribution with the passage of time after the shock wave is generated. The behavior of the system in these respects is essentially the same as in the case of the pure water simulations described above.

Figure 7

Evolution of density distribution of the system in the solvated benzene dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

Figure 10

Evolution of the distance between the centers of the solvated benzene molecules. The increase in the momentum of a slab of the solvent occurs at t = 500 ps.

Figure 8

Evolution of velocity distribution of the system in the solvated benzene dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

The graph presented in Figure 9 displays the behavior of the intermolecular benzene–benzene interaction energy. The time is counted from 500 ps before the generation of the wave; thus, t = 500 ps corresponds to the moment when linear momenta of the water molecules in the leftmost fifth of the simulation cell are increased.

Figure 9

Evolution of the binding energy in the solvated benzene dimer simulations. The increase in the momentum of a slab of the solvent occurs at t = 500 ps.

It can be seen that the benzene–benzene interaction energy is well-equilibrated and rather stable at slightly higher than −2 kcal/mol. It then rises sharply at ca. 1–2 ps, after the shock wave starts to propagate (because the wave needs some time to reach the location of the dimer). Then, it eventually settles at zero; the process is apparently almost completed by ca. 7 ps after the beginning of the wave propagation.

Additional information can be gained from examining the data in Figure 10. This graph shows how the distance between centers of the benzene molecules changes with the passage of time. It oscillates slightly around ca. 5 Å once the system is equilibrated. Then it decreases slightly at 1–2 ps after the generation of the wave and then starts to increase until the dimer dissociates. The distance reaches ca. 8 Å at 7 ps after the wave starts propagation, so the dissociation is decided by that time.

Given the above information, we can reconstruct the following physical picture. The dimer is equilibrated in the bound configuration. After the front of the shock wave reaches its location, the benzene molecules are pushed closer together by the rapidly increasing pressure (this is why the intermolecular distance decreases and the interaction energy shows a high peak). The pressure then decreases, and the potential energy stored in the compression is released, pushing the benzene molecules apart and leading to complete dissociation of the dimer. Therefore, the actual process of the dimer destruction takes place not at the moment of the shock wave arrival, but later, when the pressure dissipates.

Simulations of Uracil Dimer

Results of these simulations are given in Figures 11–14. The size of the equilibrated cell was 345.6 Å × 25.9 Å × 25.9 Å.

Figure 11

Evolution of density distribution of the system in the solvated uracil dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

Figure 14

Evolution of the distance between the centers of the solvated uracil molecules. The increase in the momentum of a slab of the solvent occurs at t = 500 ps.

Evolution of the binding energy in the solvated uracil dimer simulations. The increase in the momentum of a slab of the solvent occurs at t = 500 ps.

Once again, the behavior of the density and velocity distributions (Figure 11 and 12, respectively) is very similar to that for the pure water simulations and is not very different from the result for the benzene dimer.

Figure 12

Evolution of velocity distribution of the system in the solvated uracil dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

The evolution of the uracil–uracil interaction energy and distance between the centers of the rings is also qualitatively quite similar to that for the benzene dimer. However, there are some differences. The average interaction energy for the equilibrated state is greater in magnitude and is between −6 kcal/mol and −7 kcal/mol. The dip in the intermolecular distance is somewhat more pronounced. Furthermore, there is a partial decrease in the intermolecular distance at ca. 10 ps after generation of the shock wave, and there is a corresponding decrease in the interaction energy, even though the dimer does dissociate eventually. The reason for this partial coming back of the dimer is probably in the fact that the uracil–uracil interaction is stronger than the benzene–benzene one. As the molecules start to move away from each other, there is some resistance to this motion from the surrounding solvent that has to rearrange. While the same phenomenon occurs for the benzene–benzene dimer, the uracil–uracil attraction is greater, and so it manages to temporarily and partially restore the dimer before it dissociates completely.

However, the qualitative picture is the same, and it appears to be clear that dimers of small molecules dissociate once the wave of the locally very high pressure passes and the stored potential energy can be released in the motion of the molecules away from each other.

Simulations of XIAP–Caspase-9 Dimer Fragment

The size of the simulation cell for the protein dimers was different from that for pure water and small molecular dimers. The dimensions of the equilibrated cell for the XIAP–caspase-9 complex were 225.3 Å × 45.6 Å × 45.6 Å. Qualitatively, the behavior of the density and velocity distributions was similar to that in the previous cases, except that the smaller size along the X-axis lead to a somewhat shorter time frame of the passing of the shock wave through the system. (See Figures 15 and 16.)

Figure 15

Evolution of density distribution in the solvated XIAP–caspase-9 dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

Figure 16

Evolution of velocity distribution in the solvated XIAP–caspase-9 dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

However, the evolution of the dimer was more complex in this case (and the next protein dimer case). The size of the molecules is much larger than that of the benzene or uracil systems; thus, complete dissociation of the complexes within the short time frame of the simulations is not possible. The evolution of the XIAP–caspase-9 interaction energy is shown in Figure 17. After the peak of ca. 140 kcal/mol, which corresponds to the compression created by the passing shock wavefront, the energy oscillates and seems to settle at a value near −60 kcal/mol, which is not greatly higher than the equilibrium value of approximately −80 kcal/mol for this complex. However, the physical picture is rather far from the complex returning to its original state. The graph in Figure 17 shows the interaction energy between the molecules and not the total potential energy of the solute. In the cases of the benzene and uracil dimers, the difference between the two was not fundamentally significant, since these solute molecules are small and do not have any significant possibility of conformational change. However, the protein fragments simulated in this part of the project have a much greater ability to experience changes in the intramolecular conformations. The equilibrated value of the intermolecular energy is approximately −80 kcal/mol, with the total solute energy immediately before the increase in the linear momentum occurs being −445.8 kcal/mol. At the end of the simulations (t = 520 ps in Figure 17), the intermolecular energy is −51.4 kcal/mol, but the total potential energy of the solute is +149.6 kcal/mol. Therefore, deformation of the solutes permits the storage of additional energy: (149.6 + 51.4) – (−445.8 + 80) = +566.8 kcal/mol. This energy must dissipate eventually, and it is probable that this can lead to either direct dissociation of the dimer or to a serious deformation that makes binding less favorable.

Figure 17

Evolution of the binding energy in the solvated XIAP–caspase-9 simulations. The increase in the momentum of a slab of the solvent occurs at t = 500 ps.

Changes in the geometry of the protein fragments can be summarized as follows: the longer chain B is clearly undergoing the most significant conformational perturbations. It is essentially denaturing, and the binding becomes weak and limited to a relatively small fragment, with the remaining part of the chain being pulled away from chain A and solvated independently. A more complete understanding of this process would likely require simulation of a bigger part of the protein and a longer simulation time, but the general character of the processes that are set in motion by the passage of the shock wavefront is shown by the reported data.

Simulations of c-Myc–Max Pair Fragment

The dimensions of the equilibrated cell for the c-Myc–Max complex were 227.6 Å × 45.5 Å × 45.5 Å. The general behavior of the shock wave density (Figure 18) and velocity (Figure 19) distributions was the same as in the previous cases.

Figure 18

Evolution of density distribution in the solvated c-Myc–Max dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

Figure 19

Evolution of velocity distribution in the solvated c-Myc–Max dimer simulations. The time positions in the legend are given relative to the moment of the increase in the linear momentum.

The solute–solute intermolecular binding energy of this zipper-type complex was much lower at the start (see Figure 20), with its equilibrated value being approximately −290 kcal/mol. The compression caused by passage of the shock wave leads to an increase (reduction of magnitude) in this energy of ∼160 kcal/mol. At the same time, the solute–solute energy is more negative 20 ps after generation of the wave. At the end of the equilibration, the intermolecular interaction energy is −297.6 kcal/mol. At the same time, the overall solute potential energy is −696.5 kcal/mol. At the end of the 520 ps simulations, the interaction energy of the dimer is −436.5 kcal/mol, while the total energy of the solutes is −382.9 kcal/mol. Thus, the additional potential energy stored as a result of the solute deformation is (−382.9 + 436.5) – (−696.5 + 297.6) = 452.5 kcal/mol, which is somewhat smaller than in the XIAP–caspase-9 case but qualitatively similar to that result. Once again, one can expect this energy to play a role in either direct dissociation of the dimer or in creating altered binding conditions for the complex.

Figure 20

Evolution of the binding energy in the solvated c-Myc–Max simulations. The increase in the momentum of a slab of the solvent occurs at t = 500 ps.

Evolution of the dimer geometry is such that, toward the end of the simulations, chain A is stretched and chain B is compressed. There is no complete or near-complete dissociation of the complex, but apparently this conformational change leads to a possibility of tighter intermolecular contact and higher magnitude of the binding energy, while the deformation in the chains can still lead to dissociation of the dimer at a later time.

A general comment should be made at this point. One must be careful when assessing applicability of ultrasound to potential anticancer medical drugs. It has been pointed out[13] that ultrasound and resulting cavitation and microbubble collapse can actually reduce the efficiency of some drugs and can disrupt inhibitor binding or, in some cases, even lead to drug decomposition. This is why it is important to evaluate not only the effect of the cavitation and shockwaves on the complex that is to be disrupted on purpose, but also on the complex with the inhibitor that should be promoted.

Moreover, ultrasound irradiation has been shown to enhance the efficiency of some drugs, for example, to increase the damage to cancer cell DNA caused by cisplatin.[14] At the same time, ultrasound-induced dissociation can lead to the opposite result. Therefore, it is important to consider the balance between the two trends and to simulate effect of the ultrasound irradiation and shockwave propagation on stability of all the complexes relevant in the context of the problem at hand.

Conclusions

Simulations of the effects of passing a shock wavefront through areas containing dimers of small molecules and protein fragments have been carried out. The shock wave was emulated by adding extra liner momentum to some of the solvent water molecules. The wave could have various origins, but the main motivation for the reported research was in studying the results of the passage of shock waves generated in the collapse of cavitation bubbles caused by irradiation of low-intensity ultrasound.

In all of the cases, the dimers were severely affected by the wavefront. Dimers of smaller molecules dissociated within picoseconds after the creation of the shock waves, and protein pairs, were significantly distorted, although complete dissociation did not take place within the time of the simulations. At the same time, hundreds of kcal/mol of deformation energy were stored in changing of the protein geometries, and it could be reasonably expected that this energy would be later released, resulting in direct dissociation of the dimers or in dissociation following severe changes in the binding environment.

It should be emphasized that the dissociation or deformation tendencies were developed not during the passage of the wavefront as such, but rather later, when the local increase of pressure became less pronounced and expansion of the systems became possible. Thus, the current results are in qualitative agreement with the application of a crude dimer model.[8]

Since different complexes undergo different changes as the shock wavefront passes, it is hoped that further studies of specific complexes and their behavior under similar conditions will permit one to predict and design cases for synergy of medical drugs (inhibitors) with low-intensity ultrasound irradiation.