BACKGROUND
solvation effects play an important role in protein folding and stability.
likewise, the processes of protein:protein and protein:ligand binding are accompanied by effects such as desolvation and rearrangement of solvent molecules.
these solvation effects can be calculated by explicit solvent models, such as molecular dynamics simulations using a sphere of water molecules  <cit> .
for certain large-scale applications, however, this explicit solvent treatment is too time-consuming.
in protein design, an enormous number of amino acid sequences need to be considered  <cit> .
likewise, the screening of a large library of ligand molecules as a function of their binding affinity to a protein is a costly procedure and requires efficient methods.
to avoid these problems, implicit solvent models are often used, which yield significant computational efficiency  <cit> .
they do not consider the solvent degrees of freedom explicitly, but treat the solvent as a continuous medium having the average properties of the real solvent.
empirical methods, such as the solvent accessible surface area  model  <cit> , often provide simple and quick ways of evaluating the solvation energy with an accuracy comparable to theoretical models.
asa models have become widely accepted within available implicit solvent treatments and have been used successfully in many applications, such as protein molecular dynamics  <cit> , structure prediction  <cit>  and protein:ligand binding  <cit> .
in the asa approach, the solvation free energy of a solute is expressed as a sum of atomic contributions, weighted by their solvent-exposed area.
the contribution of each atom is quantified by a surface coefficient, which reflects the hydrophobicity or hydrophilicity of the particular atom type.
apart from non-polar contributions to the solvation free energy, such as the entropy cost for cavity formation and van der waals interactions, electrostatic contributions play an important role.
due to the fitting of the asa model to experimental data, the electrostatic contribution is partly incorporated into the parameters.
however, especially when using a small number of atom types, it is necessary to additionally calculate a screening energy, which accounts for the shielding of protein-protein electrostatic interactions by the high dielectric solvent.
a simple approach is to add a term that reduces the electrostatic interactions between protein atoms by a constant factor, ε, which plays the role of a dielectric constant.
to balance the two components of the solvent model, an overall weight α is applied to the asa term.
this combined model is known as the coulomb/accessible surface area  model  <cit> .
more accurate approaches for screening energy calculations are the generalized born  model  <cit>  and the poisson-boltzmann  model  <cit> .
a disadvantage of these methods, however, is that they are not easily pairwise decomposable  <cit> : the energy is not expressed as a sum over atom pairs.
further, they are relatively time-consuming compared with surface area based models.
the first set of atomic solvation parameters distinguishing between different atom types was developed by eisenberg and mcclachlan in  <dig>  <cit> .
they used octanol to water transfer energies for  <dig> amino acids to derive solvation coefficients for  <dig> atom types.
subsequently, a number of studies have been devoted to the parameterization of the atomic surface coefficients.
they differ in the assignment of atoms to characteristic groups and in the experimental data that were used to fit the coefficients.
ooi et al.
<cit>  used  <dig> different atomic coefficients, fitted to experimental free energies of solvation of small organic molecules.
fraternali and van gunsteren  <cit>  restricted the atom types to only two: one for carbon, representing the hydrophobic effect, and one for both nitrogen and oxygen, representing the hydrophilic effect.
these two parameters were optimized, such that the hydrophobic and hydrophilic solvent-accessible surface areas in proteins obtained from md simulations matched those measured from the corresponding x-ray structures.
later, several more highly-parameterized asa models  <cit>  were developed, which use up to  <dig> different atom types and large training sets of experimental solvation free energies for diverse organic molecules.
in previous work  <cit> , we optimized the casa model for side chain placement and mutagenesis.
we modified the atomic solvation parameters of fraternali  <cit>  by including an additional surface coefficient for atoms in charged groups, and by optimizing the dielectric constant ε and the weight α of the surface term.
the model was optimized and tested using sidechain reconstruction calculations, protein solvation energies, and stability changes for point mutations involving the insertion or removal of a charged sidechain.
in this paper, we continue to explore the performance of the casa model, pursuing two directions.
first, we consider whether a specific treatment of aromatic groups can lead to an improved parameterization.
second, we consider a broader set of test calculations than before.
we again consider stability changes, but we include a wider variety of mutations, most of which do not affect charged groups.
we consider the calculation of protein:ligand binding free energy changes due to point mutations, a very important application.
so far, most ligand binding studies using asa models used a large set of atomic surface coefficients.
pei et al  <cit> , for example, used  <dig> atom types and reproduced the binding free energies of a test set of  <dig> protein:ligand complexes with a standard error of  <dig>  kcal/mol.
our own model is much less heavily parameterized, but yields a comparable accuracy.
finally, we perform automated protein design for eight small proteins of the sh <dig> family.
computational protein design is another area that requires good implicit solvent models  <cit> .
this approach can help engineer new proteins as well as predict protein structure.
it considers a given backbone structure of a protein and predicts the amino acid sequences that fold into it  <cit> .
this information can be used either to identify mutations that stabilize a given structure or to assign the given 3d structure to new sequences with yet undetermined protein structures.
the protein design procedure is applied here to a small test set of eight sh <dig> proteins and the performance of the optimized solvation parameters is compared to that of the earlier parameter set.
the newer casa parameterization derived here treats aromatic groups as a specific group.
all the model parameters are then reoptimized from scratch.
most studies so far have used experimental transfer free energies from octanol or cyclohexane to water for small model molecules to derive solvation parameters.
however, organic solvents are only a crude approximation of the protein interior  <cit> .
a more recent study by zhou et al.
<cit>  used a database of protein mutation experiments to develop atomic solvation parameters.
with these parameters, the binding free energies of  <dig> protein-protein complexes were predicted with an rms deviation of  <dig>  kcal/mol.
likewise, lomize et al.
<cit>  achieved very good agreement with experimental data using atomic solvation parameters based on protein stabilities.
here, we use a similar approach and derive our newer solvation parameters from experimental protein and peptide stability changes.
we employ a procedure that attempts to match the computed stability changes to the experimental values, using a set of  <dig> mutations.
simultaneously, the model for the unfolded reference state of a protein is optimized: for each amino acid type, a preferred unfolded conformation is chosen from a large library of tripeptide fragments obtained from various proteins.
the final parameter set yields improved performance for protein stabilities, as expected.
the newer and the earlier parameters yield comparable, good performance for ligand binding and protein design.
given the importance of implicit solvent models in the fields of structure-based drug design, prediction of protein structure and protein design, both of our optimized casa models should be valuable tools for a wide range of applications.
RESULTS
casa parameter optimization and stability calculations
in recent work, we optimized and tested the casa model, using just three atomic categories and performing sidechain reconstruction and protein stability calculations.
we refer to the corresponding set of atomic surface coefficients as the "modified fraternali", or mf parameters .
indeed, the starting point for the earlier optimization was the parameter set of fraternali & van gunsteren  <cit> .
here, we want to test the casa model further, and to explore whether a more extensive parameterization will improve its performance.
our previous casa model  <cit>  included only polar, unpolar and ionic atom coefficients.
the current protocol also distinguishes aromatic atoms as a separate group.
this refined atom assignment may be expected to improve the solvent model, as it has been shown that aromatic groups have solvation properties different from aliphatic groups  <cit> .
with this modification, we reoptimized all the model parameters from scratch.
we scanned a range of values for the four different surface coefficients; additionally, three different values for the dielectric constant ε in the coulomb screening term  were tested.
a data set of  <dig> experimental stability mutations was used to adjust the parameters: these data correspond to  <dig> different helical peptides and  <dig> different proteins.
the protein mutations involve mainly ionized residues, while the peptide mutations  cover the full range of amino acids.
in contrast to an earlier study by lomize et al.
<cit> , which derived atomic solvation parameters based on stability data for the protein interior, and in contrast to our previous study  <cit> , this set of mutations involves mainly solvent-exposed residues .
mf: modified fraternali parameters; phia: parameters optimized here.
to compute the stability of a protein, it is necessary to construct a model for the unfolded state.
here, we assume that sidechains do not interact with each other in the unfolded state, and we describe the environment of each sidechain with a simple tripeptide model.
each sidechain thus interacts only with the local tripeptide backbone and the solvent.
for each amino acid type, one preferred structure was chosen from a large set of tripeptide conformations obtained from various proteins.
since the choice of the preferred unfolded reference structure for each amino acid type also depends on the parameterization of the solvent model, the parameter optimization had to proceed iteratively :  a preferred tripeptide structure was chosen for each amino acid using the current set of solvation parameters ; the solvation parameters were adjusted to best reproduce the experimental stability changes using the current set of reference structures.
this was done based on the rms deviation of the computed stability changes from the experimental values.
starting from the previously optimized solvation parameters  <cit> , the optimization procedure converged after  <dig> iteration cycles.
table  <dig> shows the set of surface coefficients chosen from the top-ranking results.
we refer to the new set as the phia set .
compared with the earlier, mf parameters, the phia set shows a clear difference between the aromatic and the unpolar coefficient, a more positive ionic coefficient and a slightly more negative unpolar coefficient.
the experimental and computed stability changes are compared in figure  <dig>  the correlation between the two sets is modest, 43%.
nevertheless, the mean error is reasonable, with an rms deviation of  <dig>  kcal/mol for the complete set of experimental values and a mean unsigned error of  <dig>  kcal/mol .
the rank ordering of the data is characterized by a spearman rank correlation of 29%  <cit> ; the probability of obtaining this value by chance is less than  <dig>  .
the phia results represent a significant improvement over the performance of the mf parameters, which give an rms deviation of  <dig>  kcal/mol and a mean unsigned error of  <dig>  kcal/mol.
the dielectric constant for the new, phia set of surface coefficients was set to  <dig>  a higher dielectric constant of  <dig> gives only insignificant improvement , while a lower value of  <dig> leads to noticeably poorer agreement with experiment .
in addition to the criterium of minimal deviation from the experimental data, care was taken to select a parameter set that makes sense physically.
the ordering of the coefficients should follow the expected preference of solvent exposure of each atom group: unpolar < aromatic < polar < ionized.
further, the coefficients of different groups should be sufficiently separated from each other.
the surface coefficient of the unpolar atom group was allowed to be slightly negative, as this led to better agreement with the experimental data.
to ensure that this does not lead to unphysical behaviour, the surface coefficient was tested on a propane dimer system.
studies on methane and neopropane association in water report an association free energy of - <dig>  and - <dig>  kcal/mol, respectively  <cit> .
using values between  <dig>  and - <dig>  kcal/mol/Å <dig> for the nonpolar surface coefficient, the computed association free energy varies from - <dig>  to + <dig>  kcal/mol.
the nonpolar coefficient of - <dig>  kcal/mol/Å <dig> used here gives an acceptable association free energy of - <dig>  kcal/mol.
to verify that the parameters were not overly biased by the optimization procedure, we performed cross-validation tests .
part of the data  were omitted from the optimization process, then used to test the error level.
this led to very similar parameters and errors compared to the original parameter optimization.
specifically, one cross-validation run led to the following optimized parameters : aromatic: - <dig>  ; ionic: - <dig>  ; polar: - <dig>  ; nonpolar: - <dig>  .
the optimal dielectric was  <dig>  as before.
the mean and rms errors for the omitted data were  <dig>  and  <dig>  kcal/mol, respectively, compared to  <dig>  and  <dig>  for the optimization data.
the mean and rms errors with the phia parameters for the omitted data were similar:  <dig>  and  <dig>  kcal/mol.
the other cross-validation run led to the exact same atomic coefficients and dielectric constant, and to mean and rms errors for the omitted data of  <dig>  and  <dig>  kcal/mol .
thus the optimized parameters and the error levels are similar with and without cross-validation, showing that the resulting parameters are fairly robust.
the model performance shows a moderate dependence on the system considered.
the peptide models taken alone lead to a slightly better agreement with the experimental data , which might be due to their simple helical structure compared with a large protein system .
the performance of the uncharged mutations is better than that of the charged ones .
this occurs partly because the proportion of charged mutations is higher for the proteins .
compared with the results of lomize et al.
<cit> , which gave an rms deviation from experimental stability changes of only  <dig>  kcal/mol, the performance here is less good.
in their study, however, the solvation parameters were fitted to a more restricted set of experimental data.
only buried, uncharged residues in α-helices and β-sheets of proteins were used.
in contrast, in our study, both charged and uncharged, solvent-exposed and buried, and protein and peptide mutations were included.
therefore, it can be expected that our more varied set of mutations results in a higher deviation from experiment.
binding affinities
the performance of both the earlier, mf and the new, phia solvation parameters was assessed on a large set of experimental binding affinities and resistance mutations.
these data include  <dig> mutations in six different ligand-protein systems:  tyrosyl-trna synthetase in complex with tyrosine ,  aspartyl-trna synthetase in complex with aspartate ,  lysozyme in complex with the antibody hyhel- <dig>  and  the complex of the transmembrane glycoprotein cd <dig> with the gp <dig> component of the hiv virus ;  the bpti:trypsin complex;  the chymotrypsin:bpti complex.
a seventh system was also studied: the tyrosine kinase abl in complex with the drug imatinib.
for this system, information is only available on mutations leading to resistance to the drug imatinib  <cit> , while for the other four systems, more precise values of the binding affinities are available.
the energy function was slightly different from the one used for the stability mutations above.
here, the dielectric constant  ε was set to  <dig> and the weight α of the surface area term was set to  <dig>  .
these values gave improved performance for the binding affinities.
for the first four systems, , the calculation of binding affinities followed a very simple protocol that only optimizes the conformation of the side chain in the mutated position, while the rotamers of the remaining side chains are left unchanged apart from a slight minimization.
further, the starting structures of the ligand-bound state and the ligand-free state are assumed to be identical except for the side chain at the mutated position.
although this is a rather simple approach, the resulting computed binding affinities are in reasonable agreement with the experimental values.
the mean unsigned error for the differences in binding affinities upon mutation is  <dig>  kcal/mol with the phia parameters, and  <dig>  kcal/mol with the mf parameters .
thus, the performance of the two casa variants for this application is similar.
the phia results for all mutations are given in additional file  <dig>  the performance is slightly poorer when charged residues are mutated.
the more complicated systems such as protein-antibody binding  or protein-protein binding  give slightly higher deviations from experiment .
the mean unsigned errors for the small molecule-ligands , combined, is  <dig>  kcal/mol with phia, while the protein-ligand systems , combined, give a mean unsigned error of  <dig>  kcal/mol .
aexcluding the bpti complexes, which were computed with a slightly different protocol and for which gb data are not available.
bfor gb-ace and gb-hct, some mutations were excluded due to unfavorable vdw contacts; the numbers of mutations  are in parentheses.
cbpti:trypsin and :chymotrypsin complexes.
results were also computed for two other systems, the bpti:trypsin  <cit>  and bpti:chymotrypsin  <cit>  complexes.
<dig> mutations at position  <dig> in bpti were studied in each case.
the native residue is a lysine.
one mutation was excluded  due to a large van der waals contribution to the affinity change , leaving  <dig> mutations.
with the simple protocol used above, the agreement with experiment was poorer, with a mean error of  <dig>  kcal/mol for the  <dig> mutations.
a slightly different protocol was then tried.
the entire bpti was minimized, instead of just the mutated position .
this led to improved agreement, with mean errors of  <dig>  and  <dig>  kcal/mol for bpti:trypsin and bpti:chymotrypsin, respectively: about the same level as for the lysozyme:antibody complex.
this illustrates the need for more extensive structural relaxation with some mutations.
more generally, these two examples show that, not surprisingly, with the simple energy functions and conformational exploration used here, a certain amount of system-specific parameter fitting and adjustment can be necessary.
for the resistance mutations in abl-kinase, good qualitative agreement was achieved, in so far as all the mutated proteins gave less favorable binding affinities than the native protein .
for positions  <dig> and  <dig>  the computed difference between the mutant and native binding free energies may be too small; the other three values, however, are clearly consistent with the observed resistance to imatinib binding.
in kcal/mol.
there are several mutations that are badly predicted, including r59a in the cd4:gp <dig> system and d101g and k96m in the lysozyme:antibody system .
the first two cases involve a large to small mutation, and all three involve removal of a net charge; these processes which might require a more extensive rearrangement of the surrounding residues, as in the bpti complexes .
in the simple protocol employed here, the rotamers of the side chains in the vicinity of the mutation are not reoptimized, and the slight minimization carried out here might not be sufficient to model a realistic mutated protein conformation.
for the k96m case, we tried the same protocol as for the bpti cases, applying  <dig> steps of powell minimization to the entire lysozyme protein, instead of just the mutated sidechain.
this led to a computed binding free energy change of + <dig>  kcal/mol, much closer to the experimental value of  <dig>  kcal/mol.
for one other case, the d78a mutation in tyrrs, we employed a much more extensive conformational search: all the rotamers close to the mutation site were explored using a stochastic search strategy; the error in the binding free energy decreased from  <dig>  kcal/mol to  <dig>  kcal/mol .
these examples suggest that when the sidechain charge and/or volume changes substantially, a more sophisticated conformational sampling is required.
more work in this direction is underway.
on the whole, however, we obtain fair agreement with experimental data using a very simple method.
the overall mean unsigned error for the  <dig> mutations is  <dig>  kcal/mol; the rms error is  <dig>  kcal/mol.
the correlation between the computed and measured data is 73%.
the rank ordering of the data is characterized by a spearman rank correlation of 70%  <cit> ; the probability of obtaining this value by chance is less than  <dig>  .
our errors are only slightly higher than those reported by pei et al., who obtained an rms error of  <dig>  kcal/mol for the binding affinities of  <dig> protein-ligand complexes.
likewise, the accuracy obtained here appears comparable to that achieved by guérois et al.
<cit>  for their charged mutations, although this subset was not reported separately.
recently, handel et al.
<cit>  predicted the stabilities of more than  <dig> mutants to within  <dig> kcal/mol.
partly, this improvement might be due to their more sophisticated, all-atom, force field, which increases accuracy but makes their calculations almost an order of magnitude more expensive, compared with the charmm19/casa level .
further, our data set contains a higher percentage of charged mutations, which makes predictions more demanding.
as an additional reference, the performance of the optimized casa models was compared to that of a gb/sa model, which should provide a more accurate treatment of the electrostatic contribution to the solvation free energy  <cit> .
only the  <dig> mutations for systems  were studied.
two variants of the gb model were tested:  gb-hct  <cit>  in combination with the amber, all-atom force field  <cit>  and  gb-ace  <cit>  with the charmm <dig> force field  <cit> .
for the gb-hct variant, the overall quality of the results is close to the casa level, with a slightly higher mean error .
the gb-ace variant gives notably higher mean errors for all protein-ligand systems.
this might be due to the gb-ace parameterization which is not optimized for this specific application.
the gb-hct parameters were optimized previously for computational sidechain placement and protein mutagenesis  <cit> .
some differences between gb-hct and casa might be due to the different force field treatments, as amber uses explicit hydrogens on all atoms, while charmm <dig> uses implicit hydrogens for unpolar atoms.
the qualitative agreement with the experimental resistance mutations for abl:imatinib for the two gb variants is comparable to that obtained using the optimized casa models .
protein design
in protein design, the amino acid side chains are mutated and sequences are selected to optimize the folding free energy, using a heuristic search algorithm.
in our previous protein design study  <cit> , we obtained good results for  <dig> different globular proteins using the mf solvation parameters.
here, we consider a subset of those proteins, consisting of  <dig> sh <dig> domains.
these proteins are used to assess the performance of the new, phia solvation parameters for protein design.
the protein design calculations were carried out using our proteins@home distributed computing platform, with the help of volunteers in several countries.
proteins@home is discribed in more detail elsewhere  <cit> .
for each protein,  <dig>  sequences were generated and amino acid identities relative to the native sequence were calculated.
here, we give  the full identities of the designed sequences, even though cys, gly, and pro were not allowed to mutate; see tables  <dig> and  <dig>  for further analysis, we selected two subsets of the computed sequences:  the  <dig> sequences with the highest identity scores relative to the native protein ;  the  <dig> sequences with the best folding free energies .
for these two subsets, the mean sequence identities relative to the corresponding native sequence are given in table  <dig>  we will not discuss the mean identities over all  <dig>  computed sequences, because the values lie very close to those obtained for the low energy sequences.
for the low energy subset, the average sequence identity obtained with the phia parameters is  <dig> %, slightly lower than the value obtained earlier with the mf parameters,  <dig> %  <cit> .
for the eight proteins, the identities obtained with phia range from  <dig> % to  <dig> %.
the relative performance of the mf and phia parameters depends on the protein: since hck and c-src are better predicted with phia, while 1gcqb, crk, abl, and especially csk are better predicted with mf.
considering the subset of high-scoring sequences, the phia parameters give results of roughly the same quality as mf: the mean identities for this subset are  <dig> % with phia, compared to  <dig> % with mf.
with both the phia and the mf parameters, the identity scores obtained for the  <dig> proteins lie within the range of published average identity scores for redesigned proteins  <cit> .
in a protein design study by jaramillo et al.
<cit>  sequence optimizations for  <dig> sh <dig> domains were performed and resulted in an average sequence identity of  <dig> %.
our energy-ranked phia sequences lie well above this score, with a sequence identity of  <dig> % averaged over all  <dig> proteins.
recently, saunders et al.
<cit>  used a refined protein design method and reported sequence identities as high as 37% for  <dig> globular proteins.
considering the full sequence identities of our best-scoring sequences, both our parameter sets give results that lie close to this value.
pokala and handel  <cit>  used an all-atom force field and a gb/sa solvent model to redesign  <dig> proteins and achieved somewhat higher sequence identities, between  <dig>  and  <dig> %.
this approach, however, includes a negative design criterion, which constrains the surface amino acid composition of the proteins to be native-like.
it also leads to an increase in computational effort of about two orders of magnitude, compared with the casa solvent model and a united-atom force field such as charmm <dig>
next, we addressed the question whether the sequences computed with the new solvation parameters ressemble naturally occuring sequences.
as a reference, we created a set of natural sequences in the following way.
for each of the  <dig> proteins, sequences from the swissprot database with a sequence identity of more than 60% relative to the respective native sequence were retrieved.
these sequences were then combined into a large set consisting of  <dig> sequences.
table  <dig> compares the average blosum <dig> scores obtained for the natural sequence set with those for the complete set of computed sequences and for our two sets of high-ranking sequences .
a weighting scheme was applied to the blosum scores to account for the variability of each amino acid within the natural sequences.
this scheme gives a greater weight to amino acids that are conserved within the natural set.
for a frequency of an amino acid of more than 80% at a certain position within the natural sequence set, a weight of  <dig> was assigned; for a frequency between  <dig> and 80%, a weight of  <dig>  was assigned; for a frequency of less than 50%, a weight of  <dig>  was assigned.
the blosum <dig> scores show that the computed sequences cover a considerable part of the natural sequence set.
figure  <dig> shows that the computed scores overlap with the range of natural sequences.
for most proteins, the weighted average blosum score lies below the correspinding score for the natural sequence set, but is within the range of the natural sequences.
in one case , the weighted average blosum score even exceeds the value of the natural set.
this indicates that the computed sequences do indeed show native-like characteristics and behave like distant homologues of the native sequences.
discussion
simple, efficient, solvent models are of great importance in protein modelling and structural bioinformatics.
here, we have continued to explore the performance of the casa solvent model, in two directions.
first, we considered a variant of increased complexity, where aromatic atoms are treated as a separate group.
our previous approach  <cit>  did not distinguish between unpolar and aromatic atoms.
indeed, the solvation properties of aromatic groups are rather different from other nonpolar groups found in proteins.
for this variant, we reparameterized the model completely, leading to the phia parameter set.
second, we applied the casa model to a wider set of applications than previously.
we considered protein stability changes associated with point mutations, including all amino acid types .
we also considered protein:ligand binding, an especially important application.
finally, we performed complete protein redesign with the new parameters.
for the new, phia parameterization, four different atom types were considered for the atomic surface coefficients: unpolar, aromatic, polar and ionized atoms.
the solvation parameters were fitted to a set of  <dig> experimental stability changes for protein and peptide mutations.
protein stabilities were calculated using an unfolded reference state modelled by a collection of tripeptide structures.
these reference structures were taken from a large library of structural fragments from six different proteins.
starting from the earlier, mf solvation parameters  <cit> , an iterative procedure was employed, which optimizes, at each cycle, both the solvation parameters and the model for the unfolded reference state of a protein.
atomic parameters were chosen that gave a minimal deviation from the experimental data and also represent the expected relative hydrophobicities of the atom groups.
the selected parameter set gives a mean unsigned error of  <dig>  kcal/mol for the  <dig> stability mutations, compared to  <dig>  kcal/mol with the mf parameters.
cross-validation tests gave similar parameter values and similar error levels.
for the binding calculations, over  <dig> experimental mutations in  <dig> different protein-ligand systems were used, including both small molecule ligands and protein-protein complexes.
the calculated differences in binding free energy are in reasonable agreement with the experimental data, with both casa variants.
the mean unsigned error is  <dig>  kcal/mol with the phia parameters and  <dig>  kcal/mol with the mf parameters.
it was also shown that the optimized casa models are not inferior to methods such as gb-ace/sa or gb-hct/sa, which treat the electrostatic contribution to solvation more accurately.
two additional protein-protein complexes  required a slightly modified protocol to give comparable error levels.
protein design was carried out for eight sh <dig> domain proteins and the performance of the phia parameters was compared with that of the mf parameters.
on average, slightly lower sequence identities with the native sequence were achieved using the new, phia parameters.
differences, however, are small, and the performance depends on the particular protein.
on the whole, both parameter sets give sequences of comparable quality, which are competitive with recently published results for designed proteins.
further, the computed sequences were found to have the character of naturally occuring, distant homologues of the native sequences.
CONCLUSIONS
overall, the casa model performs well for a wide variety of applications.
the phia parameters, specifically optimized here for protein stability, give a distinct improvement over the earlier, mf parameters for this application.
for ligand binding and protein design, the specific treatment of aromatic groups in the phia parameterization did not lead to an improvement in performance.
rather, the two parameter sets perform well for both applications, with the exact relative performance dependening on the particular system.
both variants provide an efficient tool for the computational engineering of ligands and proteins.

 solvation effects play an important role in protein folding and stability. likewise, the processes of protein:protein and protein:ligand binding are accompanied by effects such as desolvation and rearrangement of solvent molecules. these solvation effects can be calculated by explicit solvent models, such as molecular dynamics simulations using a sphere of water molecules. for certain large-scale applications, however, this explicit solvent treatment is too time-consuming. in protein design, an enormous number of amino acid sequences need to be considered. likewise, the screening of a large library of ligand molecules as a function of their binding affinity to a protein is a costly procedure and requires efficient methods. to avoid these problems, implicit solvent models are often used, which yield significant computational efficiency. they do not consider the solvent degrees of freedom explicitly, but treat the solvent as a continuous medium having the average properties of the real solvent. empirical methods, such as the solvent accessible surface area model, often provide simple and quick ways of evaluating the solvation energy with an accuracy comparable to theoretical models. asa models have become widely accepted within available implicit solvent treatments and have been used successfully in many applications, such as protein molecular dynamics, structure prediction and protein:ligand binding. in the asa approach, the solvation free energy of a solute is expressed as a sum of atomic contributions, weighted by their solvent-exposed area. the contribution of each atom is quantified by a surface coefficient, which reflects the hydrophobicity or hydrophilicity of the particular atom type. apart from non-polar contributions to the solvation free energy, such as the entropy cost for cavity formation and van der waals interactions, electrostatic contributions play an important role. due to the fitting of the asa model to experimental data, the electrostatic contribution is partly incorporated into the parameters. however, especially when using a small number of atom types, it is necessary to additionally calculate a screening energy, which accounts for the shielding of protein-protein electrostatic interactions by the high dielectric solvent. a simple approach is to add a term that reduces the electrostatic interactions between protein atoms by a constant factor, ε, which plays the role of a dielectric constant. to balance the two components of the solvent model, an overall weight α is applied to the asa term. this combined model is known as the coulomb/accessible surface area model. more accurate approaches for screening energy calculations are the generalized born model and the poisson-boltzmann model. a disadvantage of these methods, however, is that they are not easily pairwise decomposable: the energy is not expressed as a sum over atom pairs. further, they are relatively time-consuming compared with surface area based models. the first set of atomic solvation parameters distinguishing between different atom types was developed by eisenberg and mcclachlan in. they used octanol to water transfer energies for amino acids to derive solvation coefficients for atom types. subsequently, a number of studies have been devoted to the parameterization of the atomic surface coefficients. they differ in the assignment of atoms to characteristic groups and in the experimental data that were used to fit the coefficients. ooi et al. used different atomic coefficients, fitted to experimental free energies of solvation of small organic molecules. fraternali and van gunsteren restricted the atom types to only two: one for carbon, representing the hydrophobic effect, and one for both nitrogen and oxygen, representing the hydrophilic effect. these two parameters were optimized, such that the hydrophobic and hydrophilic solvent-accessible surface areas in proteins obtained from md simulations matched those measured from the corresponding x-ray structures. later, several more highly-parameterized asa models were developed, which use up to different atom types and large training sets of experimental solvation free energies for diverse organic molecules. in previous work, we optimized the casa model for side chain placement and mutagenesis. we modified the atomic solvation parameters of fraternali by including an additional surface coefficient for atoms in charged groups, and by optimizing the dielectric constant ε and the weight α of the surface term. the model was optimized and tested using sidechain reconstruction calculations, protein solvation energies, and stability changes for point mutations involving the insertion or removal of a charged sidechain. in this paper, we continue to explore the performance of the casa model, pursuing two directions. first, we consider whether a specific treatment of aromatic groups can lead to an improved parameterization. second, we consider a broader set of test calculations than before. we again consider stability changes, but we include a wider variety of mutations, most of which do not affect charged groups. we consider the calculation of protein:ligand binding free energy changes due to point mutations, a very important application. so far, most ligand binding studies using asa models used a large set of atomic surface coefficients. pei et al, for example, used atom types and reproduced the binding free energies of a test set of protein:ligand complexes with a standard error of kcal/mol. our own model is much less heavily parameterized, but yields a comparable accuracy. finally, we perform automated protein design for eight small proteins of the sh family. computational protein design is another area that requires good implicit solvent models. this approach can help engineer new proteins as well as predict protein structure. it considers a given backbone structure of a protein and predicts the amino acid sequences that fold into it. this information can be used either to identify mutations that stabilize a given structure or to assign the given 3d structure to new sequences with yet undetermined protein structures. the protein design procedure is applied here to a small test set of eight sh proteins and the performance of the optimized solvation parameters is compared to that of the earlier parameter set. the newer casa parameterization derived here treats aromatic groups as a specific group. all the model parameters are then reoptimized from scratch. most studies so far have used experimental transfer free energies from octanol or cyclohexane to water for small model molecules to derive solvation parameters. however, organic solvents are only a crude approximation of the protein interior. a more recent study by zhou et al. used a database of protein mutation experiments to develop atomic solvation parameters. with these parameters, the binding free energies of protein-protein complexes were predicted with an rms deviation of kcal/mol. likewise, lomize et al. achieved very good agreement with experimental data using atomic solvation parameters based on protein stabilities. here, we use a similar approach and derive our newer solvation parameters from experimental protein and peptide stability changes. we employ a procedure that attempts to match the computed stability changes to the experimental values, using a set of mutations. simultaneously, the model for the unfolded reference state of a protein is optimized: for each amino acid type, a preferred unfolded conformation is chosen from a large library of tripeptide fragments obtained from various proteins. the final parameter set yields improved performance for protein stabilities, as expected. the newer and the earlier parameters yield comparable, good performance for ligand binding and protein design. given the importance of implicit solvent models in the fields of structure-based drug design, prediction of protein structure and protein design, both of our optimized casa models should be valuable tools for a wide range of applications.  casa parameter optimization and stability calculations in recent work, we optimized and tested the casa model, using just three atomic categories and performing sidechain reconstruction and protein stability calculations. we refer to the corresponding set of atomic surface coefficients as the "modified fraternali", or mf parameters. indeed, the starting point for the earlier optimization was the parameter set of fraternali & van gunsteren. here, we want to test the casa model further, and to explore whether a more extensive parameterization will improve its performance. our previous casa model included only polar, unpolar and ionic atom coefficients. the current protocol also distinguishes aromatic atoms as a separate group. this refined atom assignment may be expected to improve the solvent model, as it has been shown that aromatic groups have solvation properties different from aliphatic groups. with this modification, we reoptimized all the model parameters from scratch. we scanned a range of values for the four different surface coefficients; additionally, three different values for the dielectric constant ε in the coulomb screening term were tested. a data set of experimental stability mutations was used to adjust the parameters: these data correspond to different helical peptides and different proteins. the protein mutations involve mainly ionized residues, while the peptide mutations cover the full range of amino acids. in contrast to an earlier study by lomize et al., which derived atomic solvation parameters based on stability data for the protein interior, and in contrast to our previous study, this set of mutations involves mainly solvent-exposed residues. mf: modified fraternali parameters; phia: parameters optimized here. to compute the stability of a protein, it is necessary to construct a model for the unfolded state. here, we assume that sidechains do not interact with each other in the unfolded state, and we describe the environment of each sidechain with a simple tripeptide model. each sidechain thus interacts only with the local tripeptide backbone and the solvent. for each amino acid type, one preferred structure was chosen from a large set of tripeptide conformations obtained from various proteins. since the choice of the preferred unfolded reference structure for each amino acid type also depends on the parameterization of the solvent model, the parameter optimization had to proceed iteratively: a preferred tripeptide structure was chosen for each amino acid using the current set of solvation parameters; the solvation parameters were adjusted to best reproduce the experimental stability changes using the current set of reference structures. this was done based on the rms deviation of the computed stability changes from the experimental values. starting from the previously optimized solvation parameters, the optimization procedure converged after iteration cycles. table shows the set of surface coefficients chosen from the top-ranking results. we refer to the new set as the phia set. compared with the earlier, mf parameters, the phia set shows a clear difference between the aromatic and the unpolar coefficient, a more positive ionic coefficient and a slightly more negative unpolar coefficient. the experimental and computed stability changes are compared in figure the correlation between the two sets is modest, 43%. nevertheless, the mean error is reasonable, with an rms deviation of kcal/mol for the complete set of experimental values and a mean unsigned error of kcal/mol. the rank ordering of the data is characterized by a spearman rank correlation of 29%; the probability of obtaining this value by chance is less than. the phia results represent a significant improvement over the performance of the mf parameters, which give an rms deviation of kcal/mol and a mean unsigned error of kcal/mol. the dielectric constant for the new, phia set of surface coefficients was set to a higher dielectric constant of gives only insignificant improvement, while a lower value of leads to noticeably poorer agreement with experiment. in addition to the criterium of minimal deviation from the experimental data, care was taken to select a parameter set that makes sense physically. the ordering of the coefficients should follow the expected preference of solvent exposure of each atom group: unpolar < aromatic < polar < ionized. further, the coefficients of different groups should be sufficiently separated from each other. the surface coefficient of the unpolar atom group was allowed to be slightly negative, as this led to better agreement with the experimental data. to ensure that this does not lead to unphysical behaviour, the surface coefficient was tested on a propane dimer system. studies on methane and neopropane association in water report an association free energy of - and - kcal/mol, respectively. using values between and - kcal/mol/Å for the nonpolar surface coefficient, the computed association free energy varies from - to + kcal/mol. the nonpolar coefficient of - kcal/mol/Å used here gives an acceptable association free energy of - kcal/mol. to verify that the parameters were not overly biased by the optimization procedure, we performed cross-validation tests. part of the data were omitted from the optimization process, then used to test the error level. this led to very similar parameters and errors compared to the original parameter optimization. specifically, one cross-validation run led to the following optimized parameters: aromatic: -; ionic: -; polar: -; nonpolar: -. the optimal dielectric was as before. the mean and rms errors for the omitted data were and kcal/mol, respectively, compared to and for the optimization data. the mean and rms errors with the phia parameters for the omitted data were similar: and kcal/mol. the other cross-validation run led to the exact same atomic coefficients and dielectric constant, and to mean and rms errors for the omitted data of and kcal/mol. thus the optimized parameters and the error levels are similar with and without cross-validation, showing that the resulting parameters are fairly robust. the model performance shows a moderate dependence on the system considered. the peptide models taken alone lead to a slightly better agreement with the experimental data, which might be due to their simple helical structure compared with a large protein system. the performance of the uncharged mutations is better than that of the charged ones. this occurs partly because the proportion of charged mutations is higher for the proteins. compared with the results of lomize et al., which gave an rms deviation from experimental stability changes of only kcal/mol, the performance here is less good. in their study, however, the solvation parameters were fitted to a more restricted set of experimental data. only buried, uncharged residues in α-helices and β-sheets of proteins were used. in contrast, in our study, both charged and uncharged, solvent-exposed and buried, and protein and peptide mutations were included. therefore, it can be expected that our more varied set of mutations results in a higher deviation from experiment. binding affinities the performance of both the earlier, mf and the new, phia solvation parameters was assessed on a large set of experimental binding affinities and resistance mutations. these data include mutations in six different ligand-protein systems: tyrosyl-trna synthetase in complex with tyrosine, aspartyl-trna synthetase in complex with aspartate, lysozyme in complex with the antibody hyhel- and the complex of the transmembrane glycoprotein cd with the gp component of the hiv virus; the bpti:trypsin complex; the chymotrypsin:bpti complex. a seventh system was also studied: the tyrosine kinase abl in complex with the drug imatinib. for this system, information is only available on mutations leading to resistance to the drug imatinib, while for the other four systems, more precise values of the binding affinities are available. the energy function was slightly different from the one used for the stability mutations above. here, the dielectric constant ε was set to and the weight α of the surface area term was set to. these values gave improved performance for the binding affinities. for the first four systems,, the calculation of binding affinities followed a very simple protocol that only optimizes the conformation of the side chain in the mutated position, while the rotamers of the remaining side chains are left unchanged apart from a slight minimization. further, the starting structures of the ligand-bound state and the ligand-free state are assumed to be identical except for the side chain at the mutated position. although this is a rather simple approach, the resulting computed binding affinities are in reasonable agreement with the experimental values. the mean unsigned error for the differences in binding affinities upon mutation is kcal/mol with the phia parameters, and kcal/mol with the mf parameters. thus, the performance of the two casa variants for this application is similar. the phia results for all mutations are given in additional file the performance is slightly poorer when charged residues are mutated. the more complicated systems such as protein-antibody binding or protein-protein binding give slightly higher deviations from experiment. the mean unsigned errors for the small molecule-ligands, combined, is kcal/mol with phia, while the protein-ligand systems, combined, give a mean unsigned error of kcal/mol. aexcluding the bpti complexes, which were computed with a slightly different protocol and for which gb data are not available. bfor gb-ace and gb-hct, some mutations were excluded due to unfavorable vdw contacts; the numbers of mutations are in parentheses. cbpti:trypsin and:chymotrypsin complexes. results were also computed for two other systems, the bpti:trypsin and bpti:chymotrypsin complexes. mutations at position in bpti were studied in each case. the native residue is a lysine. one mutation was excluded due to a large van der waals contribution to the affinity change, leaving mutations. with the simple protocol used above, the agreement with experiment was poorer, with a mean error of kcal/mol for the mutations. a slightly different protocol was then tried. the entire bpti was minimized, instead of just the mutated position. this led to improved agreement, with mean errors of and kcal/mol for bpti:trypsin and bpti:chymotrypsin, respectively: about the same level as for the lysozyme:antibody complex. this illustrates the need for more extensive structural relaxation with some mutations. more generally, these two examples show that, not surprisingly, with the simple energy functions and conformational exploration used here, a certain amount of system-specific parameter fitting and adjustment can be necessary. for the resistance mutations in abl-kinase, good qualitative agreement was achieved, in so far as all the mutated proteins gave less favorable binding affinities than the native protein. for positions and the computed difference between the mutant and native binding free energies may be too small; the other three values, however, are clearly consistent with the observed resistance to imatinib binding. in kcal/mol. there are several mutations that are badly predicted, including r59a in the cd4:gp system and d101g and k96m in the lysozyme:antibody system. the first two cases involve a large to small mutation, and all three involve removal of a net charge; these processes which might require a more extensive rearrangement of the surrounding residues, as in the bpti complexes. in the simple protocol employed here, the rotamers of the side chains in the vicinity of the mutation are not reoptimized, and the slight minimization carried out here might not be sufficient to model a realistic mutated protein conformation. for the k96m case, we tried the same protocol as for the bpti cases, applying steps of powell minimization to the entire lysozyme protein, instead of just the mutated sidechain. this led to a computed binding free energy change of + kcal/mol, much closer to the experimental value of kcal/mol. for one other case, the d78a mutation in tyrrs, we employed a much more extensive conformational search: all the rotamers close to the mutation site were explored using a stochastic search strategy; the error in the binding free energy decreased from kcal/mol to kcal/mol. these examples suggest that when the sidechain charge and/or volume changes substantially, a more sophisticated conformational sampling is required. more work in this direction is underway. on the whole, however, we obtain fair agreement with experimental data using a very simple method. the overall mean unsigned error for the mutations is kcal/mol; the rms error is kcal/mol. the correlation between the computed and measured data is 73%. the rank ordering of the data is characterized by a spearman rank correlation of 70%; the probability of obtaining this value by chance is less than. our errors are only slightly higher than those reported by pei et al., who obtained an rms error of kcal/mol for the binding affinities of protein-ligand complexes. likewise, the accuracy obtained here appears comparable to that achieved by guérois et al. for their charged mutations, although this subset was not reported separately. recently, handel et al. predicted the stabilities of more than mutants to within kcal/mol. partly, this improvement might be due to their more sophisticated, all-atom, force field, which increases accuracy but makes their calculations almost an order of magnitude more expensive, compared with the charmm19/casa level. further, our data set contains a higher percentage of charged mutations, which makes predictions more demanding. as an additional reference, the performance of the optimized casa models was compared to that of a gb/sa model, which should provide a more accurate treatment of the electrostatic contribution to the solvation free energy. only the mutations for systems were studied. two variants of the gb model were tested: gb-hct in combination with the amber, all-atom force field and gb-ace with the charmm force field. for the gb-hct variant, the overall quality of the results is close to the casa level, with a slightly higher mean error. the gb-ace variant gives notably higher mean errors for all protein-ligand systems. this might be due to the gb-ace parameterization which is not optimized for this specific application. the gb-hct parameters were optimized previously for computational sidechain placement and protein mutagenesis. some differences between gb-hct and casa might be due to the different force field treatments, as amber uses explicit hydrogens on all atoms, while charmm uses implicit hydrogens for unpolar atoms. the qualitative agreement with the experimental resistance mutations for abl:imatinib for the two gb variants is comparable to that obtained using the optimized casa models. protein design in protein design, the amino acid side chains are mutated and sequences are selected to optimize the folding free energy, using a heuristic search algorithm. in our previous protein design study, we obtained good results for different globular proteins using the mf solvation parameters. here, we consider a subset of those proteins, consisting of sh domains. these proteins are used to assess the performance of the new, phia solvation parameters for protein design. the protein design calculations were carried out using our proteins@home distributed computing platform, with the help of volunteers in several countries. proteins@home is discribed in more detail elsewhere. for each protein, sequences were generated and amino acid identities relative to the native sequence were calculated. here, we give the full identities of the designed sequences, even though cys, gly, and pro were not allowed to mutate; see tables and for further analysis, we selected two subsets of the computed sequences: the sequences with the highest identity scores relative to the native protein; the sequences with the best folding free energies. for these two subsets, the mean sequence identities relative to the corresponding native sequence are given in table we will not discuss the mean identities over all computed sequences, because the values lie very close to those obtained for the low energy sequences. for the low energy subset, the average sequence identity obtained with the phia parameters is %, slightly lower than the value obtained earlier with the mf parameters, %. for the eight proteins, the identities obtained with phia range from % to %. the relative performance of the mf and phia parameters depends on the protein: since hck and c-src are better predicted with phia, while 1gcqb, crk, abl, and especially csk are better predicted with mf. considering the subset of high-scoring sequences, the phia parameters give results of roughly the same quality as mf: the mean identities for this subset are % with phia, compared to % with mf. with both the phia and the mf parameters, the identity scores obtained for the proteins lie within the range of published average identity scores for redesigned proteins. in a protein design study by jaramillo et al. sequence optimizations for sh domains were performed and resulted in an average sequence identity of %. our energy-ranked phia sequences lie well above this score, with a sequence identity of % averaged over all proteins. recently, saunders et al. used a refined protein design method and reported sequence identities as high as 37% for globular proteins. considering the full sequence identities of our best-scoring sequences, both our parameter sets give results that lie close to this value. pokala and handel used an all-atom force field and a gb/sa solvent model to redesign proteins and achieved somewhat higher sequence identities, between and %. this approach, however, includes a negative design criterion, which constrains the surface amino acid composition of the proteins to be native-like. it also leads to an increase in computational effort of about two orders of magnitude, compared with the casa solvent model and a united-atom force field such as charmm next, we addressed the question whether the sequences computed with the new solvation parameters ressemble naturally occuring sequences. as a reference, we created a set of natural sequences in the following way. for each of the proteins, sequences from the swissprot database with a sequence identity of more than 60% relative to the respective native sequence were retrieved. these sequences were then combined into a large set consisting of sequences. table compares the average blosum scores obtained for the natural sequence set with those for the complete set of computed sequences and for our two sets of high-ranking sequences. a weighting scheme was applied to the blosum scores to account for the variability of each amino acid within the natural sequences. this scheme gives a greater weight to amino acids that are conserved within the natural set. for a frequency of an amino acid of more than 80% at a certain position within the natural sequence set, a weight of was assigned; for a frequency between and 80%, a weight of was assigned; for a frequency of less than 50%, a weight of was assigned. the blosum scores show that the computed sequences cover a considerable part of the natural sequence set. figure shows that the computed scores overlap with the range of natural sequences. for most proteins, the weighted average blosum score lies below the correspinding score for the natural sequence set, but is within the range of the natural sequences. in one case, the weighted average blosum score even exceeds the value of the natural set. this indicates that the computed sequences do indeed show native-like characteristics and behave like distant homologues of the native sequences. discussion simple, efficient, solvent models are of great importance in protein modelling and structural bioinformatics. here, we have continued to explore the performance of the casa solvent model, in two directions. first, we considered a variant of increased complexity, where aromatic atoms are treated as a separate group. our previous approach did not distinguish between unpolar and aromatic atoms. indeed, the solvation properties of aromatic groups are rather different from other nonpolar groups found in proteins. for this variant, we reparameterized the model completely, leading to the phia parameter set. second, we applied the casa model to a wider set of applications than previously. we considered protein stability changes associated with point mutations, including all amino acid types. we also considered protein:ligand binding, an especially important application. finally, we performed complete protein redesign with the new parameters. for the new, phia parameterization, four different atom types were considered for the atomic surface coefficients: unpolar, aromatic, polar and ionized atoms. the solvation parameters were fitted to a set of experimental stability changes for protein and peptide mutations. protein stabilities were calculated using an unfolded reference state modelled by a collection of tripeptide structures. these reference structures were taken from a large library of structural fragments from six different proteins. starting from the earlier, mf solvation parameters, an iterative procedure was employed, which optimizes, at each cycle, both the solvation parameters and the model for the unfolded reference state of a protein. atomic parameters were chosen that gave a minimal deviation from the experimental data and also represent the expected relative hydrophobicities of the atom groups. the selected parameter set gives a mean unsigned error of kcal/mol for the stability mutations, compared to kcal/mol with the mf parameters. cross-validation tests gave similar parameter values and similar error levels. for the binding calculations, over experimental mutations in different protein-ligand systems were used, including both small molecule ligands and protein-protein complexes. the calculated differences in binding free energy are in reasonable agreement with the experimental data, with both casa variants. the mean unsigned error is kcal/mol with the phia parameters and kcal/mol with the mf parameters. it was also shown that the optimized casa models are not inferior to methods such as gb-ace/sa or gb-hct/sa, which treat the electrostatic contribution to solvation more accurately. two additional protein-protein complexes required a slightly modified protocol to give comparable error levels. protein design was carried out for eight sh domain proteins and the performance of the phia parameters was compared with that of the mf parameters. on average, slightly lower sequence identities with the native sequence were achieved using the new, phia parameters. differences, however, are small, and the performance depends on the particular protein. on the whole, both parameter sets give sequences of comparable quality, which are competitive with recently published results for designed proteins. further, the computed sequences were found to have the character of naturally occuring, distant homologues of the native sequences.  overall, the casa model performs well for a wide variety of applications. the phia parameters, specifically optimized here for protein stability, give a distinct improvement over the earlier, mf parameters for this application. for ligand binding and protein design, the specific treatment of aromatic groups in the phia parameterization did not lead to an improvement in performance. rather, the two parameter sets perform well for both applications, with the exact relative performance dependening on the particular system. both variants provide an efficient tool for the computational engineering of ligands and proteins. 
this protocol uses a set of experimental stability mutations and a large set of tripeptide conformations to find the best tripeptide structures and solvation parameters. both variants provide an efficient solvent treatment for the computational engineering of ligands and proteins. we test a set of model parameters optimized earlier, and we also carry out a new optimization in this work, using as a target a set of experimental stability changes for single point mutations of various proteins and peptides. for ligand binding, we obtain reasonable agreement with a set of experimental mutation data, with a mean unsigned error of kcal/mol with the new parameters and kcal/mol with the earlier ones. likewise, the new parameters perform well for the design of sh domain proteins where an average of % sequence identity relative to the native sequences was achieved. the optimization procedure is general, and could be used with other force fields. the computation of stability changes requires a model for the unfolded state of the protein. protein structure prediction and computational protein design require efficient yet sufficiently accurate descriptions of aqueous solvent. the performance of the casa model is assessed by two further applications: calculation of protein-ligand binding affinities and computational protein design. in our approach, this state is represented by tripeptide structures of the sequence ala-x-ala for each amino acid type x. we followed an iterative optimization scheme which, at each cycle, optimizes the solvation parameters and a set of tripeptide structures for the unfolded state. overall, the two casa variants explored here perform very well for a wide variety of applications. using the optimized parameters, we obtain a mean unsigned error of kcal/mol for the stability mutations. we show that the optimized casa model is not inferior to the generalized born/surface area model for the prediction of these binding affinities. further, it was shown that the computed sequences have the character of naturally-occuring homologues of the native sequences.
 protein structure prediction and computational protein design require efficient yet sufficiently accurate descriptions of aqueous solvent. we continue to evaluate the performance of the coulomb/accessible surface area implicit solvent model, in combination with the charmm molecular mechanics force field. we test a set of model parameters optimized earlier, and we also carry out a new optimization in this work, using as a target a set of experimental stability changes for single point mutations of various proteins and peptides. the optimization procedure is general, and could be used with other force fields. the computation of stability changes requires a model for the unfolded state of the protein. in our approach, this state is represented by tripeptide structures of the sequence ala-x-ala for each amino acid type x. we followed an iterative optimization scheme which, at each cycle, optimizes the solvation parameters and a set of tripeptide structures for the unfolded state. this protocol uses a set of experimental stability mutations and a large set of tripeptide conformations to find the best tripeptide structures and solvation parameters.  using the optimized parameters, we obtain a mean unsigned error of kcal/mol for the stability mutations. the performance of the casa model is assessed by two further applications: calculation of protein-ligand binding affinities and computational protein design. for these two applications, the previous parameters and the ones optimized here give a similar performance. for ligand binding, we obtain reasonable agreement with a set of experimental mutation data, with a mean unsigned error of kcal/mol with the new parameters and kcal/mol with the earlier ones. we show that the optimized casa model is not inferior to the generalized born/surface area model for the prediction of these binding affinities. likewise, the new parameters perform well for the design of sh domain proteins where an average of % sequence identity relative to the native sequences was achieved. further, it was shown that the computed sequences have the character of naturally-occuring homologues of the native sequences.  overall, the two casa variants explored here perform very well for a wide variety of applications. both variants provide an efficient solvent treatment for the computational engineering of ligands and proteins. 