CHARMM c30b1 sasa.doc



File: SASA ]-[ Node: Top
Up: (commands.doc) -=- Next: Syntax\n


                   The SASA implicit solvation model


   Please read the section 'Syntax of the SASA command' before using
SASA.


                   Characteristics of the SASA model

   The SASA model is a fast implicit solvation model that is useful to
simulate structured peptides and miniprotein motifs [1]. It is based on
the solvent accessible surface area and uses only two solvation
parameters. It approximates the solvent accessible surface area with a
simple analytical function that is easily derivable, therefore the
speed. The model accounts for the electrostatic screening between solute
charges by using a distant dependent dielectric function and by
neutralizing the formal charges (Asp, Glu, Arg, Lys, and the termini) as
in the EEF1 model, see [2].

   The SASA model has been successfully applied to peptides, removing
the major artifacts of in vacuo simulations and even obtaining
reversible folding [3]. Latest benchmarks indicate that a simulation
with SASA is only about 50% slower than an in vacuo simulation.

                          Theoretical aspects

   We assume - as is usually done - that the potential energy of the
system consisting of the solute and the solvent can be decomposed in
three parts: the intra-solute potential energy U(X), the intra-solvent
potential energy V(Y) and the interaction potential energy of the solute
and the solvent W(X,Y), where X denotes the degrees of freedom of the
solute and Y the degrees of freedom of the solvent. Integrating out all
solvent degrees of freedom one gets the potential of mean force W(X),
also called the effective energy. It can be written as the sum of the
intra-solute potential energy U(X) and the mean solvation term DW(X)
that describes the solvent induced effects. This is a rigorous result
from statistical mechanics. For more details please refer for instance
to the review [4].

   Any implicit solvation model based on the solvent accessible surface
area loosely relies on the idea that the major contribution to the free
energy of solvation of the solute is determined by the surface of the
solute that is accessible to the solvent. Following this idea, it is
assumed that the mean solvation term is a sum of atomic contributions
proportional to the solvent accessible surface area of each atom. The
constants of proportionality are the solvation parameters. This
approximation includes the free energy cost for the transfer of an atom
from the interior of a protein to the solvent. However, it doesn't
account for the screening of solute charges.

                  Technical aspects of the SASA model

   There are three important aspects any implicit solvation model that
is based on the solvent accessible surface area has to take care of:
(A) how to calculate the solvent accessible surface area, (B) how to
account for the screening of solute charges, and (C) how many solvation
parameters to use.

         (A) Calculation of the solvent accessible surface area

   In the SASA model the solvent accessible surface area is
approximated by a probabilistic approach. For details please refer to
[5]. The base is a simple formula for the probability to hit the
accessible surface of atom i if N atoms are present by choosing
randomly a spot on the solvation shell of atom i. However, this formula
is only true under the assumption that all the atoms are distributed
randomly. This is of course not the case, due to the fact alone that
the van der Waals spheres must not penetrate each other. Now instead of
elaborating in sophisticated probability calculations the pragmatic
approach of [6] is adapted where the original formula is parametrized by
including two sets of parameters: a set of atom type parameters and a
set of connectivity parameters. The atom type parameters depend on the
atom type and correct for systematic errors primarily due to
hybridization. The connectivity parameters distinguish bound atoms from
more distant ones. Also, the radii used to calculate the approximated
solvent accessible surface were optimized for this purpose. (These radii
are used only for the calculation of the approximated area and they are
different from the radii used for the CHARMM van der Waals energy term.)
The optimal values for these three sets of parameters were taken from
[6]. For proteins subsets of the sets given in [6] are sufficient. The
SASA model uses the same subsets that are taken in [7], where 14
different atom types are used.

   The only internal degree of freedom considered by this approximation
is the interatomic distance. Therefore the major defect is the absence
of a better correlation between exact and approximated areas upon
changes in internal coordinates like dihedral angles.

                    (B) Screening of solute charges

   For the screening of solute charges the SASA model uses a distant
dependent screening function, eps(r)=2r, and neutralizes the charged
groups of polar amino acids in exactly the same way as it is done in the
EEF1 model proposed by Lazaridis and Karplus [2].

                        (C) Solvation parameters

   By default, the SASA model uses only two non-vanishing solvation
parameters: one for hydrophobic and one for hydrophilic groups. The
solvation of explicit hydrogen atoms is neglected. However, this can
be changed by the user (see below). The solvation parameters were
adjusted to the EEF1-modified CHARMM 19 polar hydrogen parameter set
by Philippe Ferrara in a trial and error approach. The criterion was to
minimize the root mean square deviation from the native state for six
small proteins by performing molecular dynamics simulations of 1ns at
300K [1]. For this parametrization process the SASA standard setup was
used (apart from the solvation parameters that were of course varied).
See below for the SASA standard setup.

                        Implementation in CHARMM

   To every CHARMM atom type the SASA model assigns a solvation
parameter (zero by default for explicit hydrogen atoms and non-zero for
hydrophobic and hydrophilic groups), a radius optimized for the
approximation of the solvent accessible surface area and an atom type
parameter. Additionally, a connectivity parameter is assigned to every
pair depending on whether the two atoms of the pair form a 1-2 or a more
distant pair. These parameters are called the SASA parameters.

   When initializing SASA, i.e., by invoking the SASA command, all the
SASA parameters are printed to the CHARMM output file. A value of
-999.000 means that this specific parameter hasn't yet been determined
for SASA and therefore the corresponding CHARMM atom type can't be used
by default for the SASA calculations - the user would have to assign a
meaningful value. (Currently, SASA does not support the following CHARMM
atom types by default: HA, HT, LP, CT, CM, NP, OH2, OM, OT, OS, and FE.)
All SASA parameters can be changed by the user. For the corresponding
syntax, see below.

   To evaluate the approximative formula for the solvent accessible
surface area, all pairs within the cutoff are required. Most, but not
all of these pairs, are included in the nonbond pair list in CHARMM. The
missing pairs are stored in a new pair list, the SASA pair list. More
precisely the SASA pair list contains all pairs that are excluded from
the nonbond pair list and that don't belong to the fixed exclusions (as
given in the topology file) in the case of a positive nonbond exclusion
mode and that are within the cutoff. This list assures correct operation
of the SASA model for any nonbond exclusion mode (-5 to 5). This means
that for any nonbond exclusion mode from 1 to 5, the SASA energy and its
derivatives are identical. The same is true for any nonbond exclusion
mode from -5 and 0. Of course there is a difference between the nonbond
exclusion modes 0 and 1 corresponding to the fact that the fixed
exclusions are treated differently: For an exclusion mode from -5 to 0
they are included in the nonbond pair list, opposed to an exclusion mode
from 1 to 5 where they are excluded from the nonbond pair list.

                         Range and limitations

   The SASA model has been applied to structured peptides, see for
instance [3]. However, you currently can't expect it to work for large
proteins mainly for two reasons. Firstly, it has been parametrized for
small proteins [1] and secondly, the dielectric function does not take
different environments into account, i.e., it does not distinguish
whether or not the interacting partial charges are buried or on the 
protein surface.

                                 Caveat

   The most distant pairs that can appear in the SASA pair list are 1-4
pairs. Since the 1-2, 1-3, and 1-4 pairs are usually very close in space,
there is no check for the pairs in the SASA pair list whether or not
these pairs are within the nonbonding cutoff. Hence, if one runs a
simulation with a cutoff so short that 1-4 pairs are outside the cutoff,
the SASA mean solvation energy and its derivatives are no longer
consistent with the cutoff. See also below for the SASA standard setup.

   Please note that the radii used for the calculation of the
approximated solvent accessible surface area are labelled 'van der Waals
radii' in the SASA output of the CHARMM output file. This is consistent
with the terms used in [6]. However, these radii are different from and
do not replace in any way the CHARMM default van der Waals radii used to
calculate the van der Waals energy term.

   Two additional files are needed to use SASA (taken from EEF1):

(1) toph19_eef1.inp : This is a modification of toph19.inp where ionic
		      side chains and termini are neutralized and contains
		      an extra parameter type (CR).
(2) param19_eef1.inp: This is a modification of param19.inp which includes
                      the extra parameter type (CR).

These files can be found in test/data/.

* Menu:

* Syntax::         Syntax of the SASA command
* References::     Useful references    
* Example::        Input file 




File: SASA ]-[ Node: Syntax
Up: Top -=- Next: References -=- Previous: Top\n


                       Syntax of the SASA command


   There is only one SASA command:


SASA atom-selection [S<atom-type> <real>] [R<atom-type> <real>] [P<atom-type> <real>] [fco <real>] [nco <real>]


atom-selection:== (see *note select:(select.doc).)

<atom-type>:== a CHARMM atom type from param19_eef1.inp, i.e., one of
               the following: H HC HA HT
                              LP
                              CT C CH1E CH2E CH3E CR1E CM
                              N NR NP NH1 NH2 NH3 NC2
                              O OC OH1 OH2 OM OT OS
                              S SH1E
                              FE

   This command sets up the SASA model for a simulation. All the values
have to be given on one command line. Invoking the SASA command a second
time reinitializes all values either to the default or to the user
specified values.

atom-selection        This determines the atoms to be used for the SASA
                      calculations. Any atom not selected here is
                      treated by SASA as if not existent: Not only is
                      the solvation free energy not calculated but also
                      the decrease in the solvent accessible surface
                      area of the selected atoms due to the not selected
                      atoms is not considered.

S<atom-type> <real>   Changes the solvation parameter of the CHARMM atom
                      type <atom-type> from the default value to <real>.

R<atom-type> <real>   Changes the radius used by SASA (for the
                      calculation of the approximated solvent accessible
                      surface area) of the CHARMM atom type <atom-type>
                      from the default value to <real>.

P<atom-type> <real>   Changes the atom type parameter of the CHARMM atom
                      type <atom-type> from the default value to <real>.

fco <real>            Changes the connectivity parameter for 1-2 pairs
                      from the default value to <real>.

nco <real>            Changes the connectivity parameter for more
                      distant than 1-2 pairs from the default value to
                      <real>.

   Please also have a look at the following examples. The SASA standard
setup looks like this:


nbond nbxmod 5 atom rdiel shift vatom vdistance vshift -
      cutnb 8.0 ctofnb 7.5 ctonnb 6.5 eps 2.0 e14fac 0.4 wmin 1.5

sasa selection (.not. hydrogen) end


   This command sets up SASA in its standard form. If you would like to
change the solvation parameters, for instance, use the following syntax:


nbond nbxmod 5 atom rdiel shift vatom vdistance vshift -
      cutnb 8.0 ctofnb 7.5 ctonnb 6.5 eps 2.0 e14fac 0.4 wmin 1.5

sasa selection (.not. hydrogen) end -
     SN -0.07  SNR  -0.07  SNH1 -0.07  SNH2 -0.07  SNH3 -0.07  SNC2 -0.07  SO   -0.07  SOC -0.07 SOH1 -0.07 -
     SC  0.014 SCH1E 0.014 SCH2E 0.014 SCH3E 0.014 SCR1E 0.014 SS    0.014 SSH1E 0.014


   This command sets up SASA and changes the solvation parameters of the
hydrophilic groups from the default of -0.06 to -0.07 and of the
hydrophobic groups from the default of 0.012 to 0.014. Note that you
have to specify every group separately. (Do not use these values. They
were chosen arbitrarily just for the sake of illustration of the
syntax.)

   The nonbond options are the default nonbond options from the default
param19.inp file with the exception of the eps value that is set to 2
instead of 1 and rdiel is used instead of cdiel. The default solvation
parameters are -0.06 for hydrophilic groups (N, NR, NH1, NH2, NH3, NC2,
O, OC, and OH1) and 0.012 for hydrophobic groups (C, CH1E, CH2E, CH3E,
CR1E, S, and SH1E) and 0.0 for explicit hydrogen atoms. Please note that
the param19_eef1.inp file has different default nonbond options
(especially the cutoffs are different) that were not used to parametrize
SASA and therefore should not be used or used only with care in
simulations with SASA if consistency is important!

   A second set of parameters has been tested: -0.144 for hydrophilic
groups, 0.024 for hydrophobic groups and 0.0 for explicit hydrogen atoms
with eps(r)=r. It is not the default. It seems to work better for large
proteins but no results have been published up to date (December 2001).

   If you select all atoms for SASA instead of excluding the (explicit)
hydrogens, all atoms will be considered for calculating the solvent
accessible surface area of each atom. However, since the solvation
parameter for hydrogen atoms is zero by default, the solvation energy of
the hydrogen atoms is also zero by default. If you insist on including
the solvation energy due to the solute hydrogen atoms, you have to
assign a non-zero solvation parameter to the hydrogen atoms by yourself,
using 'SH <real> SHC <real>' in the CHARMM input file. Be aware that
this is not the default.

   If you run a simulation with several molecules (so that you have more
than one segment identifier), make sure that you invoke the SASA command
after the generation of the last segment since any molecule generated
after the last use of the SASA command is not included in the SASA
calculations. You don't have to invoke the SASA command after every
generation of a segment, it is sufficient to use the SASA command once
after the generation of the last segment.

   The SASA solvation energy is stored in the variable 'SASL'. Use
'?SASL' in the CHARMM input file to access the value.



File: SASA ]-[ Node: References
Up: Top -=- Next: Example -=- Previous: Syntax\n


                               References

[1] Ferrara, P.; Apostolakis, J.; Caflisch, A.; Evaluation of a Fast
    Implicit Solvent Model for Molecular Dynamics Simulations; Proteins
    2002; 46; 24-33.

[2] Lazaridis, T.; Karplus. M.; Effective Energy Function for Proteins
    in Solution; Proteins 1999; 288; 477-487.

[3] Ferrara, P.; Caflisch, A.; Folding Simulations of a Three-Stranded
    Antiparallel Beta-Sheet Peptide; Proc. Natl. Acad. Sci. USA; 2000;
    97; 10780.

[4] Roux, B.; Simonson, T.; Implicit Solvent Models; Biophysical
    Chemistry; 78; 1999; 1-20.

[5] Wodak, S. J.; Janin, J.; Analytical Approximation to the Solvent
    Accessible Surface Area of Proteins; Proc. Natl. Acad. Sci. USA;
    1980; 77; 1736.

[6] Hasel, W.; Hendrickson, T. F.; Clark, S. W.; A Rapid Approximation
    to the Solvent Accessible Surface Area of Atoms; Tetrahedron
    Computer Methodology 1988; Vol. 1; No. 2; 103-116.

[7] Fraternali, F.; van Gunsteren, W. F.; An Efficient Mean Solvation
    Force Model for Use in Molecular Dynamics Simulations of Proteins
    in Aqueous Solution; J. Mol. Biol. 1996; 256; 939.




File: SASA ]-[ Node: Example
Up: Top -=- Next: Top -=- Previous: References\n


                                Examples


   Check test/c29test/sasa.inp for a more complex example with
minimization, equilibration and dynamics. Here is a short version:

* Example input file for the SASA implicit solvation model.
*



! --- Begin generation procedure ---

open read card name toph19_eef1.inp unit 30
read rtf card unit 30
close unit 30

open read card name param19_eef1.inp unit 30
read parameter card unit 30
close unit 30

open read card name filename.crd unit 30
read sequence coor unit 30
close unit 30

generate main warn setup

! --- End generation procedure ---



! --- Begin reading coordinates ---

open read card name filename.crd unit 30
read coordinate card unit 30 
close unit 30

! --- End reading coordinates ---



! --- Begin setting up SASA ---

! Use the SASA standard setup.

nbond nbxmod 5 atom rdiel shift vatom vdistance vshift -
      cutnb 8.0 ctofnb 7.5 ctonnb 6.5 eps 2.0 e14fac 0.4 wmin 1.5

sasa selection (.not. hydrogen) end

! --- End setting up SASA ---



! --- Begin minimization ---

minimize sd   nstep 300 nprint 20 tolgrad 0.1
minimize conj nstep 200 nprint 20 tolgrad 0.1

! --- End minimization ---



stop

CHARMM .doc Homepage


Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology
FDA/CBER/OVRR Biophysics Laboratory
Modified, updated and generalized by C.L. Brooks, III
The Scripps Research Institute