%% LyX 2.1.3 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english,noae,round]{article}
\usepackage{mathpazo}
\usepackage[T1]{fontenc}
\usepackage[letterpaper]{geometry}
\geometry{verbose,tmargin=2.5cm,bmargin=2.5cm,lmargin=2.5cm,rmargin=2.5cm}
\usepackage{color}
\usepackage{babel}
\usepackage{amsbsy}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage[authoryear]{natbib}
\usepackage[unicode=true,
 bookmarks=true,bookmarksnumbered=false,bookmarksopen=false,
 breaklinks=false,pdfborder={0 0 1},backref=false,colorlinks=true]
 {hyperref}
\hypersetup{pdftitle={An introduction to CHNOSZ},
 pdfauthor={Jeffrey M. Dick},
 citecolor=blue}
\usepackage{breakurl}

\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
<<echo=F>>=
  if(exists(".orig.enc")) options(encoding = .orig.enc)
@
\newenvironment{lyxcode}
{\par\begin{list}{}{
\setlength{\rightmargin}{\leftmargin}
\setlength{\listparindent}{0pt}% needed for AMS classes
\raggedright
\setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt}
\normalfont\ttfamily}%
 \item[]}
{\end{list}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
%\VignetteIndexEntry{An introduction to CHNOSZ}

% so DOIs in bibliography show up as hyperlinks
\newcommand*{\doi}[1]{\href{http://dx.doi.org/#1}{doi: #1}}

\makeatother

\begin{document}

\title{An introduction to CHNOSZ}


\author{Jeffrey M. Dick}

\maketitle
<<echo=FALSE>>=
options(width=90)
@


\section{About}

This document will orient you to the basic functionality of CHNOSZ,
a package for the R software environment. R is a powerful language
and also very fun to use. Don't worry if you're new to it; just plow
through the examples below and you'll start to get the hang of it.
If you want a more structured approach to learning the language, there
are some excellent guides in the Manuals section of the \href{http://www.r-project.org}{R Project page}.

The package was developed since 2006 to support a research project
on the thermodynamic properties of proteins. Since that time, the
functions in the package have expanded to include calculation of the
thermodynamic properties of reactions, and especially the construction
of equilibrium chemical activity diagrams for both inorganic and organic
systems. The development of the package since 2009 has focused on
the calculation of the equilibrium chemical activities of large numbers
of proteins with applications to interpretation of metagenomic data
and protein abundances in a variety of settings.

The database and functions are flexible in their use, allowing one
to model the relative stabilities of proteins, minerals or aqueous
species using very similar commands. Examples below are intended to
demonstrate basic usage to new users.


\subsection{Outline of workflow}

CHNOSZ is made up of a set of functions and supporting datasets. The
major components of the package are shown in the figure below, which
is an updated version of the flowchart from \citet{Dic08} (boxes
-- functions; ellipses -- datasets; bold text -- major user functions).

\includegraphics[bb=0bp 440bp 612bp 792bp,clip,width=0.75\columnwidth]{flowchart}

Some common usage scenarios are:
\begin{itemize}
\item using \texttt{info()} to search for species in the thermodynamic database;
\item using \texttt{subcrt()} to calculate the thermodynamic properties
of species and reactions;
\item using the sequence \texttt{basis()}, \texttt{species()}, \texttt{affinity()},
\texttt{equilibrate()}, \texttt{diagram()} to assign the basis species
that define the dimensions of chemical composition in a system, define
the species of interest for relative stability calculations, calculate
the affinities of formation reactions of the species of interest under
reference (non-equilibrium) conditions, calculate the equilibrium
chemical activities, and finally plot the results;\footnote{\texttt{equilibrate()} appeared in version 0.9-9 of CHNOSZ. In previous
versions, the equilibrium calculations were invoked by calls to \texttt{diagram()}.}
\item using \texttt{revisit()} to calculate/plot statistics of the chemical
activities of the species of interest and \texttt{findit()} to search
for combinations of activities of basis species, temperature and/or
pressure that optimize those statistics. These features, first appearing
in version 0.9-3 of the package, are covered briefly toward the end
of this document.
\end{itemize}
The functions are designed with an interactive setting in mind; you
can use CHNOSZ without having to write your own scripts. The examples
in this vignette are meant to portray a simple interactive session.
However, as you become more familiar with CHNOSZ and R, you will probably
find it helpful to save sequences of function calls that produce interesting
results. The results can then be reproduced on demand by yourself
or others with whom you might share your scripts.


\subsection{Installing and loading CHNOSZ}

If you have just installed R, and you are online, installing the CHNOSZ
package should be as simple as selecting ``Install packages from
CRAN'' or similar menu item in the R GUI or using the following command
to start the package installation process. (If you are not online,
you instead have to tell R to install the package from a local package
file.)

<<install_CHNOSZ,eval=FALSE>>=
install.packages("CHNOSZ")
@

Then load the CHNOSZ package to make its functions available in your
working session.

<<library_CHNOSZ>>=
library(CHNOSZ)
@

Then load the \texttt{thermo} object, which contains the thermodynamic
database and is also where your system settings will be stored.

<<data_thermo>>=
data(thermo)
@

The rest of this document assumes that the CHNOSZ package and data
are loaded.


\section{Thermodynamic database}


\subsection{Warning about internal consistency of thermodynamic data}

All thermodynamic data and examples are provided on an as-is basis.
It is up to you to check not only the accuracy of the data, but also
\emph{the suitability of the data AND computational techniques} for
your problem. By combining data taken from different sources, it is
possible to build an inconsistent and/or nonsensical calculation.
An attempt has been made to provide a primary database (OBIGT.csv)
that is internally consistent, but no guarantee can be made. Where
possible, data with known or suspected inconsistencies have been placed
into a secondary database (OBIGT-2.csv) that should be regarded as
experimental. If there is any doubt about the accuracy or suitability
of data for a particular problem, please \emph{consult the primary
sources} (which can be located using \texttt{browse.refs()}; see Section
\ref{sub:thermo_refs}). Do not assume that by adding any species
to your calculation (or to any of the examples), you will necessarily
obtain a reasonable answer. Do not assume that the examples are correct,
or that they can be applied to your problem. As with the data, please
\emph{compare the construction and output of the examples to the primary
sources}, cited in the reference list in each help page. Examples
without a reference (and some with references) demonstrate experimental
features of CHNOSZ. 


\subsection{\texttt{info()} part I}

So you want to know what are the standard molal thermodynamic properties
and equations of state parameters of aqueous ethylene? Look no further
than the \texttt{info()} function, which provides a convenient interface
to retrieve entries from the thermodynamic database packaged with
CHNOSZ.

<<info_ethylene>>=
info("ethylene")
@

There are two species named ``ethylene'' in the database. Normally,
\texttt{info()} gives preference to aqueous species if they exist,
so in this case we find that aqueous ethylene is species number 88
in the database. Let's display this entry, now by giving the species
index to the function.

<<info_88>>=
info(88)
@

If you were instead interested in the properties of the gas, you could
run:

<<info_ethylene_gas>>=
info("ethylene","gas")
@

\texttt{info()} itself is used by other functions in the package.
It prints output to the screen, but also returns a numeric value if
it finds a species matching the search term. So, we can retrieve the
properties of aqueous acetic acid without having to type in the species
ID number.

<<info_acetic>>=
aadata <- info(info("acetic acid"))
print(aadata)
@


\subsection{\texttt{\label{sub:thermo_refs}thermo\$refs}}

The thermodynamic data and other parameters used by the functions,
as well as system definitions provided by the user in an interactive
session, are stored in a list object called \texttt{thermo}.

<<summary_thermo>>=
summary(thermo)
@

Within this list, the thermodynamic database is contained in a data
frame (an R object that is like a matrix with named columns), \texttt{thermo\$obigt},
and the references to the original sources of thermodynamic data in
the literature are listed in \texttt{thermo\$refs}. Many of the authors
who are responsible for these data would be grateful if you cite them
whenever these data are used in publications! Use the \texttt{browse.refs()}
function without any arguments to show citation information for all
of the references in a browser window. You can include a species index
number to open the URL(s) associated with that entry in the database
(this requires an Internet connection).

<<browse.refs,eval=F>>=
browse.refs(88)
@

\texttt{browse.refs: opening URL for SH90 (E. L. Shock and H. C. Helgeson,
\citeyear{SH90})}


\subsection{\texttt{info()} part II}

Want to know what acids are in the database?

<<info_acid>>=
info("acid")
@

Here, \texttt{info()} couldn't find an exact match to a name, so it
performed a fuzzy search. That's why ``uracil'' and ``metacinnabar''
show up above. If you really just want species whose names include
the term ``acid'', you can add a placeholder character to narrow
the search. (Note: don't use an underscore (``\_'') here because
that character is reserved for names of proteins. Any other character
will do; here we use a space.)

<<info_spaceacid>>=
info(" acid")
@

The names of species other than proteins use (almost) exclusively
lowercase letters. \texttt{info()} can also be used to search the
text of the chemical formulas as they are entered in the database;
the symbols for the elements always start with a capital letter. The
example below lists the formulas of aqueous species, then minerals,
that contain the symbol commonly used to represent the hydroxide group.

<<info_OH>>=
info("(OH)")
@


\section{Proteins}


\subsection{\texttt{protein()}}

There are few things more fun than calculating the standard molal
Gibbs energy of formation from the elements at 25 $^{\circ}$C and
1 bar of a protein using group additivity. And there are few proteins
whose thermodynamic properties are more well studied than lysozyme
from the egg of the chicken.

<<protein_LYSC>>=
ip <- iprotein("LYSC_CHICK")
aa <- ip2aa(ip)
aa2eos(aa)
@

What happened there? Well, the first line found the row number (6)
of \texttt{thermo\$protein} that contains the amino acid composition
of LYSC\_CHICK. The second line extracted as a data frame. The third
line used amino acid group contributions \citep{DLH06} to calculate
the standard molal thermodynamic properties and equations of state
parameters of the aqueous protein species. There are other functions
available for calculating e.g. the chemical formula of the protein.

<<formula_LYSC>>=
pf <- protein.formula(aa)
as.chemical.formula(pf)
@


\subsection{\texttt{info()}}

Most of the time you probably won't be using the \texttt{iprotein()}
function. That's because \texttt{info()} recognizes the underscore
character as being an essential part of the name of a protein. The
names of proteins in CHNOSZ are mostly consistent with those used
in \href{http://www.uniprot.org}{Swiss-Prot/UniProtKB}.

<<info_LYSC>>=
si <- info("LYSC_CHICK")
info(si)
@

When CHNOSZ is first loaded, the thermodynamic properties and parameters
of the proteins are not present in \texttt{thermo\$obigt}. Therefore,
the first call to\texttt{ info()} just above had a side effect of
adding the computed properties and parameters to \texttt{thermo\$obigt}.


\section{Reaction properties}


\subsection{A single species}

A major feature of CHNOSZ is the ability to calculate standard molal
properties of species and reactions as a function of temperature and
pressure. The function used is called \texttt{subcrt()}, which takes
its name (with modification) from the well known SUPCRT package \citep{JOH92}.
\texttt{subcrt()}, like \texttt{info()}, has the name of a species
(including proteins) as its first argument (it also works if you give
it the numeric index of the species in the database). If no reaction
coefficients are given, the function calculates the standard molal
properties of the indicated species on a default temperature-pressure
grid.

<<subcrt_water>>=
subcrt("water")
@

The columns in the output are temperature ($^{\circ}$C), pressure
(bar), density of water (g cm$^{-3}$), logarithm of the equilibrium
constant (only meaningful for reactions; see below), and standard
molal Gibbs energy and enthalpy of formation from the elements (cal
mol$^{-1}$), and standard molal entropy (cal K$^{-1}$ mol$^{-1}$),
volume (cm$^{3}$ mol$^{-1}$) and heat capacity (cal K$^{-1}$ mol$^{-1}$). 

Compared to other species available in CHNOSZ, the equations for calculating
the properties of liquid $\mathrm{H_{2}O}$ are quite complex. The
package uses a Fortran subroutine taken from SUPCRT for these calculations.
See \texttt{help(water)} for more information.


\subsection{A reaction}

To calculate the properties of a reaction, enter the stoichiometric
reaction coefficients as a second argument to \texttt{subcrt()}. Reactants
have negative coefficients, and products have positive coefficients.
The function call below also shows the specification of temperature.

<<subcrt_C2H5OH>>=
subcrt(c("C2H5OH","O2","CO2","H2O"),c(-1,-3,2,3),T=37)
@

For historical reasons (i.e., the prevalence of the use of oxygen
fugacity in geochemical modeling; \citealp{And05}), $\mathrm{O_{2}}$
breaks the general rule in CHNOSZ that species whose states are not
specified are given the aqueous designation if it is available in
the thermodynamic database. If you want to specify the physical states
of the species in the reaction, that's possible too. For example,
we can ensure that dissolved $\mathrm{O_{2}}$ instead of the gaseous
form is used in the calculation. 

<<subcrt_C2H5OH_O2aq>>=
subcrt(c("C2H5OH","O2","CO2","H2O"),c(-1,-3,2,3),c("aq","aq","aq","liq"),T=37)
@

A useful feature of \texttt{subcrt()} is that it emits a warning if
the reaction is not balanced. Let's say you forgot to account for
oxygen on the left-hand side of the reaction\footnote{This example is motivated by the unbalanced reaction found at the
\href{http://en.wikipedia.org/wiki/Ethanol_metabolism}{Wikipedia entry on ethanol metabolism}
on 2010-09-23 and still present as of 2011-08-15: ``Complete Reaction:
C$_{2}$H$_{6}$O(Ethanol)$\rightarrow$C$_{2}$H$_{4}$O(Acetaldehyde)$\rightarrow$C$_{2}$H$_{4}$O$_{2}$(acetic
Acid) $\rightarrow$Acetyl-CoA$\rightarrow$3H$_{2}$O+2CO$_{2}$''.}.

<<subcrt_C2H5OH_unbal>>=
subcrt(c("C2H5OH","CO2","H2O"),c(-1,2,3),T=37)
@

The function still reports the results of the calculations, but use
them very cautiously (only if you have a specific reason for writing
an unbalanced reaction). In the next section we'll see how to use
another feature of CHNOSZ to automatically balance reactions.


\section{Basis species}


\subsection{What are basis species?}

\textbf{Basis species} are a minimal number of chemical species that
represent the compositional variation in a system. Operationally,
a \textbf{system} is the combination of basis species and species
of interest which is set up by the user to investigate a real-life
system. The basis species are akin to thermodynamic components, but
can include charged species. 

There are at least two reasons to define the basis species when using
CHNOSZ. First, you might want to use them to automatically balance
reactions. Second, they are required for making chemical activity
diagrams. Let's start with an example that \emph{doesn't} work.

<<basis_not,eval=FALSE>>=
basis(c("CO2","H2O","NH3","H2S","H+"))
@
\begin{lyxcode}
Error~in~put.basis(basis,~mystates)~:~

~~the~stoichiometric~matrix~must~be~square~and~invertible

In~addition:~Warning~messages:

1:~basis:~5~compounds~(~CO2~H2O~NH3~H2S~H+~)~

2:~basis:~6~elements~(~C~H~N~O~S~Z~)~
\end{lyxcode}
A limitation of CHNOSZ is that the number of basis species must be
equal to the number of elements, plus one if charge is present. This
way, any possible species of interest made up of these elements can
be compositionally represented by a linear combination of the basis
species. Now let's write a working basis definition.

<<basis>>=
basis(c("CO2","H2O","NH3","O2","H2S","H+"))
@

First basis definition! Note the column names, which give CHNOSZ its
name. These represent the elements in the commonly-occurring amino
acids, together with charge, denoted by ``Z''.


\subsection{Auto-balancing a reaction}

Now that the basis species are defined, try the unbalanced reaction
again.

<<subcrt_C2H5OH_auto>>=
subcrt(c("C2H5OH","CO2","H2O"),c(-1,2,3),T=37)
@

Here, \texttt{subcrt()} detected an unbalanced reaction, but since
the missing element was among the elements of the basis species, it
added the appropriate amount of $\mathrm{O_{2}}_{\left(gas\right)}$
to the reaction before running the calculations. You can go even further
and eliminate $\mathrm{CO_{2}}$ and $\mathrm{H_{2}O}$ from the function
call, but still get the same results.

<<subcrt_C2H5OH>>=
subcrt(c("C2H5OH"),c(-1),T=37)
@

What if you were interested in the thermodynamic properties of the
reaction of ethanol to acetaldehyde, but didn't want to balance the
reaction yourself (and you also didn't know how the formulas of the
species are written in the database)?

<<subcrt_ethanol>>=
subcrt(c("ethanol","acetaldehyde"),c(-1,1),T=37)
@

Notice how 2 H's needed to be added to the right-hand side of the
reaction; in our definition of basis species this comes out to $\mathrm{H_{2}O}-0.5\mathrm{O_{2}}$.
With a different choice of basis species, but the same elements, the
reaction might look quite different. As an example, suppose you had
amino acids in mind. The first line below, \texttt{data(thermo)},
is a quick way to reset the \texttt{thermo} object to its original
state, in order to clear the current system definition.

<<subcrt_acetaldehyde>>=
data(thermo)
basis(c("glutamic acid","methionine","isoleucine","lysine","tyrosine","H+"))
subcrt(c("ethanol","acetaldehyde"),c(-1,1),T=37)
@

In this case, the function finds that 2 H's are the compositional
equivalent of $0.5\mathrm{C_{6}H_{13}NO_{2}}-0.125\mathrm{C_{6}H_{14}N_{2}O_{2}}-0.250\mathrm{C_{9}H_{11}NO_{3}}$.
It's pretty easy for the computer to figure that out using matrix
operations, but probably isn't something you'd want to do by hand.
You might complain that this reaction is not likely to represent an
actual metabolic process ... as always, the challenge (and fun) of
coming up with a useful basis definition is in relating the species
to observable quantities.


\subsection{It works for proteins too!}

Let's set the basis definition again, this time using a keyword that
refers to a preset combination of basis species commonly encountered
in the documentation for CHNOSZ. Then we will use \texttt{subcrt()}
to calculate the thermodynamic properties of a reaction to form a
protein from the basis species.

<<subcrt_LYSC>>=
data(thermo)
basis("CHNOS+")
subcrt("LYSC_CHICK",1,T=25)
@

Note that using the keyword argument in \texttt{basis()} also set
the logarithms of activities (or fugacity in the case of $\mathrm{O_{2}}_{\left(g\right)}$)
to nominal values. While these settings do not affect the results
of the \texttt{subcrt()} calculation (which normally returns only
the standard molal properties of the reaction), they are essential
for calculating the relative stabilities of the species of interest.

If the protein is not available in CHNOSZ's own database, the amino
acid composition of the protein can be retrieved from the UniProtKB
(if the computer is connected to the internet). N.B. The following
example is not evaluated when compiling this vignette in case the
R package checks are run without internet access.

<<subcrt_ALAT1,eval=FALSE>>=
aa <- uniprot.aa("ALAT1_HUMAN")
add.protein(aa)
subcrt("ALAT1_HUMAN",1,T=25)
@


\section{Activity diagrams}


\subsection{Carbonate speciation (Bjerrum diagram)}

The sequence of commands \texttt{basis}-\texttt{species-affinity}-\texttt{{[}equilibrate{]}-diagram},
with various arguments, can be used to create a wide variety of diagrams.
The first two lines below configure the basis species and the aqueous
species. Some of the elements in the basis species are not in the
species of interest but that's OK (the opposite wouldn't be). We want
to make an activity diagram as a function of a single variable, and
configure the program to do this by the single \texttt{pH} argument
to \texttt{affinity()}. The \texttt{ylim} in the first \texttt{diagram()}
is optional; it ``zooms in'' the $y$-axis (by default, the $y$-axis
is expanded to contain all the values on the lines). Also the last
two lines are optional, unless you really do want to see the effect
of temperature (I do!)

\setkeys{Gin}{width=0.6\textwidth}

<<Bjerrum_diagram,fig=TRUE,width=4,height=4>>=
basis("CHNOS+")
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH=c(4, 12))
e <- equilibrate(a)
diagram(e, ylim=c(-6, 0))
a <- affinity(pH=c(4, 12), T=150)
e <- equilibrate(a)
diagram(e, add=TRUE, col="red")
@

\setkeys{Gin}{width=1.0\textwidth}

This just shows the speciation (relative abundances) of the aqueous
carbonate species. Calculating the solubility of a carbonate mineral,
or of $\mathrm{CO_{2}}$ gas, is possible, but more involved; try
the following:

\setkeys{Gin}{width=0.6\textwidth}

<<demo_solubility, results=hide, fig=TRUE>>=
demo("solubility", ask=FALSE)
@

\setkeys{Gin}{width=1.0\textwidth}

See the code of the demo (look for \texttt{demo/solubility.R} in the
directory where CHNOSZ is installed) to change the calculation from
calcite to $\mathrm{CO_{2}}$.


\subsection{Stability diagram for proteins}

Suppose that we are asked to calculate the relative stabilities of
some proteins from different organisms. We will use part of a case
study from \citet{Dic08}. \emph{Methanocaldococcus jannaschii} is
a hyperthermophilic methanogen known to live at higher temperatures
than \emph{Methanococcus voltae} (also a methanogen) and \emph{Haloarcula
japonica }(a halophile). These archaeal organisms produce cell-surface
glycoproteins (a.k.a. surface-layer proteins).

After defining the basis species we can define the \textbf{species
of interest}, i.e. those proteins whose relative stabilities we wish
to calculate. 

<<CSG_species>>=
basis("CHNOS")
species(c("SLAP_ACEKI", "CSG_METJA", "CSG_METVO", "CSG_HALJP"))
@

Note the output: the matrix denotes the coefficients of each of the
basis species in the formation reaction for one mole of each of the
species of interest. The \textbf{formation reaction} is the chemical
reaction to form one mole of a species of interest (as a product)
from a combination of basis species (as reactants and/or products,
depending on the stoichiometric constraints). The formation reactions
generally are \emph{not} statements about the mechanisms of reactions.
The species definition also includes reference values for the chemical
activities of the species of interest. 

Now we are all set up to calculate the chemical affinities of the
formation reactions. The chemical affinity is the negative of the
Gibbs energy change of a reaction per unit of reaction progress; it
is calculated in CHNOSZ using $\boldsymbol{A}=2.303RT\log(K/Q)$ ($R$
-- gas constant, $T$ -- temperature, $K$ -- equilibrium constant,
$Q$ -- activity product).

\texttt{affinity()} can accept arguments describing the range of chemical
conditions we're interested in. The names of the arguments can refer
to the basis species. Here, we vary the logarithm of the fugacity
of oxygen. The chemical activities of the other basis species are
taken to be constants equal to the values shown above.

<<CSG_affinity>>=
a <- affinity(O2=c(-90,-70))
@

Now we can use \texttt{equilibrate()} to calculate the equilibrium
activities of the proteins and \texttt{diagram()} to plot them. \texttt{normalize=TRUE}
invokes the normalization of the chemical formulas of the proteins
by the lengths (number of residues), which in most cases is desirable
for calculating equilibrium activities of proteins. We'll also specify
where the legend should be placed on the plot.

\setkeys{Gin}{width=0.6\textwidth}

<<CSG_diagram,fig=T,width=4,height=4>>=
e <- equilibrate(a, normalize=TRUE)
diagram(e, legend.x="bottomleft", ylim=c(-6, -2))
@

\setkeys{Gin}{width=1.0\textwidth}

Notably, the protein from the organism found at the highest temperatures
is relatively stable at more reduced conditions.


\subsection{More proteins, more dimensions}

Now let's add some bacterial surface-layer proteins. They are in some
way functional analogs (but not \href{http://en.wikipedia.org/wiki/Homology_(biology)}{homologs})
of the archaeal cell-surface glycoproteins. This plot is made using
the ``maximum affinity method'' to calculate the predominant species;
alternatively we could insert a call to \texttt{equilibrate()} in
order to calculate the equilibrium activities of all the species,
but the resulting predominance diagram would be identical.

\setkeys{Gin}{width=0.6\textwidth}

<<PredominanceDiagram,fig=T,width=4,height=4>>=
species(c("SLAP_ACEKI", "SLAP_GEOSE", "SLAP_BACLI", "SLAP_AERSA"))
basis(c("NH3", "H2S"), c(-1, -10))
a <- affinity(O2=c(-85, -70), H2O=c(-5, 0))
diagram(a, normalize=TRUE)
@

\setkeys{Gin}{width=1.0\textwidth}

Equilibrium predominances for proteins as a function of two chemical
activities! If you don't like the colors in the plot, don't worry...
the colors can be changed by using the \texttt{col} argument of \texttt{diagram()}.
This example hints at the multidimensional nature of the stability
problem. Note how the order of predominance fields at $\log a_{\mathrm{H_{2}O}}=0$
matches the order of proteins with highest equilibrium activities
in the previous diagram. Interpreting the meaning of low activities
of $\mathrm{H_{2}O}$ in these calculations remains a challenge.

Why did we increase the activity of $\mathrm{NH_{3}}$ and decrease
that of $\mathrm{H_{2}S}$? It was done here in order to increase
the size of the equilibrium predominance fields of the bacterial proteins.
This behavior is a result of the elemental makeup of the proteins:
the bacterial proteins under consideration are, per residue, more
nitrogen-rich and sulfur-poor than their archaeal counterparts (except
for CSG\_HALJP, which has no sulfur). CHNOSZ has a function to display
the compositional makeup of the proteins, per residue (indicated by
\texttt{normalize=TRUE}), in terms of the basis species.

<<residue.info>>=
protein.basis(species()$name, normalize=TRUE)
@


\subsection{A mineral example}

This example is modeled after a figure on p. 246 of \citet{BJH84}
for the system HCl-H2O-CaO-CO2-MgO-(SiO2) at 300 $^{\circ}$C and
1000 bar.

\setkeys{Gin}{width=0.6\textwidth}

<<Bowers,fig=T,width=4,height=4>>=
basis(c("HCl","H2O","Ca+2","CO2","Mg+2","SiO2","O2","H+"),
  c(999,0,999,999,999,999,999,-7))
species(c("quartz","talc","forsterite","tremolite","diopside",
  "wollastonite","monticellite","merwinite"))
a <- affinity("Mg+2"=c(-12,-4),"Ca+2"=c(-8,0),T=300,P=1000)
diagram(a)
@

\setkeys{Gin}{width=1.0\textwidth}

The 999's in the assignment of logarithms of activities of basis species
could be any number -- these settings do not affect the outcome of
the calculation. This is so because 1) $\mathrm{HCl}$, $\mathrm{CO_{2}}$
and $\mathrm{O_{2}}$ have zero stoichiometric coefficients in the
species, 2) the activities of $\mathrm{Ca}^{+2}$ and $\mathrm{Mg}^{+2}$
correspond to the axes of the diagram, and their ranges are taken
from the call to \texttt{affinity()}, and 3) $\mathrm{SiO_{2}}$ is
the immobile (conserved) component. Also note that ``Mg+2'' and
``Ca+2'' are not valid names of objects in R, but we can use them
as names of arguments by putting them in quotation marks in the call
to \texttt{affinity()}. 

Here, the scales of the axes here depend on the pH setting. This calculation
is therefore logically different from the formulation used by \citet{BJH84},
where the axes are $\log\left(a_{\mathrm{Mg^{+2}}}/\sigma_{\mathrm{Mg}^{+2}}a_{\mathrm{H^{+}}}^{2}\right)$
and $\log\left(a_{\mathrm{Ca^{+2}}}/\sigma_{\mathrm{Ca}^{+2}}a_{\mathrm{H^{+}}}^{2}\right)$
(where $\sigma$ is a function of the solvation of the subscripted
species). However, the geometry of the stability fields in the diagram
produced here is consistent with the previous work.

In just a few lines it's possible to make a wide variety of activity
diagrams for organic and inorganic species. Try it for your favorite
system!


\section{Where to go from here}

You can explore the package documentation through R's help system;
just type \texttt{help.start()} at the command line and select CHNOSZ
in the browser window that comes up. Besides this document, there
are other vignettes on topics of relative abundances and chemical
activities of proteins, relative protein stability in a hot spring,
and Gibbs energy minimization.

As a more visual way to get an idea of the types of calculations available
in CHNOSZ, try running the examples in the help files for individual
functions. A good one to try out might be \texttt{diagram()}; you
can run all of the examples there with a single command:

<<example_diagram,eval=FALSE>>=
example(diagram)
@

Or you can use the following to run \emph{all} of the examples provided
in the documentation for the package. You will see a lot of text fly
by on the screen, as well as a variety of plots. The examples will
take about 5--10 minutes to run, depending on your machine. 

<<examples,eval=FALSE>>=
examples()
@

There are even more examples that can be accessed by \texttt{demo()}
(or \texttt{demos()} to run all of them):

<<demo,eval=FALSE>>=
demo("findit")
@

If you want to add to or modify the thermodynamic database, read the
instructions at the top of the help page for \texttt{thermo}:

<<helpthermo,eval=FALSE>>=
help(thermo)
@

Have fun!


\section{More activity diagrams}

The following pages contain activity diagrams created with more complex
series of command to show what ``real life'' usage of CHNOSZ might
look like. Also, examples using some functions not covered above (buffer,
revisit, findit) are included. To save space, the output from the
commands (other than the plots) is hidden.

\clearpage


\subsection{Protein Eh-pH}

This plot, showing the relative stabilities of extracellular $\alpha$-amylase
on an Eh-pH diagram is similar to Fig. 13 of \citet{DLH06}. The basis(``CHNOSe'')
includes an electron, so Eh can be used as a variable. The arguments
in \texttt{diagram()} include setting the balanced component to $\mathrm{CO_{2}}$;
\texttt{normalize=FALSE} is left as the default so that whole protein
formulas, not normalized by protein length, are used in the calculations.

Note that \texttt{add.obigt()} could be added at the very beginning
to load the properties of the methionine sidechain group used by \citet{DLH06}
(the default values in the current version of CHNOSZ are taken from
\citealp{LD12}), but for this example it doesn't make a whole lot
of difference.

Toward the end of the script, points are added for Eh-pH values from
soils \citep{BKM60} and hot springs in Yellowstone ($+$) \citep{SWM+05}
and Iceland ($\blacktriangle$) \citep{SA02}. The symbols identifying
the latter two sources were swapped in the figure caption of \citet{DLH06}.
Finally legends are drawn to identify the lines and symbols. The \texttt{describe.property()}
function of CHNOSZ is used to generate the temperature notation (italic
$T$, and degree sign in the units).

\setkeys{Gin}{width=0.5\textwidth}

<<amylaseplot, fig=TRUE, results=hide, width=5, height=5>>=
basis("CHNOSe")
basis(c("NH3", "H2S"), c(-6, -3))
species(c("AMY_BACSU", "O08452_PYRFU"))
a <- affinity(pH=c(2, 10), Eh=c(-1, 1))
diagram(a, balance="CO2", lwd=2, fill=NULL, cex.names=1.5)
a <- affinity(pH=c(2, 10), Eh=c(-1, 1), T=100)
diagram(a, balance="CO2", lwd=2, fill=NULL, col="red", names=NULL, add=TRUE)
water.lines()
water.lines(T=398.15, col="red")
BKMdat <- read.csv(system.file("extdata/cpetc/BKM60_Fig7.csv", package="CHNOSZ"))
points(BKMdat$pH, BKMdat$Eh, pch=20)
points(c(8.5, 6.2, 4.2, 6.7), c(0.018, 0.223, 0.022, 0.067), pch=3, col="red")
points(c(8.24, 7.17, 8.11), c(-0.44, -0.41, -0.49), pch=17, col="red")
ltext <- c(describe.property("T", 25), describe.property("T", 100))
legend("topright", legend=ltext, lty=1, col=c("black", "red"))
ltext <- c("soils [BKM60]", "yellowstone [SWMP05]", "iceland [SA02]")
legend("bottomleft", legend=ltext, pch=c(20, 3, 17))
@

\setkeys{Gin}{width=1.0\textwidth}

\clearpage


\subsection{Subcellular proteins}

Localizations and abundances of proteins from YeastGFP \citep{HFG+03,GHB+03}
are used here to calculate an abundance-weighted average of amino
acid compositions of proteins in different subcellular compartments
of yeast. The relative stabilities of these 23 model proteins are
calculated as a function of the logarithm of oxygen fugacity ($\log f_{\mathrm{O_{2}}_{\left(g\right)}}$).
This figure is similar to Fig. 3 of \citet{Dic09}. Differences in
positions of the lines can be attributed to updated parameters for
the methionine sidechain group \citep{LD12} that are used in the
current version of CHNOSZ.

\setkeys{Gin}{width=0.7\textwidth}

<<yeastplot, fig=TRUE, results=hide, width=7, height=5>>=
locations <- yeastgfp()
gfp <- yeastgfp(locations)
aa <- more.aa(gfp$protein, "Sce")
for(i in 1:length(locations)) {
  avgaa <- aasum(aa[[i]], gfp$abundance[[i]], average=TRUE, protein=locations[i])
  add.protein(avgaa)
}
basis("CHNOS+")
species(locations, "Sce")
a <- affinity(O2=c(-82, -65))
e <- equilibrate(a, loga.balance=0, normalize=TRUE)
mycolor <- topo.colors(length(locations))
diagram(e, names=locations, ylim=c(-5, -3), legend.x=NA,
  col=mycolor, lwd=2)
dp <- describe.property(c("T", "P"), c(25, 1))
db <- describe.basis(ibasis=(1:6)[-5])
legend("topright", legend=c(dp, db), bty="n")
@

\setkeys{Gin}{width=1.0\textwidth}

Notable features include: proteins in the early Golgi, endoplasmic
reticulum (ER) and vacuole are chemically the most stable relative
to those in other locations in the cell; proteins in the vacuole are
very oxidized; the stability of proteins decreases going from early
Golgi to Golgi to late Golgi; the least stable proteins (or most energetic)
are found in the microtubules and bud neck.

\clearpage


\subsection{Buffers}

Chemical activity buffers permit calculating the activities of basis
species from the activities of the species of interest, instead of
the other way around. In CHNOSZ there are two ways to perform buffer
calculations: by assigning the name of a buffer in \texttt{thermo\$buffer}
to the basis species, or by using the \texttt{what} argument of \texttt{diagram()}
to solve for the activity of the indicated basis species. The former
method is more versatile (multiple activities can be buffered, e.g.
both $\mathrm{S_{2}}$ and $\mathrm{O_{2}}$ by pyrite-pyrrhotite-magnetite,
and the buffers have effects on equilibrium activity diagrams) while
the latter is more convenient (no need to set up the buffer in \texttt{thermo\$buffer}
-- the buffers come from the species of interest). The plot below,
based on Fig. 6 of \citet{SS95}, shows values of $\log f_{\mathrm{H_{2}}_{\left(g\right)}}$
buffered by minerals and in equilibrium with different activities
of organic species that are calculated using these two methods.

\begin{small}
\setkeys{Gin}{width=1.0\textwidth}

<<bufferplot, fig=TRUE, results=hide, width=8, height=3.5>>=
layout(matrix(1:2, nrow=1), widths=c(2, 1))
b.species <- c("Fe", "CO2", "H2O", "N2", "H2", "H2S", "SiO2")
b.state <- c("cr1", "gas", "liq", "gas", "gas", "aq", "aq")
b.logact <- c(0, 1, 0, 0, 0, 0, 0)
basis(b.species, b.state, b.logact)
xlim <- c(0, 350)
thermo.plot.new(xlim=xlim, ylim=c(-4, 4), xlab=axis.label("T"), ylab=axis.label("H2"))
bufferline <- function(buffer, ixlab) {
  basis("H2", buffer)
  a <- affinity(T=xlim, P=300, return.buffer=TRUE, exceed.Ttr=TRUE)
  lines(a$vals[[1]], a$H2)
  text(a$vals[[1]][ixlab], a$H2[ixlab], buffer)
}
bufferline("FeFeO", 20)
bufferline("QFM", 38)
bufferline("PPM", 102)
bufferline("HM", 51)
basis("H2", 0)
for(logact in c(-6, -10, -15)) {
  species(c("formaldehyde", "HCN"), logact)
  a <- affinity(T=xlim, P=300)
  d <- diagram(a, what="H2", lty=c(2, 3), add=TRUE)
  text(a$vals[[1]][13], mean(sapply(d$plotvals, c)[13, ]), logact)
}
plot.new()
legend("topleft", legend = c(describe.property("P", 300), describe.basis(ibasis=c(2,4)),
  "minerals", "HCN", "formaldehyde"), lty=c(NA, NA, NA, 1, 2, 3), bg="white")
@

\setkeys{Gin}{width=1.0\textwidth}
\end{small}

\clearpage


\subsection{revisit}

\texttt{revisit()} computes some summary statistics (default is the
coefficient of variation) about the equilibrium chemical activities
of species calculated by \texttt{equilibrate()}. In this example,
the coefficient of variation of the activities of the amino acids
is plotted as a function of $\log f_{\mathrm{O_{2}}_{\left(g\right)}}$
and $\log a_{\mathrm{CO_{2}}_{\left(aq\right)}}$.

\begin{small}
\setkeys{Gin}{width=0.35\textwidth}

<<revisit, fig=T, results=hide, width=4, height=4>>=
basis("CHNOS")
species(c("isoleucine", "tyrosine", "glutamic acid", "methionine", "aspartic acid"))
a <- affinity(CO2=c(2, 5), O2=c(-72, -70))
e <- equilibrate(a, balance=1)
r <- revisit(e)
title(main=paste("CV minimum =", round(r$optimum, 2)))
@

\setkeys{Gin}{width=1.0\textwidth}
\end{small}

The \texttt{balance=1} in \texttt{equilibrate()} means the relative
stabilities are calculated using the formation reactions written per
mole of amino acid (not conserving e.g. $\mathrm{CO_{2}}$ which would
be the default behaviour in this system). Note the star showing the
conditions where the coefficient of variation is minimized. Let's
look at the fractional equilibrium abundances of the amino acids at
these conditions.

\begin{small}
\setkeys{Gin}{width=0.6\textwidth}

<<revisit_alpha, fig=T, results=hide, width=6, height=4>>=
basis(c("CO2", "O2"), c(r$x, r$y))
a <- affinity()
par(mfrow=c(1, 2))
e <- equilibrate(a, balance=1)
d <- diagram(e, alpha=TRUE, names=aminoacids(1, species()$name))
plot.new()
legend("topleft", describe.basis(basis()), bg="white")
@

\setkeys{Gin}{width=1.0\textwidth}
\end{small}

This hypothetical metastably equilibrated mixture has very little
methionine and aspartic acid. Can we find where the relative abundances
of the amino acids have a more even distribution?

\clearpage


\subsection{findit}

\texttt{findit()} performs a gridded search, on successively smaller
hypercubes, for the conditions that optimize one of the statistics
returned by \texttt{revisit()}. In this example containing five species,
with \texttt{balance=1}, four compositional variables are the maximum
that can be considered ($\mathrm{NH_{3}}$ is excluded because it
has a reaction coefficient of 1 in all of the formation reactions);
the hypercube in this case is a tesseract. The grid resolution is
10, so the equilibrium chemical activities of the amino acids are
computed at $10^{4}$ discrete combinations of chemical potential
of the basis species with each iteration.

\begin{small}
\setkeys{Gin}{width=0.4\textwidth}

<<findit, fig=T, results=hide, width=4, height=4>>=
basis("CHNOS")
species(c("isoleucine", "tyrosine", "glutamic acid", "methionine", "aspartic acid"))
f <- findit(list(CO2=c(-5, 5), O2=c(-85, -65), H2S=c(-10, 5), H2O=c(-10, 0)), 
  niter=5, res=10, balance=1)
@

\setkeys{Gin}{width=1.0\textwidth}
\end{small}

After 5 iterations, what are the fractional equilibrium abundances
of the amino acids? Note that, during its operation, \texttt{findit()}
updates the activities of the basis species so we don't have to set
them manually.

\begin{small}
\setkeys{Gin}{width=0.6\textwidth}

<<findit_alpha, fig=T, results=hide, width=6, height=4>>=
a <- affinity()
par(mfrow=c(1, 2))
e <- equilibrate(a, balance=1)
d <- diagram(e, alpha=TRUE, names=aminoacids(1, species()$name))
plot.new()
legend("topleft", describe.basis(basis()), bg="white")
@

\setkeys{Gin}{width=1.0\textwidth}
\end{small}

We found a combination of chemical activities of basis species that
lowered the variation of the equilibrium activities of the amino acids.
Woohoo!

A more complete analysis might include more amino acids, different
balancing constraints or chemical activity (or even temperature) limits,
and would probably increase \texttt{niter} and/or \texttt{res} in
\texttt{findit()}. Note that the results of a high-dimensional optimization
such as this one may be misleading if the resolution is not high enough,
or the initial hypercube is too small. Even after a successful simulation,
it remains a matter of discussion what the optimized chemical activities
of the basis species and the amino acids signify.


\section{Document information}

Revision history:
\begin{itemize}
\item 2010-09-30 Initial version
\item 2011-08-15 Add \texttt{browse.refs()}; modifying database hint changed
to \texttt{help(thermo)}
\item 2012-06-16 Add ``More activity diagrams''
\item 2015-05-14 Add warning about internal consistency of thermodynamic
data
\end{itemize}
R session information:

<<session_info>>=
sessionInfo()
@

\bibliographystyle{plainnat}
\bibliography{vig}

\end{document}
