Use Steepest Descent to Build a Better Water Model for Computer-Aided Drug Design
A great part of my research is devoted to improving the accuracy of binding affinity prediction of drug candidates. You may think that the medicines we take are all made of substances with overwhelmingly complex chemical formulas. On the contrary, orally active drugs are often relatively small molecules (molecular mass <500 daltons) with rather simple structures (10 or fewer rotatable bonds), according to the famous Lipinski’s rule.
Below shows one of the most important drugs ever developed in human history: penicillin. More accurately, it is penicillin G, one of the many derivatives of penicillin (derivatives are a group of compounds with small structural tweaks in organic chemistry). The penicillin drugs have been extensively used as antibiotics, and have saved numeruous lives over the years.
Oftentimes, the very first step for such a small molecule to serve as a drug is to interact with certain types of proteins in the human body. Once the drug candidate interacts with proteins, it has a shot at triggering a cascade of physiological events aimed to cure diseases. If the molecule does not even interact with the proteins, most likely it won’t do much about the targeted disease. There is a term for this type of interactions between the protein and small drug-like compounds called “binding”, and we use another term called “binding affinity” to quantitatively measure how tight the protein and drug candidates interact with each other. An excellent introduction on protein-small molecule binding, drug design, and the measurement of binding affinities can be found here. Of course, to make a good drug, a strong binding is only the beginning of a long story: there are many other aspects a medicinal chemist needs to pay attention to, such as absorption, distribution, metabolism, and excretion (ADME). But like I mentioned, binding affinity measurement is a critical component in the early stage of drug discovery process.
One way of measuring the binding affinities of chemical compounds and identifying lead compounds (the starting point of a drug, promising but not good enough, plenty of room to improve) is through experiment. Millions of compounds can be screened withn a few days in big pharmaceutical companies, and this step is highly automated nowadays, mostly done by robotic arms.
Alternatively, binding affinities can be calculated through computational modeling and simulations. If good and fast enough predictions can be generated in silico, it can save a lot of cost and time of synthesizing the compounds and testing them using biochemical assays. In addition, computational work can provide a lot of details about how the drug-like molecule and protein interaction with each other at atomic level, which is usually lacking in experiment.
Most computational predictions can be boiled down to two aspects: accuracy and precision. Accuracy refers to the closeness of a measured value to a standard or known value. Precision refers to the closeness of two or more measurements to each other (definition). Your measurements can be accurate but not precise, precise but not accurate, neither, or both. For example, the water density at 1 atm is 1.000 kg/liter at room tempature. If someone measures the water density under the same temperature and pressure conditions and reports 0.807, 0.806 and 0.808 kg/liter, I would consider those data points as precise but not accurate. For comparison, 1.208, 1.008 and 0.814 are realtively more accurate because the averaged value is very close to our target value, but still not precise enough. Smaller variance among measurements are expected given that this is not such a hard experiment to do. In contrast, data sets including 1.005, 1.008 and 0.997 are more in the realm of “accurate and precise”.
In the real-world modeling and simulations, you can generate very rough predictions using oversimplified model in no time, but the outcomes are likely to be neither accurate nor precise. Or, you could use sophisticated model and increase the sampling size, but you could still end up with not totally satisfying results, due to the limitations of your model and parameters. And that’s sort of where we are right now on computing binding affinities using complex physics-based methods — the awkward space between the third and the fourth targets: good on precision but still sturuggling a lot on accuracy. The scatter plot shows the correlation between the computed and experimental binding affinties of a series of host-guest complexes. In brief, host-guest molecules are model systems we used to study protein-ligand binding.
We often use Linear regression model to fit the predicted binding affinities to the experimental values. In the scatter plot below, the solid black line is the “one-to-one” ratio line, which indicates a perfect binding prediction. The red dashed line is the trendline generated by the linear regression model for the values predicted using water model 1, and the green dashed line is the trendline for water model 2 predicted values. As you can see here, in both cases, there is a strong linear relationship between the predicted and measured values. However, the green spheres (affinities predicted using water model 2) are collectively closer to the one-to-one line than the red squares (affinities predicted using water model 1), which suggests a more accurate prediction.
So how did we go from water model 1 to water model 2? Basically we used gradients to guide our parameter search. We wanted to simulate an aqueous environments because drugs will be dissolved in our body fluid, and an aqueous environment can be approximated as thousands of or more water molecules which look like this:
Whether we like it or not, we end up with spending most of the time simulating the water, other than the stuff we are actually interested in, because the water molecules take up most of the space in our simulation box. But, water is important for molecular dynamics simulations, and it is really, really hard to get it right.
The simplest three-site (one oxgyen atom and two hydrogen atoms, no virtual sites) rigid water model has several parameters, including van der Waals parameters of the oxygen atom, the partial charges, the bond length and the angle parameters. To optimize the water parameters, we computed the first-order gradients of the binding affinities with respect to the parameters we wanted to refine, and the gradients can tell us how much and to what direction the computed binding affinity will change, once we make a small perturbation of the water parameters. We then used those information to reduce the RMSE value (the root-mean-square deviation). Another round of simulations were carried out, using the perturbed parameters, and gradients were computed once again from new simulation data. Those steps can be done iteratively until the error could not be reduced further by changing parameters. And the backbone of this optimization procedure is gradient descent, i.e. the steepest descent.