I'm impressed. The methodology is far more advanced than I thought was possible at this time. They are clearly using some advanced simulation techniques...
In silico drug discovery in 3D
Snippet from the methods section...
In Silico Modeling. The predict algorithm and methodology were recently reported elsewhere (14, 21). Here, we shall give only a brief overview of the method. predict is a de novo GPCR modeling methodology that combines the properties of the protein sequence with those of its membrane environment, without relying on the rhodopsin (or bacteriorhodopsin) x-ray structure. The predict algorithm searches through the receptors' conformation space for the most stable 3D structure(s) of the TM domain of the GPCR protein within the membrane environment. To ensure that the final model represents the most stable conformation, the method simultaneously optimizes several thousand alternative conformations of the receptor (denoted as “decoysâ€). The final model is accepted only if it is significantly more stable than the majority of the decoys.
The algorithm solves this complex search and optimization task efficiently by using a reduced representation of the protein-membrane system, which balances computational efficiency and accuracy. In this representation, each side chain is represented by two to four virtual atoms (22), allowing for an efficient search through rotamer space (23) while retaining a low dimensional representation of the system. The reduced representation is expanded to an all-atom model toward the end of the modeling process. The algorithm also takes into account, in a simplified way, the presence of the membrane environment and the different character of the membrane lipophilic core and the polar head group region. These components, as well as various protein–protein interactions, are introduced into the modeling procedure by means of the energy function (see below).
Following are the main steps in the predict algorithm (details in refs. 14 and 21). First, an extended sequence around the TM domain (but longer than it) is identified by using known methods, such as a combination of hydrophobicity and sequence conservation patterns (24). A 2D grid is used to construct thousands of alternative packing geometries that cover the protein's conformation space. Hydrophobic moments are used to rotate the helices so that the bundle presents a hydrophobic surface toward the membrane. Each decoy then undergoes a series of optimization steps, including optimization of helix orientation, helix vertical alignment (relative to the other helices and relative to the membrane/water boundary), helix position, and helical tilt angles. Each change in any of these factors is followed by stochastic simulated annealing optimization of the side-chain rotamers in the vicinity of the change. Optimization of the tilt angles is attempted along all possible three- and four-helical arches. Finally, the optimized models are ranked according to their predict energy score. Models with energy scores significantly lower than other decoys are considered solutions. Similarity clustering is then used to reduce the number of solutions. The lowest energy representative of the largest cluster is the final model. This model, still in a reduced representation, is then expanded to an all-atom model maintaining the specific side-chain rotamers that were optimized by predict.
The energy function used for optimizing the 3D conformations and for scoring the models includes two terms, an intraprotein residue–residue interaction term, and a single-residue term, reflecting its interaction with the membrane,
where Resi and Resj are the two interacting residues, and their interaction energy Eint(Resi, Resj) is defined as
where εij are the Miyazawa and Jernigan (25) contact energies between residue i and j, fij is a distance function with the general shape of a “soft†Lennard–Jones potential (a 6–4 potential in agreement with the multiatomic nature of the virtual “atoms†used in this representation, unlike the atomistic 12–6 Lennard–Jones function), λarom is an aromatic-clustering factor highlighting aromatic–aromatic interactions, λcat is a cation–π interaction factor reflecting what is recognized as an important noncovalent binding interaction in α-helical peptides (26), and λpolar is a polar–polar interaction factor that emphasizes their contribution in agreement with studies that point to specific polar interactions implicated in driving TM helix association (27, 28), especially in the hydrophilic core of GPCRs. The protein–membrane interaction term Emembrane(Resi, Zi) is a function of the chemical character of Resi and its position Zi in the direction normal to the membrane plane. The value of Zi determines whether the residue is interacting with the lipid core or with the polar head groups, adjusting the interaction accordingly.
The predict optimization algorithm is used as part of a four-step modeling process: (i) coarse modeling. predict searches through the entire protein conformation space, evenly covered by ≈1,500 decoys, to identify regions of stability; (ii) fine modeling. predict is used a second time to comb the neighborhood of the most stable “coarse†models, optimizing ≈5,000 decoy structures in the vicinity of each model. This step allows the algorithm to rapidly focus on regions of stability in the protein's energy landscape and to efficiently identify the most stable “fine†model; (iii) molecular dynamics refinement. The resulting all-atom model is minimized and then subjected to up to 300-ps molecular dynamics simulations with charmm (29) and the charmm22 force field (30). Multiple constraints are applied during the simulations to ensure that the model does not deviate significantly from the predict model. These refinement dynamics introduce helical kinks and relax the side-chain conformations; and (iv) virtual protein–ligand complex. A protein–ligand complex is carefully constructed through molecular dynamics, mimicking the experimental cocrystallization process, which locks the target in a ligand-bound conformation.
