UNCONSTRAINED LIGAND EXPLORATION

What it does

The unconstrained ligand exploration and binding site search ready-made script has been designed to predict the interactions and different binding sites between a protein and a ligand. Starting with conformations where the ligand is in the protein’s surface, a non-bias exploration will migrate the ligand around the surface and search for binding sites.

How it does

The control file mixes large and small ligand translation and rotation perturbations. Furthermore, the perturbation direction can be kept for few steps (waitfor parameter). At each iteration the algorithm performs: translation and rotation of the ligand, a move in a given normal modes, side chain prediction (including the side chains of the ligand) and overall minimization.

Changing parameters

We recommend playing around with different parameters that will adapt better the simulation to your specific system.

Translation range: the translation range will control how much the ligand will move in the perturbation step. Usually we combine long [5-8 Å] and small [~1Å] translations. Larger translations will reduce the number of steps to cover the surface but will also eliminate precision. Furthermore, when perturbing very small ligands, like an ion or diatomic molecule, a large translation might result into “tunneling” inside the protein. You will find the ligand crossing the surface in one step, obviously an artifact from wrong parameter selection. Thus, for small ligand the largest perturbation should be around ~3Å. Also, if large translations are used (in particular with a large waitfor (see below), you might move the ligand into the bulk solvent where he might remain most of the simulation.

Waitfor: this parameter will keep the direction (randomly chosen) for the ligand translation for several steps. It is ideal to use with small translation ranges [2-3 Å]. PELE basically interprets the direction search like a cone rather than a line. Thus, we add some variability. Increasing the waitfor, to values of 4-6, will facilitate the exploration of difficult cavities. Basically whenever the perturbation direction points to one of this cavities, the ligand will try to enter it for the waitfor number of steps, before changing into a different random direction. Please note that to use waitfor the steered parameter must be equal to 1.

Tries: in the ligand perturbation step we can try different perturbations and choose the best one. The criterion to select the best one is the total energy evaluations. Thus, for those systems that are electrostatically driven, increasing the number of tries will possibly facilitate the binding cavity search. Increasing the tries will also difficult that the ligand will escape into solution. However, the more tries the less easy will be to enter a difficult channel and to provoke the induced fit since the ligand will prefer to stay outside to avoid steric repulsion.

lcom_con: this parameter will constraint the ligand’s center of mass to its final position after the perturbation, making it more difficult to move along the final minimization. This might facilitate the entrance in small cavities (when combined, for example, with a low tries). It will also make more difficult the tunneling of the surface above described. Use it like:

	lcom_con XX &
	

where XX is the force constant (in kcal/mol) to its current position.

ANM frequency and mode: controlling the frequency and modes in ANM will facilitate the cavity search and binding in some cases. It is highly recommended that before running the ligand dynamics simulations you have some knowledge of the ANM normal modes (which can be done with one of our scripts: “Normal mode exploration”). It is possible that one mode will facilitate the opening and closing of the active site. Then you can tell PELE to sample always this mode (instead of a random selection), in consecutive positive and negative directions, for a number of steps, for example ten steps. This is accomplished by:

	anm_altm_freq 10 &
	anm_altm_type 4 &
	lananm mode XX  &
	

where XX is the chosen mode. Using anm_altm_type 3 you will select modes randomly (the default). Try, for example, to combine lower perturbations with larger frequency or opposite to adjust the best conformational changes to your system. You can also include partially other modes in the backbone perturbation by mixing them:

	lanmanm mix_modes 0.6 &
	

Here we use 60% of the mode we have chosen and 40% a random mix of the reminder modes. Of course, if you select a mode it should be in the list of total number of modes to include (for the top 6 modes use: lanmanm neig 6 &)

Expected outcome

The outcome will include a series of structures with different total and OPLS binding energies. See for example the movies in this section below. In our test cases, there is a very good correlation of the best (lowest) binding energy with the binding site. Typically we use 8-24 processors for ~24 hours to complete a search (depending on the size and ligand complexity). For difficult (closed) active sites (with access channels…), it might be necessary to run ~80 procs for ~48 hours. Alternatively, migrating the ligand with constraints might be an option in these difficult cases.

Caution

If a promoting mode is crucial for the binding (such as in kinases) it is important to accurately sample it, possibly choosing anm_altm_type 4 &. Also, when you look at the trajectories for the ligand entrance, be careful with the tunneling events described above. You might need to print every accepted step for small ligands to check that no tunneling is happening.

References

We are currently writing a larger benchmark manuscript. However initial applications have been already published.

1. Lucas, F.M. & Guallar, V. Single vs. multiple ligand pathways in globins: A computational view [Biochimica et biophysica acta, In Press (2013)].