Restraints and their corresponding free energies of attachment and release

The BAT code employs essentially the same set of restraints as previously described for the APR method²⁴. Both the protein and the ligand are subject to two types of restraint. One type restrains the position and orientation of the molecule in the frame of reference of the simulation box. These translational and rotational (TR) restraints are constructed with added length, angle, and dihedral potential energy terms defined between atoms of the molecule (protein or ligand) and three dummy atoms, termed N1, N2 and N3, whose locations are fixed in the lab frame (Fig. 2). These keep the protein and ligand restrained relative to the simulation box, and thus to each other²⁴. The other type of restraint is applied to internal degrees of freedom of the protein and ligand, so these limit conformational freedom. They are designed to reduce fluctuations during the calculation of the transfer free energy, and thus help convergence. This benefit must be weighed against the added computational cost of computing the attach and release free energies (Eq. ³). Note that the final free energy of binding should be the same, aside from numerical error, with or without the use of conformational restraints, and their use is, at least in principle, optional.

The TR restraints on the protein comprise harmonic potentials applied to one distance (D2), two angles (A3, A4), and three dihedrals (T4, T5, and T6), between the dummy atoms and three protein atoms (P1, P2, P3), which are termed the protein anchors²⁴. These protein TR restraints are present and active during the equilibration phase and throughout most of the free energy calculation. They do not affect the protein’s conformational distribution, so there is no need to compute the free energy of attaching and releasing them. The force constant for the restraint on D2 is set in the BAT.py input file with the rec_distance_force variable, and the force constants for A3, A4, T4, T5, and T6, are set with the rec_angle_force parameter. The reference values of these restraints are taken from the starting conformation. (See user manual for details).

The conformational restraints on the protein comprise three harmonic distance restraints among the three protein anchors (P1, P2 and P3), to reduce TR fluctuations and coupling between TR motions and conformational changes. In addition, harmonic restraints may be applied to the backbone ϕ and ψ angles in a user-selected range of protein residues (Fig.  4), in order to keep this section relatively rigid when transferring the ligand from binding site to bulk, particularly when using the APR method. Backbone ω angles are not restrained, as these are already rather rigid, due to the double-bond character of the peptide bond. This option is activated in the BAT.py input file using the rec_bb variable, with the chosen residue range defined by the bb_start and bb_end parameters. The spring constants for the protein distance and backbone dihedral conformational restraints are specified with the rec_discf_force and rec_dihcf_force variables, respectively.

(Top) Second BRD4 bromodomain with the restrained backbone section (part of the ZA-loop) in yellow, the rest of the protein in blue, and the ligand colored by element. The structure is from the 5uf0 cocrystal structure. (Bottom) Protein backbone showing Ramachandran torsion angles. The optional backbone restraints in BAT.py are applied to only ϕ and ψ angles.

The free energies of attaching, (ΔGp,att), and releasing (ΔGp,rel) the protein conformational restraints are calculated in the absence of the protein TR restraints (top and bottom processes in Fig.  2), using a number of simulation windows with Nw values of the restraint spring constants, between zero and its full value, corresponding to Nw windows. At each window, i, the spring constant for a given restraint, r, is defined by:

where kir is the spring constant, kfr is its full value defined in the input file, and attach_rest(i) is the multiplying factor associated with each window. These factors are defined by the attach_rest array in the BAT input file, and should go from 0 to 100. The same factors are used for all of the attach and release calculations for the protein and the ligand, which were determined and optimized in our previous APR study on the BRD4(1) protein²⁴. Note, however, that the free energy of the ligand TR release is computed semi-analytically rather than with simulations. The free energies are obtained from the trajectories of all windows, using the multistate Bennett acceptance ratio (MBAR) method⁴⁸, according to the equation:

Here the subscripts i, j, and k index the simulation windows; n indexes the Nj samples from window j, each with coordinates rn, and Nk is, analogously, the number of samples from window k; Gi and Gk are the free energies of windows i and k, respectively; β=1/kBT; Ui(rn) is the potential energy from the restraints defined in window i acting on the coordinates rn, which correspond to the nth sample from window j. Thus, Ui(rn) is given by

where R is the number of restraints being attached or released; xnr is the value in sample (or frame) n from window j of the internal coordinate (distance, angle, or torsion) corresponding to restraint r; and kir and x0,ir are, respectively, the spring constant and equilibrium value for the rth harmonic restraint in window i. The program pyMBAR⁴⁸ is invoked by BAT to solve Eq. (⁶) self-consistently for the free energies across the windows. The free energy difference between the initial (unrestrained) and final (fully restrained) states is then obtained directly from this set.

Like the protein, the ligand is subject to harmonic TR and conformational restraints²⁴. The TR restraints again comprise a distance, D1, two angles, A1 and A2, and three torsions, T1, T2, T3, defined relative to the three fixed dummy atoms N1, N2 and N3. The spring constant for D1 is set in the BAT.py input file using the lig_distance_force variable. For the angle and dihedral TR restraints, the spring constant is defined by the lig_angle_force parameter. The reference values of these restraints are taken from the initial coordinates. The ligand conformational restraints include harmonic potentials on the three distances between its anchor atoms L1, L2 and L3 (Fig. 2). In addition, essentially all dihedral angles are also restrained to make each ligand rigid and thereby accelerate convergence. For simplicity, torsions within rings are not excepted from the set of restraints, although they are not always necessary. The BAT.py script automatically assigns these restraints for each ligand. It uses the ligand’s AMBER parameter/topology (prmtop) file (Fig. 5) to identify all proper dihedral terms not involving a hydrogen atom, and assigns a restraint to one arbitrarily chosen dihedral term for each central bond. The spring constants for the ligand’s internal distance and dihedral restraints are set in the BAT.py input file via the lig_discf_force and lig_dihcf_force parameters, respectively, and their reference values are taken from the starting coordinates. Fig. 5 illustrates the assignment of 14 conformational dihedral restraints for the ligand from cocrystal structure 5uf0.

Example of ligand dihedral restraints, for PDB ID 5uf0. (Left) Section of AMBER parameter/topology file listing all the ligand dihedrals that do not include a hydrogen atom. Each row lists two torsions in terms of five indices; the first four map to specific atoms and the fifth maps to the associated force field parameters. Dihedrals restrained in the BAT procedure are highlighted in purple font, with redundant ones in black and improper dihedrals in red. (Right) Ligand from 5uf0 with restrained torsions highlighted with purple bonds. Cyan: carbon. White: hydrogen. Red: oxygen. Blue: nitrogen.

The free energies of attaching and releasing the ligand restraints may be separated into conformational and TR parts:

During the attachment stage (att), the ligand is in the binding site of the restrained protein. The conformational restraints are applied first, yielding the free energy change ΔGl,conf,att for making the ligand essentially rigid. The TR restraints are then applied, yielding the free energy change (ΔGl,TR,att) for restraining the ligand in the binding site. During the release stage (rel), ΔGl,conf,rel is computed with the ligand in a separate simulation box with no TR restraints present. The values of ΔGl,conf,att, ΔGl,TR,att, and ΔGl,conf,rel are calculated the same way as the protein conformational restraints, using MBAR (Eqs. (⁵)–(⁷)), with simulation windows having intermediate values of the harmonic spring constants, also defined by the attach_rest input array. The final term in Eq. (⁹), ΔGl,TR,rel, is calculated by numerical quadrature of the following integral, which is based on Euler angles and spherical coordinates:

Here C∘ is the standard concentration, 1 M = 1/1661Å3, and r, θ and ϕ are the distance D1, angle A1, and dihedral T1, respectively²⁴. In the last term on the right, which integrates over ligand orientation, Θ is the angle A2, Φ is the dihedral T2, and Ψ is the dihedral T3. These are three Euler angles which define the orientation of the ligand in space. The harmonic potential applied to the distance r has the form:

with kd the lig_distance_force spring constant and r0 its reference value. A similar expression is used for restrained angles and dihedrals:

with a a given angle/dihedral, ka the lig_angle_force spring constant, and a0 its reference value. The value of r0 is the reference distance, D1, from dummy atom N1 to ligand atom L1, in the bound state; this is always sets to 5.00 Å by construction (“Anchor atoms and dummy atoms”).