What's the suggested best practice to simulate peptide solubility? I'm considering to use Gromacs Tutorial Lysozyme in Water as the model.
Later, I'm considering to use Radius of Gyration as measure of solubility?
Is it a reasonable approach? Or is there a better way?
I'm not sure why you would use radius for gyration as a measure of solubility. Radius of gyration roughly measures how "unfolded" (or more accurately, how "stretched out") a polymer is. Although this may correspond to a higher solvent accessible surface area (which can itself be calculated and is useful in some cases), it has nothing to do with solubility.
If you want solubility, you'll likely want the free energy of solvation, or, as it is called when water is the solvent, the free energy of hydration. Alchemical free energy methods are one family of tools that have been useful for calculating these quantities.
This paper (full text here) compares experimental and simulated hydration free energies for several hundred small molecules. It describes the protocol used, and files to reproduce the simulations are available in the SI and on GitHub.
References within that paper, in particular in the section titled "Hydration and solvation free energies have a range of applications", may be useful to point you situations more analogous to the specific experimental comparison you have in mind.
The problem with solubility of large molecules, is that you need to know the chemical potential in the liquid and the solid. You can use Free Energy calculations to get the residual contribution to the chemical potential in the liquid, but, the residual contribution to the chemical potential of the solid is really hard. It can be done, but, not accurately. MD attempt at solubility
Macroscopic modelers do solubility calculations routinely by using an activity coefficient model to calculate the excess contribution to the chemical potential. here the excess is the total contribution to the chemical potential, minus the reference state (pure liquid) - similar, but not quite the same as a residual contribution. They also run into problems with the solid chemical potential, so they fall back on experiment.
I will save the derivation, but in the end what you need is the activity coefficient of of the solute in the solution and the experimental (or predicted) enthalpy of fusion of the solute.
Thus you solve:
$$\ln (x_i) = \frac{\Delta H_m}{R} \left(\frac{1}{T_m} - \frac{1}{T} \right) + RT \ln \gamma_i(x_i)$$
(note I need to verify this, my memory is hazy and Wikipedia is bad at this stuff. I will look up my own code later, it works, so... I trust it more than Wikipedia. I may have plus and/or minus signs wrong here).
Here we have the experimental enthalpy of melting (fusion) and some models activity coefficient(macroscopic equation of state, free energy calculation, divine intervention ...). note that the activity coefficient depends on the composition, so the composition (solubility) is on the left and right hand side. Thus, we have to solve this iteratively, which is why it is nice for macroscopic models, which takes fractions of seconds, and, not so nice for molecular models which take hours(have not actually done this yet using MD, it is on my list of things to do). I find that a simple secant method is fine when doing it using macroscopic equations of state. Simply updating will unfortunately lead to divergence sometimes.
It is possible to calculate the activity coefficient of the solute using molecular simulation (molecular dynamics or Monte Carlo) free energy of solvation calculations if you understand how to relate that output value to the activity coefficient, but is far from trivial... it is much easier to do wrong than right.
If using molecular simulations, you must understand that you are calculating a free energy with a different reference state than the typical activity coefficient models use (they are based on the Lewis-Randall chemical potential reference state (pure liquid)). COSMO models have to deal with the derivation for solubility using a chemical potential model that is not Lewis-Randall, so check that out to see how to use your models excess or residual chemical potential properly, or if using, for instance, a Henry law model, derive it carefully from scratch, and have many colleagues look it over.
Another option is to set up slabs of the solid in the solution, and then simply run MD or MC. the solids will slowly dissolve into the solution, and at some appropriate spot in the bulk, you measure the concentration of the solute. This is also not so trivial as it seems, but look to literature for guidance. A good question is "what is the solid (crystal) form of my drug molecule?". The answer is frequently unknown... Also, does the FF model the crystal solid well? probably not.
Try to correlate something to experimental results and use this fit to predict solubility. This type of rational design strategy makes me a bit queasy, but, sometimes science doesn't have a quantitative strategy that is computationally feasible, and we have to fall back on parameter fitting. Who knows, maybe there is a good correlation between radius of gyration and free energy, or solubility. You could throw in hydrogen bonding sites into the regression as well, the list of features never ends...
