text:
\section*{Introduction}

Drug Discovery~\cite{Lyne2002} is a process that includes several phases, from virtual \emph{in silico} simulations to \emph{in-vitro} and \emph{in-vivo} experimentation.
Molecular Docking~\cite{Meng2011} (MD) is an important step of  the  drug discovery process which aims at calculating the preferred position and shape of one molecule to a second when they are bound to each other.
The computational resources needed to address MD related tasks are usually quite large, as problems of this kind  involve several degrees of freedom and rapidly growing dimensionality.

In recent years, the field of Quantum Computing~\cite{nielsen_chuang_2010} (QC) has undergone significant developments  from  both hardware and algorithms points of view.
In particular, we are living in the so-called NISQ~\cite{Preskill2018quantumcomputingin} era in which it is not clear whether the quantum devices currently available are actually capable of producing better or comparable results compared to classical methods.

In this work, we explored the possibility of  using a QC technique called Quantum Annealing~\cite{7055969,babej2018coarsegrained,Marchand2019} (QA)  in the virtual screening phase of MD.
In particular, we focused on the Molecular Unfolding (MU) process, which is the first step in geometric molecular docking techniques.
MU aims at finding the molecular configuration that maximizes its volume, or equivalently, that maximizes the internal distances between atoms inside the molecule.

Our aim was to develop a  Quantum Molecular Unfolding approach and execute it on the   latest quantum annealing hardware (D-Wave Advantage and 2000Q) in order to  understand the capabilities of these devices and compare their performances to  state-of-art molecular docking methods.

The remainder of the paper is structured as follows.
In the first section, the topic of molecular docking is introduced, followed by a section where the problem of Molecular Unfolding is discussed. In  Sec.3, concepts of quantum annealing are illustrated while in  Sec.4  a binary optimization formulation of the problem is given.  Later in the fifth section, the dataset and   pre-processing steps employed before the unfolding phase are presented. This is followed in Sec.6 by  a review  of  state of the art classical algorithms employed to compare solutions obtained with  the quantum devices. Finally, in the last two sections, results regrading the embedding  on the QPUs are presented  as well as a detailed performance comparison between the different methods.

\section{Preliminaries on Molecular Docking}

In-silico approach to drug discovery~\cite{Lyne2002} can be viewed as a multi-tiered process  that encompasses several sequential computational techniques  with the aim of screening virtual libraries of the order of billions of compounds for the most suitable molecules to forward to later experiments.

Molecular Docking~\cite{Meng2011} (MD) is an essential step in virtual screening which is used  to simulate the atomic interactions of  a ligand  inside a protein binding site, in order to highlight their possible biochemical reactions and predict whether a stable complex could be formed.
\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.45\textwidth]{./pictures/lock_and_key_reale.JPG}
	\caption{Visual representation of Molecular Docking  for a ligand inside a
		protein pocket, similar to the relation between locks and keys.}
	\label{fig:background:Lock and Key}
\end{figure}
\\
\\

The MD process is divided in two main tasks:
\begin{itemize}
	\item Detection of three-dimensional \emph{poses} of the ligand, i.e. valid ligand conformations, positions and orientations inside the active site  of the protein (usually  called \emph{pocket}).

	\item Ranking of the poses  via a scoring function.  Usually the lower  the docking score, the better  the resulting  \emph{binding affinity}.
\end{itemize}

Since  electrons inside atoms are repulsive to each other, their interaction  affects both molecular shape and  reactivity. Therefore,  by controlling the molecular shape, it is possible to  predict the ligand  reactivity to a protein's pocket as well as the energy cost of such configuration. This direct connection between molecular conformation and binding affinity is what enables geometrical scoring functions~\cite{Gadioli2020}.

In our approach, the docking process considers the pocket as a rigid structure, while the ligand is a flexible set of atoms. Furthermore, from a strictly geometrical interpretation, the ligand can  be seen as a set of chemical bonds (or edges) with a fixed length, where only a  subset of edges are  \textit{rotatable}.

Rotatable bonds (or torsionals) are bonds that split the molecule in two nonempty disjointed fragments, when virtually removed.  Consequently, such fragments  can rotate independently from each other around the axis of the rotatable bond.
These  concepts are graphically reported in Fig.\ref{fig:background:Fragments and rotatable bonds.}, where the rightmost rotatable bond is splitting the molecule in left and right fragment.
Finally,  note that the number of fragments in a molecule is usually double the number of  rotatables.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.4\textwidth]{./pictures/fragments_and_rotables.JPG}

	\caption{Fragments and rotatable bonds}
	\label{fig:background:Fragments and rotatable bonds.}
\end{figure}
\FloatBarrier

In general, it is possible to point out three main phases in MD: \emph{Ligand expansion}, \emph{Intial Placement} and\emph{ Shape Refinement} inside the pocket.

Here we focus on  the ligand expansion phase, which is essential for improving  docking. In fact,  an initial pose of the ligand that is set a priori  may introduce shape bias  affecting  the final quality of the docking.  Molecular Unfolding (MU) is the technique used for removing such initial bias  by expanding the ligand to an unfolded shape.

\section{Molecular Unfolding Problem Definition}
The objective of the MU problem is to find the unfolded shape of the ligand or,  in other terms, the torsion configuration that maximises the molecular volume expressed as   total sum of  internal distances between pairs of atoms in the  ligand. The starting point is the folded molecule defined by a  set of atoms together with its  fixed and rotatable bonds.
We can assign to each  rotatable bonds $t_{i}$ a variable corresponding to the angle $\theta_{i}$  which is responsible for fragment rotation.
As a convention, we  identify the ordered set of rotations by a vector:
\begin{equation}\label{eq:t_vec}
t=\begin{bmatrix} \theta _{1}, \theta _{2},  \ldots  , \theta _{n} \end{bmatrix}
\end{equation}
where each torsion  $\theta_i$ around the bond's axis can assume  values $\left[ 0, 2\pi \right)$.
Given a molecule, it is possible to construct the associated graph as shown in Figure~\ref{fig:test_molecule2} where  edges represent bonds (which  can not contract).  Torsional bonds are depicted in red in the image below.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.4\textwidth]{pictures/aspirina3d.png}\ \ \
	\includegraphics[width=0.45\textwidth]{pictures/orig_aspirina2_num_red.png}
	\caption{ Illustration of a simple organic molecule. The molecular structure presents 4 rotatable bonds, highlighted  in red in the graph representation.}
	\label{fig:test_molecule2}
\end{figure}
\FloatBarrier

The objective for the MU problem can be expressed in mathematical terms as follows: given a molecule,  find the  torsional configuration
\begin{equation}
\centering
\vect{t^{unfold}} = [\theta_1^{unfold},\, \dots , \theta_M^{unfold}]
\end{equation}

such that the following quantity is maximised

\begin{equation} \label{eq:D}
\centering
D({\bf \vect{t}}) = \sum_{\substack{a,b \in \mathcal{M} \\ a \neq b}}  D_{ab}({\Theta})^2\,
\end{equation}

The function $ D_{ab}( \Theta) $ denotes the distance between two different atoms $ a $ and  $ b $ belonging to the molecule $ \mathcal{M} $ while $ D( \bf{ \vect{t} }) $ is  the sum of all square  distances between  atoms.

The rationale behind this choice is that the objective function is simple and relies just on the geometry of the molecule.  Each distance  $ D_{ab}( \Theta) $  only depends  on the angles induced by those rotatable bonds that appear in the shortest path connecting atom $a$ to $ b$ in the graph.

It is worth noticing that  it is not necessary to calculate all the pairwise distances that are expressed in equation \ref{eq:D}. In fact, it is possible to apply  simplifications to  $ D_{ab}(\Theta) $, as expressed by the following two conditions:

\begin{condenum}
	\item The shortest path between two atoms must contain at least one torsional. This is due to the fact that if two atoms belong to the same rigid fragment,  their relative position  never changes. \label{cond: conditions4Distance1}
	\item The shortest path connecting two atoms must have at least three edges. This is a consequence of the fact that atoms connected by only one or two  bonds  have always the same distance, even when  bonds are rotated. \label{cond: conditions4Distance2}
\end{condenum}

A useful mathematical representation for rotations are rotation matrices, $R(\theta_i)$, which are function  of the  angle $ \theta_i $  associated to torsion $ T_i $ .
If two atoms are connected via multiple torsional, the rotation matrix  $R({\Theta})= R(\theta_i,\theta_j, \dots ,\theta_k)  $,
characterizing their relationship is
an ordered multiplication of single torsion rotation matrices
\begin{equation} \label{eq:R}
R({\Theta})=R(\theta_i,\theta_j, \dots, \theta_k) = R(\theta_i)\times R(\theta_j) \times ... \times R(\theta_k)
\end{equation}
Each rotation  $R(\theta_i)$ is expressed via a $4\times4$ rotation matrix  acting on atomic coordinates.  The extremes of a bond connecting  two atoms are identified by the coordinates of the atoms themselves.

Consider now the distance $ D_{ab}( {\Theta}) $ and
the initial positions of atoms  $ a $ and $ b  $   which are  identified with  $\vec{a_0}$
and
$\vec{b_0}$ respectively.
Relative positions are obtained by fixing one of the two atoms while  rotating the other as follows
\begin{equation}
\vec{a}= \vec{a_0} \ \, \ \ \ \ \  \  \vec{b}=R({\Theta}) \vec{b_0}
\end{equation}
Hence, we can rewrite the single contributes  in equation \ref{eq:D} using the Euclidean distance as
\begin{equation}\label{eq:Euclidean_dist}
D_{ab}( {\Theta})^2 = \| \vec{a_0} - R({\Theta}) \vec{b_0}\|^2
\end{equation}
Figure~\ref{fig:distance_according_to_rotation} shows how the    relative positions change as function of the angles assumed by the two rotatable bonds.
\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.30\textwidth]{pictures/distanza_by_rotazione.png}
	\caption{ Visual representation of the action of two rotations $\theta_1$ and $\theta_2$ on the distance expressed in eq.
		\ref{eq:Euclidean_dist}.}
	\label{fig:distance_according_to_rotation}
\end{figure}

\FloatBarrier

\section{Quantum Annealing}
Quantum annealing (QA) is a meta-heuristic technique for addressing challenging combinatorial problems.
It can be seen as a  variation of the \emph{Simulated Annealing (SA)}~\cite{Kirkpatrick1983} algorithm where quantum fluctuation, instead of thermal effects, are employed in the exploration of the configuration space.
Quantum annealing devices, like those   implemented  with \emph{superconductive qubits} by the Canadian  company \textit{D-Wave}, are able to  solve  combinatorial  problems expressed in   Quadratic Unconstrained Binary Optimization (QUBO) terms.
A generic QUBO problem  is defined as
\begin{equation}\label{eq: Ising_problem}
O(x)=\sum_{i}h_i {x}_i+ \sum_{i>j}J_{i,j} x_i x_j
\end{equation}
where $ x_i \in \left\lbrace 0,1 \right\rbrace  $  are binary variables while  $ h_i $ and  $ J_{ij} $ are parameters whose values encode the optimization task  in such a way that the  minimum of $ O(x) $ represents the solution to the optimization problem.

In general,  a given  optimization problem can be easily translated  in binary combinatorial terms as a a high-order quadratic unconstrained binary optimization (HUBO) problem.
\begin{equation}\label{eq: Ising_problem}
O_{Hubo}(x)=\sum_{i}\alpha_i {x}_i+ \sum_{i,j}\beta_{i,j} x_i x_j+  \sum_{i,j,k}\gamma_{i,j,k} x_i x_j x_k +...
\end{equation}

Fortunately, it is always possible to convert  HUBOs into QUBOs.  The trick used is to add new ancillary binary variables and  substitute  high-order terms with  sums of quadratic expressions (made up with   original binary variables and  new ones) in order to preserve  local and global minima~\cite{mengoni2020breaking}.

As an example, a cubic term like $ \pm x_1\cdot x_2\cdot x_3 $ can be divided using the ancillary binary variable $a_1$ as
\begin{equation}\label{eq: make_quadratic_rel}
\pm x_1 \cdot x_2 \cdot x_3 \rightarrow  \pm a_1 x_3 + 2 (x_1 x_2 - 2 x_1 a_1 - 2 x_2 a_1 + 3 a_1 )
\end{equation}
I practice, one can decide to build a custom function for quadratisation or exploit the functions present in D-Wave software offer such as \textit{make\_quadratic}.

\section{HUBO Formulation}
As a first step, it is necessary to formulate  the MU problem in binary optimization terms in order to employ a quantum annealing approach.

Let us consider a discretisation of the angle $ \theta_i $ associated to a rotatable bond $T_i$ into $d$ possible values
\begin{equation}
\theta_i= [\theta_i^1 , \theta_i^2 ,  \theta_i^3,..., \theta_i^d ]
\end{equation}
As a consequence,  continuous functions of the angles like sine and cosine can be  discretised as well in $d$  values
\begin{align}
\sin(\theta_i)& = [\sin(\theta_i^1) , \sin(\theta_i^2) ,  \sin(\theta_i^3) ,..., \sin(\theta_i^d)   ]\\
\cos(\theta_i)& = [\cos(\theta_i^1) , \cos(\theta_i^2) ,  \cos(\theta_i^3) ,..., \cos(\theta_i^d)   ]
\end{align}
Moreover, since torsional angles can assume  only one value at a time, it is necessary to   associate each $\theta_i$ with only one value among all the possible $\theta_i^k$. A {\emph{One-Hot Encoding}} strategy is adopted for each  torsional $T_i$, to which are  assigned  binary variables $x_{ik}$,  with  $1 \leq k \leq d$,  such that
\begin{equation}\label{eq:system_classical_constraint_onehot}
x_{ik} = \begin{cases}
1 & \text{if } \theta_i= \theta_i^k; \\
0 & \text{otherwise.}
\end{cases}
\end{equation}
Definition \ref{eq:system_classical_constraint_onehot} tells us that only one among the binary variables can be assigned to a one or truth value. This means that  their sum  is always equal to one:
\begin{equation}\label{eq:sum_classical_constraint_onehot}
\sum_{k=1}^d x_{ik} =1
\end{equation}
Note that equation \ref{eq:sum_classical_constraint_onehot} must hold for all  rotatables. In order to convert such algebraic constraint into QUBO, the following term has to be considered:
\begin{equation}\label{eq: squared-sum-hard}
\sum_{i} \left( \sum_{k=1}^d x_{ik} -1 \right)^2
\end{equation}
Constraint of eq. \ref{eq: squared-sum-hard} is  usually  called a \emph{hard constraint} because it is mandatory to satisfy  this term i.e. minimize it in the combinatorial optimization.

Regarding those functions that appear in the rotation matrix $R({\theta_i}) $   such as
$\sin(\theta_i)$  and $ \cos(\theta_i) $, they can be  expressed  as
\begin{align}\label{eq: one-hot constraint}
\sin(\theta_i)&= \sum_{k=1}^d \sin(\theta_i^k)\  x_{ik} \\
\cos(\theta_i)&= \sum_{k=1}^d \cos(\theta_i^k)\  x_{ik}
\end{align}
Therefore,  a rotation matrix $R({\theta_i}) $  becomes  a  function of all the binary variables  $x_{ik}$ needed to represent the angle $\theta_i$
\begin{equation}
R({\theta_i}) =R(x_{i1}, x_{i2}, ..., x_{id})
\end{equation}
For a generic rotation we have
\begin{align}
\nonumber
R({\Theta})& = R(\theta_i)\times R(\theta_j) \times ... \times R(\theta_k)\\
\nonumber
&=R(x_{i1}, x_{i2}, ..., x_{id})\times R(x_{j1}, x_{j2}, ..., x_{jd})  \times ...\\
&... \times R(x_{k1}, x_{k2}, ..., x_{kd})
\end{align}

Note that, for what concerns the granularity of the rotations, in order to obtain a precision of $\Delta\theta_i=~\theta_i^{k+1}-~\theta_i^k $ in  the representation of angle $\theta_i$,  the  number of variables $x_{ik}$ needed for each torsion is given by
\begin{equation} \label{eq: discretization}
d =\dfrac{2 \pi}{\Delta\theta_i}= \dfrac{2 \pi}{\theta_i^{k+1}-\theta_i^k}
\end{equation}
Given a number $ M $  of torsional bonds, the total number of binary variables is
\begin{equation}
n=d\times M=\dfrac{2 \pi}{\Delta\theta_i}\times M
\end{equation}
The general form of the HUBO optimization function can be written as
\begin{equation}\label{eq: HUBO}
O(x_{ik})=A_{const}\sum_{i} \left( \sum_{k=1}^d x_{ik} -1 \right)^2 -  \sum_{\substack{a,b}} D_{ab}( {\Theta})^2
\end{equation}
In the last equation, the second term identifies the \emph{optimization term}, which has a minus sign in front because it should be maximized (while the whole expression is minimized). The first term instead represent the \emph{hard constraint}, where the parameter $ A_{const} $ is a penalty scalar that  modulates its strength.
Finally, it is worth noticing that, given a molecule with $M$ torsionals, the highest order term appearing in eq.\ref{eq: HUBO} is $2M$, where the number 2 comes from the  squared distances.

\section{Dataset and Pre-processing}
\subsection{Ligand Dataset}
The original ligand dataset in our possession is made up of 118 molecules, with a number of atoms ranging form 20 to more that 120   and a number of torsionals going from few units to about 50.

As shown in Figure \ref{fig:proof_expansion} below,  it is possible to see that 60\% of the ligands in our database  present  10-20 fragments, which means 5-10 rotatable bonds. The number of atoms is peaked about the   40-50  interval.
For this reason, we decided to concentrate our research only on molecules that had at most  $12$ rotatable bonds and an increasing number of atoms going from $20$ to $50$.

Since our study aimed at highlighting the change in problem complexity  by increasing the number of torsionals, we considered only a subset of torsionals at a time, out of 12 possible rotatable bonds.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.45\textwidth]{pictures/number_of_fragments.JPG}

	\includegraphics[width=0.45\textwidth]{pictures/number_of_atoms.JPG}

	\caption{Analysis of the distribution for number of atoms and number of fragments.}
	\label{fig:proof_expansion}
\end{figure}
\FloatBarrier

\subsection{Pre-processing Phase}
The pre-processing   starts with the parsing of the information about the  molecule via analysis of its .MOL2 file.
A graph containing a three-dimensional description of the ligand together with bond type (bonds could be single, double, triple, amide or aromatic) is constructed.

Two  simplifications related to chemical-physical properties of ligands  are applied in the  pre-processing phase: one is due to the fact that only single and aromatic bonds can be  rotatable, while the other concerns the removal of terminal hydrogen atoms whose distances with the rest do not affect the total sum of internal distances.
Note that the removal of  terminal hydrogens may trigger the generation of conformations where the final position of  atoms do not satisfy their minimum Van Der Waals distance. Such invalid shapes are discarded in post-processing.

After this,   the  graph associated to the molucle is divided into fragments as shown in Figure \ref{fig:Rotables_influence_set_example}.
At this stage, the betweenness centrality~\cite{Freeman1977}, defined in  equation \ref{eq: betweeness},  is employed to obtain an ordering of the atoms within the  graph. By consequence, this procedure also specifies an ordering of rotatable bonds.
\begin{equation}  \label{eq: betweeness}
\mathrm{ Betweeness \ centrality\ of\ atom\ \textit{v} : }\ \
\sum_{v \neq s \neq t} \frac{\sigma_{st}(v)}{\sigma_{st}}
\end{equation}
Equation \ref{eq: betweeness} is calculated for each atom $v$ as the sum of all the possible shortest paths connecting atoms $s$ and $t$ which  also cross $v$, divided by the total number of shortest paths between $s$ and $t$.
The atom with the greatest centrality is chosen as  centre of the graph and origin of the torsional ordering.

Each fragment is identified as a set of atoms influenced by the same set of rotatable bonds,  called \emph{rotatables influence set}.
The definition of rotatables influence set is the following
\begin{equation}  \label{eq: inf_set}
\mathrm{ Rotables\ influece\ set\ :\ }\ \
I_{s} = E_{C_{a,a_{k}}} \cap E_{R}
\end{equation}
$C_{a}$ = atom\ with\ greatest\ betweeness\ centrality \\
$ E_{R}$ = Rotatable\ bonds\ \\
$E_{C_{a},a_{k}}$ = Bonds\ on\ the\ shortest\ path\ $\sigma_{C_{a},a_{k}}$ \\ \\
Atoms of the same fragment are mapped to a unique shortest path, which is composed by a set of edges belonging to the influence set.
\begin{figure}[htbp]
	\centering
	\includegraphics[width=0.45\textwidth]{pictures/frammenti_con_angoli.png}\ \ \
	\caption{Atoms belonging to the same fragment, depicted using colors, are affected by the same  rotatables influence set.}
	\label{fig:Rotables_influence_set_example}
\end{figure}

Finally, contributions appearing in equation \ref{eq:Euclidean_dist} are obtained for each fragment, where distances between  atoms  are calculated only once and  respecting  conditions \ref{cond: conditions4Distance1} and \ref{cond: conditions4Distance2}. Distances between atoms belonging to  the same fragment are fixed, hence  not considered.

\section{Model Development}

\subsection{Distance Simplification}

The first simplification   concerns the calculation of  contributions in equation \ref{eq:Euclidean_dist}, as the number of terms  increases quadratically with the number of atoms involved.
Hence,  instead of measuring the Euclidean distance between each atom $a$ in  fragment $A$ and  each  atom $b$ in  fragment $B$, only distances between  different fragments have been considered. In particular, each fragment was identified with the coordinates of the median atom inside the fragment.

Partial contributions for single fragments are summed  up together in the HUBO as in  equation \ref{eq: HUBO}.
Note that such optimization function  is a non-linear symbolic expression where  symbolic variables representing torsional angles  appear in the rotation matrix in the form of   trigonometric functions as defined in equation \ref{eq: one-hot constraint}.

\subsection{Coarse-grained Rotations}
The number of variables needed to express the MU problem in combinatorial optimization form depends on the discretisation of the angles, as shown in eq. \ref{eq: discretization}. It is therefore essential to understand how to chose an optimal  granularity for the torsionals.
In this regard, previous knowledge about the dataset was used  to obtain approximate  solutions achieving an acceptable quality.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.48\textwidth]{pictures/granularity_degradation_effect.JPG}\ \ \
	\caption{Analysis of the  degradation by transitioning to a coarse-grained rotation. }
	\label{fig:45vs1}
\end{figure}
Our dataset  is displayed in Figure \ref{fig:45vs1}, where molecules are  ordered by increasing volume. The vertical axis shows the percentage of degradation in the molecular volume when the unfolding is made with an angle greater than 1 degree. The lines drawn in the plot identify granularities of 45, 20, 10, 5 and 2 degrees. The surprising result is that the degradation is always below 4\% for  rotations at 45 degrees with average around 1.7\%.

For this reason, even though we introduced the possibility of using different angular discretisations, as shown in eq. \ref{eq: discretization},  we focuses on  a granularity of  45 degrees.
The combinatorial optimization  would only require  8 binary variables  for each rotatable bond to identify its angle while  keeping the final quality loss acceptably low.
\subsection{Threshold Approximation}

The generation of the HUBO phase  is one of the most computationally intensive, since we are using a symbolic framework. The challenges that arise  are strictly related to the enormous amount of memory used and the unparalleled implementation of the symbolic engine.
Hence, a  straight derivation of the HUBO becomes unrealistic when increasing problem size.

For this reason, a  last simplification was introduced, which  consists in eliminating the  optimization terms in \ref{eq: HUBO}  (only those not related to the hard constraints) whose   coefficient is strictly less than a certain threshold value.

The effect of such threshold approximation (also called \textit{chop}) was to speed up the computation of the HUBOs, as shown in Figure \ref{fig: HUBOtime},  and decrease the number of terms in the QUBO function in order to match the quantum annealing hardware specifications and limitations such as number of qubits (i.e. binary variables) and connectivity (i.e. quadratic terms) of the problem.
Our study doesn't go beyond the  8 rotatable bonds as  HUBO construction involving 10 torsionals exceeded the time limit  with the current setup and software.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.45\textwidth]{./pictures/unchopped/avg_creation_time.jpg}
	\caption{Average creation time for HUBOs. We were able to generate HUBOs without the threshold approximation, up to 5 torsionals. In the case of 6 and 8 rotatable bonds, the threshold was increase to the vale of 30 and 60 respectively. }
	\label{fig: HUBOtime}
\end{figure}
\FloatBarrier

\section{Altenative Search Strategies}
In the following, classical optimization strategies  used to compare   quantum annealing results  are presented.

\subsection{Random Search}
Random Search~\cite{random_search}  algorithm works by picking up randomly the candidate solution from a uniform distribution.
A big strength  of this method is  the high parallelization and the possibility of using  Bayesian techniques or probability distributions that are more suitable for a given problem.

\begin{algorithm}
  \caption{Random Search}
\label{Random Search}

	\begin{algorithmic}[1]

\Require $A \ =\ Angles,\ Poly = \ D_{ab}( {\Theta})^2 , \ T = Torsional \ Bonds \ Set$\
\Ensure  ${\gamma_{max}} : \mathbf{max}\ D_{ab}( {\gamma_{max}})^2$ \
\State $  | T | = m $ \
\State $ \Gamma \ =\ A_1 \times A_2 \times... \times A_m$ \ \Comment{Cartesian product of the set of angles - m times.}
\While{ $t \ < \ MaxTime$ } \\
\ \ \ $ \gamma = \ pick\_rand\_uniform(\ \Gamma \ ) \ $\\
\ \ \ $Evaluate\ Poly(\Theta = \gamma)$ \ \Comment{Evaluate assignments and record the optimum.}
\EndWhile
\State return $[\  \gamma_{max}, \ max \ Poly ]$\

	\end{algorithmic}
\end{algorithm}

\subsection{GeoDock Search}

The GeoDock search is a greedy algorithm which considers a single torsional at a time and makes a decision locally for each rotatable bond. The algorithm usually does not lead to an optimal solution but it can reach good approximate solutions in a reasonable time.

This algorithm is the only one that is not employing  symbolic computation. Each bond is rotated independently, starting from the most central and internal bond with respect to betweenness centrality.     Only the conformation that improves the volume is recorded as a solution. The whole process is repeated a number of times or until a plateau in the improvement is reached.

\begin{algorithm}
	\caption{GeoDock Search}

	\label{GeoDock-inspired Search}

	\begin{algorithmic}[1]
		\Require $A \ =\ Angles,\ Poly = \ D_{ab}( {\Theta})^2 , \ T = Torsional \ Bonds \ Set$\

		\Ensure  ${a_{max}} : \mathbf{max}\ D_{ab}( a_{max})^2$ \

		\State $  | T | = m $ \

		\For{ $t_i,\theta_i \ in \  T_{ordered}$ }  \Comment{Ordered by betweenness centrality.}
		\For{ $\ a \ in \ A \ $ } \\
		\ \ \ \ \ \ \ $Evaluate\ Poly(\theta_i \ = \ a )$ \ \Comment{Record intermediate optimums.}

		\EndFor
		\EndFor
		\State return $[\  \underline{a}_{max}, \ max \ Poly ]$\

	\end{algorithmic}

\end{algorithm}

\subsection{Simulated Annealing}

Simulated annealing (SA) is a meta-heuristic technique,  useful in the search of the optimum in large solution spaces.

\begin{algorithm}[H]
	\caption{Simulated Annealing }

	\label{Simulated Annealing Search}

	\begin{algorithmic}[1]
		\Require $A \ =\ Angles,\ Poly = \ D_{ab}( {\Theta})^2 , \ T = Torsional \ Bonds \ Set, \ N_{max} \ = \ Total\ Iterations$\

		\Ensure  ${a_{max}} : \mathbf{max}\ D_{ab}( a_{max})^2$ \

		\State $T\gets T_{max}$
		\State $a_{max}$ = NULL

		\While{$T>T_{min} \  and \  N<N_{max}$}
		\State $next$ = $pickNeighbour(T, best)$
		\State $\Delta E $ = $Poly(next)-Poly(a_{max})$
		\State $r$ = $randomNumberUniform(0,1)$
		\If{$\Delta E < 0$}
		\State $a_{max} = next$  \Comment{The neighbour is accepted as the new optimum.}
		\ElsIf{ $r< e^{ \frac{-\Delta E}{k_b * T}  } $}
		\State $a_{max}$ = $next$   \Comment{Moving to a neighbour.}
		\EndIf
		\State $T$ = $geometricDecrease(T)$
		\EndWhile

		\State return $[\  \underline{a}_{max}, \ max \ Poly ]$\

	\end{algorithmic}

\end{algorithm}
Since we can visualize the solution space as points belonging to a graph,  SA starts with an initial guess $s_c=s_i$ with its evaluation $f(s_i)$. If the neighbour has a lower evaluation, we take it as the current solution $s_c$.
If the neighbour has a higher evaluation, we accept it with a probability written as:
\begin{equation}\label{eq:SAacceptance}
P_{accept} = \min(1,\  e^{-\frac{Poly(neigh.)-Poly(s_c)}{T}})
\end{equation}
where T is the temperature that starts from a high value and decreases  in time according to a  \textit{cooling schedule}.
Since SA solves the QUBO formulations of our problems, this algorithm was taken as  main comparison against Quantum Annealing.

\section{Problem Complexity and Embedding}

\subsection{Resources Estimation}
\begin{figure}[t!]
	\centering
	\includegraphics[width=0.45\textwidth]{./pictures/unchopped/max_lin_terms_unch.jpg}
	\includegraphics[width=0.45\textwidth]{./pictures/unchopped/max_high_terms_unch.jpg}
	\caption{Mean number of linear and non-linear terms in the  HUBOs, obtained varying the number of torsionals and value of the approximation threshold .}
	\label{fig:capitolo6:avg_num_unapprox_hubo}

\end{figure}

The  metrics considered for evaluating the complexity of the problem are   the number of linear and high-order terms appearing in the HUBO  (plotted  in Figure \ref{fig:capitolo6:avg_num_unapprox_hubo}). The number of linear terms follows the relation  $Number\ of\ rotatables \times \ granularity \ of \ rotation$. Recall that  a granularity of 45 degrees was selected , i.e. $d=8$ in equation \ref{eq: discretization}.

The number of non-linear terms grow rapidly  going from 2 to 6 torsions. However, increasing problem size,   the  threshold approximation  allowed us to avoid the exponential growth that these numbers would otherwise follow.

Since the final objective is to be able to embed our problems into a QPU, a second threshold approximations acting on the QUBO has been applied.
The choice of the amount of chop is done by a fine tuning of the threshold, which is lowered until the embedding cannot be performed anymore.

The plots in  Figure \ref{fig:capitolo6:avg_num_approx_qubo} show an increasing trend in the number of QUBO linear terms associated with  the same threshold, for an increasing number of rotatables. Linear terms also  decrease by increasing the threshold when the  number of torsionals is fixed.
Quadratic terms have a similar behaviour, except done for the class with  8 rotatable and a value for the approximating  threshold set at 200, where more complex problems  seem to peak.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.45\textwidth]{./pictures/qubos/qubo_avg_lin_terms.png}
	\includegraphics[width=0.45\textwidth]{./pictures/qubos/qubo_avg_quad_terms.png}
	\caption{Average number of linear and quadratic terms in the  QUBOs, obtained varying the number of torsionals and value of the approximation threshold.}
	\label{fig:capitolo6:avg_num_approx_qubo}
\end{figure}

\subsection{D-Wave Embedding}

In the embedding phase, a heuristic algorithm tries to find the optimal matching between problem resources and physical D-Wave hardware.
Several  embedding are possible, so the one that uses the least physical qubits out of few trials is employed.  It's implicit in this choice   that such embedding  will also have shortest chains on average, i.e. chains of physical qubits representing a single QUBO variable.

The embeddings  shown in Figure \ref{fig:capitolo6:avg_num_embeddings} are highly representative of  the capabilities of each  QPU topology.
On average, the chains obtained with a Chimera topology are 2.09 times longer than those found using the  Pegasus topology, which therefore  returns  52\% shorter chains.

Regarding the number of qubits, Pegasus is again the winner as it enables embedding with 1.39 up to 2.4 times less qubits with respect to Chimera.
On average, Advantage is 1.96 times more efficient than 2000Q, or  51\% less expensive in terms of qubits; this will be reflected in the quality of results.
There are only two results slightly diverging from our expectations.
The  columns regarding the class $6-100$ are shorter than expected because, due to the complexity of the problem,   only  simpler instances were embedded. Hence the average is biased towards a lower value.

There other   important detail regards  the class  $8-100$ where  both  the average chain length and the average number of qubits are saturated to the maximum value for  the 2000Q Chimera topology. This happened because the device couldn't embed all the possible molecules for that class, therefore   radically increasing the average.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.45\textwidth]{./pictures/embeddings/avg_chains.png}
	\includegraphics[width=0.45\textwidth]{./pictures/embeddings/avg_qubits.png}
	\caption{Plots depicting the variation of  average chain length and  average number of qubits in the embeddings by increasing  number of rotatables and value of the approximating threshold (also called chop)}
	\label{fig:capitolo6:avg_num_embeddings}
\end{figure}

\section{Experimental Results}

A short remark should be made on the way these algorithms run.
Regarding the QA set-ups,  since the QUBO constant \emph{$A_{const}$}  modulating the strength of the hard constraint must be tuned, ten runs  for each value assigned to $A_{const}$ are performed. $A_{const}$ is evaluated as the maximum coefficient in the optimization contribution multiplied by an increasing factor. Each run is performing ten thousands forward anneals with an annealing schedule of one microsecond.
Simulated annealing is performed on the same QUBOs,   with the possibility of choosing the number of epochs and the function  which is responsible for the temperature   decrease. In our case we used 500 annealing epochs and a geometrical decrease function.

Since  executions of  quantum  and simulated annealing are repeated many times in the attempt of finding the best solution, also the other classical techniques are executed until their run-time does not exceed a  time limit which was set from the beginning.
Random search divides the solution space in N parallel processes and as long as it does not run out of time, it executes and records the maximal value found.
The Geodock Search is trying an incremental approach, which consists in performing the rotation on  an incremental number of bonds.

\subsection{Volume Gain in Time}

Plots in the Figures \ref{fig:capitolo6:distance_true} show the average trends in volume gain obtained by each algorithm on a fixed time window, increasing  rotatable bonds.
On the Y-axis,  the  maximum volume reached by every algorithm at each second is displayed.  The plots show the algorithmic behaviour up to 100 seconds since  this is the interval where  greatest changes in volume appear.

\begin{figure}[htpb]
	\includegraphics[width=.48\textwidth]{./pictures/dist_true/2_Rot_True.jpg}
	\includegraphics[width=.48\textwidth]{./pictures/dist_true/4_Rot_True.jpg}
	\includegraphics[width=.48\textwidth]{./pictures/dist_true/6_Rot_True.jpg}
	\includegraphics[width=.48\textwidth]{./pictures/dist_true/8_Rot_True.jpg}
	\caption{Comparison of average volume gained on a time window of 100 seconds for 2, 4, 6 and 8  rotatable bonds.}
	\label{fig:capitolo6:distance_true}
\end{figure}

First thing worth noticing is that  SA gave the best overall performance as it was able to provide good solutions falling inside the  $[90\%, 100\%)]$ volume gain range   with a run-time that never exceeded the 15 seconds.

It is also possible to observe that Advantage is the best among   D-Wave devices. For an increasing number of rotatable bonds,  volume achieved by 2000Q quickly dies off, while the results of Advantage remain competitive for a good number of torsionals.

When 2 torsionals are considered,  the MU problem looks like an easy task and  basically all the approaches are able to recover the optimal solution (only GeoDock gets close but fails to find optimal). Random search, SA and the Advantage QA are able to reach the  100\% volume gain within the first 5 seconds of runtime.

Increasing the problem size to 4 torsionals we can foresee a deviation of the trends: SA is again the fastest to reach a plateau while  the quantum annealers perform better than random search in the first 5 seconds. After that, random search presents a disruptive improvement which comes  with a delay due to  overheads in the memory access.
A promising result is obtained after 25 seconds, where Advantage is able to find the optimal solution among the methods. This  indicates that QA could be a  competitive approach when  solving  medium-small sized problems.

It is possible to see another  reversal of the trends when the number of torsions rises to 6, in which SA and random search establish themselves as the most reliable solvers, followed by the Advantage quantum annealer. This behaviour is confirmed at 8 rotatable bond, where performances of 2000Q become very similar to that of  the GeoDock greedy algorithm.

\subsection{Comparison between SA and  QA}

Comparing classical algorithms with quantum annealers in terms of absolute time isn't fair, since the quantum devices have run-times that are comprising of  programming time,  readout time and  delay time, besides the time spent in actually performing the annealing.

Therefore  more suitable metrics for comparing annealing solvers is given by the  Time To Solution (TTS),
\begin{equation}\label{eq:TTS}
TTS = \frac{total \ execution \ time}{occurrence \ of \ the \ best \ solution}
\end{equation}
The total execution time can be calculated as the number of anneals multiplied by  the total  access time for  single anneal.   The occurrences of the best solution simply count  the number of times the best solution is found in the annealing.
The TTS metric can be interpreted as the inverse of the probability of finding the optimal configuration in a unit time interval. Moreover, the TTS is telling us how long we have to wait on average before the annealer outputs the best solution. For this reason, better solvers have  low values of TTS.

\begin{figure}
	\centering
	\includegraphics[width=0.8\linewidth]{./pictures/velocityNstuff/tts.jpg}
	\caption{A comparison between Simulated Annealing, Advantage and 2000Q on TTSs, for increasing number of torsionals.}
	\label{fig:capitolo6:comparisons_tts}
\end{figure}

Plot in figure \ref{fig:capitolo6:comparisons_tts} show that the  TTS is very high for 2000Q, followed by Advantage, while SA has the the best TTS. Simulated annealing has a TTS several times smaller than Advantage, which means that its convergence to the optimum is much faster than in quantum annealing.

\begin{table}
	\centering

	\caption{Comparison on the TTS, for increasing number of torsionals.}
	\begin{tabular}{rrrr}
		\toprule
		torsionals &         sa &        adv &         two \\
		\midrule
		2 &   0.007171 &   0.165081 &    0.868552 \\
		4 &   4.138503 &  12.818871 &   65.482955 \\
		5 &   5.502274 &  34.508321 &   62.996198 \\
		6 &  10.384806 &  50.430784 &  204.457250 \\
		8 &  12.942243 &  43.363135 &  503.938826 \\
		\bottomrule
	\end{tabular}
	\label{tab:caption}
\end{table}%

Regarding the volumes achieved, plots in Figure~\ref{fig:capitolo6:comparisons_vol} represents the degradation in the final volume gain as the problem complexity increases. The volumes  have a  common starting point since all the three algorithm actually reach the maximum steadily and completely for problems involving 2 rotatable bonds.
As seen in the previous section, Advantage works well up to 4 torsions,  SA has an inflection on 4 bonds but then recovers and finishes as the best in absolute terms; it has a volume which is 4.65 times better than the one found by 2000Q at 8 bonds, and 1.4 times better than  Advantage.
\begin{figure}
	\centering
	\includegraphics[width=0.8\linewidth]{./pictures/velocityNstuff/volumi.jpg}
	\caption{A comparison between Simulated Annealing, Advantage and 2000Q on volumes.}
	\label{fig:capitolo6:comparisons_vol}
\end{figure}

\begin{table}
	\centering
	\caption{Comparison on the Volumes, for increasing number of torsionals.}

	\begin{tabular}{rrrr}
		\toprule
		torsionals &    sa &   two &   adv \\
		\midrule
		2 &  1.00 &  1.00 &  1.00 \\
		4 &  0.92 &  0.86 &  0.97 \\
		5 &  0.99 &  0.84 &  0.95 \\
		6 &  0.98 &  0.62 &  0.88 \\
		8 &  0.93 &  0.20 &  0.66 \\
		\bottomrule
	\end{tabular}
	\label{tab:caption}
\end{table}%
The last concept worth mentioning is the measure of how much volume in percentage  can be gained per  unit of TTS.
Practically, it is computed as the ratio between the normalized  volume gain and  TTS.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.4\textwidth]{./pictures/velocityNstuff/Normal_Vol_over_tts.jpg}
	\caption{A comparison between Simulated Annealing, Advantage and 2000Q on the normalised decreasing quality of the volumes per time to solution.}
	\label{fig:capitolo6:comparisons_velocity_normal_velocity}
\end{figure}

A negative slope in the plot of figure \ref{fig:capitolo6:comparisons_velocity_normal_velocity} can be interpreted  as an increasing TTS required to achieve the same volume when scaling up to a more difficult problem.
A higher negative slope means  faster worsening of the algorithm with the increasing difficulty, while a smaller slope is a synonym of stability. Therefore, trends   running towards the bottom of the plot  belong to  algorithms which are defective.
The SA line falls down very fast from 2 to 4 torsionals, and then starts becoming stable. The same happens to Advantage which slows down this negative acceleration before. 2000Q never stops decreasing its quality, if not for a special case between 4 and 5 torsionals.

By checking the ratio of the slopes when increasing the number of torsionals, we can understand how much the algorithms are worsening adding  degrees of freedom.
On average the smallest growth in scaling of the slopes belongs to Advantage, with 1.36, whereas SA presents a higher 1.609, meanwhile 2000Q has an average of 6.01.

Although the results obtained with  Advantage may not be the winners in absolute terms, the fact that Advantage has   good performances  with respect to SA according to this metric  is a very promising.  This result  reflects how much QA is still immature  but also how far from saturation the margin of growth for QA really is.   Improvements in the next generation of QPU hardware could definitely affect both TTS and quality of sampling in a positive way.

\begin{table}[htp]
	\centering
	\caption{Average Normalized Volume Gain over TTS slopes  for the three annealing algorithms. }
	\begin{tabular}{rrr}
		\toprule
		sa &       two &       adv  \\
		\midrule
		1.609948  & 6.0189816 &  1.36206524  \\

		\bottomrule
	\end{tabular}

	\label{tab:caption}
\end{table}%

\begin{table}[t!]
	\centering
	\caption{Normalized Volume Gain over TTS for each algorithm, at different number of rotatable bonds.}
	\begin{tabular}{rrrr}
		\toprule
		sa &       two &       adv &  torsionals \\
		\midrule
		2.001031 &  2.000690 &  2.005038 &   2 \\
		0.003190 &  0.022822 &  0.025046 &   4 \\
		0.002582 &  0.023171 &  0.009112 &   5 \\
		0.001354 &  0.005269 &  0.005776 &   6 \\
		0.001031 &  0.000690 &  0.005038 &   8 \\
		\bottomrule
	\end{tabular}

	\label{tab:caption}
\end{table}%

\section{Conclusions}

In this work, we explored the possibility of  using a quantum annealer to support the drug discovery process within the in-silico virtual screening phase.
In particular, we studied a new method regarding the process of molecular unfolding, which is the first step in geometric molecular docking techniques.
This phase aims at finding the molecule's configuration that maximizes its volume, or equivalently, that maximizes the internal distances of the atoms that compose the molecule.
We proposed our Quantum Molecular Unfolding model with the aim of executing it on D-Wave's latest hardware, Advantage and 2000Q.
The usage of the quantum annealers has the objective of understanding the capabilities of these quantum annealing devices and discovering if it is possible to improve the quality and the throughput of the ligand expansion in the state-of-art molecular docking methods.

The model is constructed starting from the identification of the rotatable bonds which are the parameters of the problem.  Once their discrete rotations are rewritten via  one-hot encoding thanks to the introduction of binary variables, the total sum of  internal atomic distances can be expressed as a HUBO. The HUBO is then transformed into a QUBO after   three steps are applied: (i) distance simplification, (ii) coarse-grained rotations  and (iii)  threshold approximation.
These three approximations enabled us to reduce the complexity of the model with the aim of  embedding and running complex instances on the  QPUs.

We compared the performances of Advantage and 2000Q, with three other optimization methods, which are parallel random optimization,  simulated annealing, and a GeoDock greedy algorithm. The outcomes of the experiments show that Advantage have good performances on medium-small  problems while  classical techniques are  better than  quantum approaches when problem size increases.   Interestingly, the results obtained by the quantum annealer, both while running on Advantage and 2000Q, are better than the GeoDock greedy approach.
In terms of the evolution of the quantum annealers, embedding our problems on Advantage, compared to 2000Q, was on average 51\% less expensive in terms of qubits and with chains 52\% shorter. Advantage significantly outperforms 2000Q in terms of time-to-solution (TTS) and volume gain when increasing the number of torsionals. Advantage also presents the best scaling of the normalized volume gain over time-to-solution when compared to both 2000Q and simulated annealing; this last result is one of the most promising as D-Wave devices are only bound to get better.

Future work can be the introduction of new dynamic thresholds approximation, which may lead to smaller embeddings but still of high quality, and  the introduction of   hybrid quantum approaches. A more difficult but also feasible and interesting task could be  extending  the quantum annealing approach to the whole docking process.
Finally, we are confident that this study will shed further light on the applications of quantum computing, thus enriching knowledge on topics that are becoming central in computational sciences.

\nocite{*}
\bibliographystyle{unsrt}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

\begin{figure}[!t]
\begin{center}
\includegraphics[width=250pt]{MvSigmav.ps}
\end{center}
\caption{ The distribution of Galactic GCs in absolute magnitude vs.
line--of--sight velocity dispersion (Pryor \& Meylan 1993). The inset
scale shows the limiting line--width (and thus spectral resolution,
${\sc R}=\lambda/\Delta\lambda$) obtainable for GCs as a function of
the line--of--sight velocity dispersion, $\sigma_v$. }
\label{fig:mvsigma}
\end{figure}

We are developing a method which will let us measure detailed chemical
abundances of Globular Clusters (GCs) from high--resolution, echelle
spectra of their integrated light. Our goal is to use the detailed
abundance patterns of GC systems in galaxies as distant as 4--5 Mpc to
constrain their formation and chemical enrichment histories
in the same way as the abundance patterns of old stars have
been used to do so in the Milky Way (\textit{e.g.}
Edvardsson \textit{et~al.}\ 1993;
McWilliam \textit{et~al.}\ 1995;
Prochaska \textit{et~al.}\ 2003;
Bensby \textit{et~al.}\ 2004ab;
Fulbright, McWilliam \& Rich 2005).

The elements which are the most useful for understanding the processes
of galaxy evolution are those which are returned to the ISM primarily
through either
high--mass, core--collapse SNe (Type II = SNII) or through
low--mass, accretion--induced SNe (Type Ia = SNIa). The key is
that the former (SNII) evolve on time--scales of order
$10^6$ years, while the latter (SNIa) generally take of order
$10^9$ years.  Elements produced by SNII ---
``$\alpha$--elements'' (\textit{e.g.} Si, S, Ca) and r--process elements
(\textit{e.g.} Eu) --- will therefore build up rapidly over the first few
megayears of star formation or in a starburst.  Elements produced
in SNIa and SNII --- Fe and Fe--peak elements (e.g. Sc, V, Cr, Mn, Fe,
Co, Ni) --- build up over many gigayears as the contribution from SNIa
increases.  The ratio [$\alpha$/Fe] relative to [Fe/H] is therefore a
rough constraint on the star formation rate.  Comparison with overall
abundance then constrains the timescale of formation.

Unfortunately, the Milky Way is currently the only galaxy for which
even a partial history of formation and enrichment can be traced in
the stellar fossil record through to the present day.  This is because
detailed abundance analysis requires high signal--to--noise ratios and
high spectral resolutions (S/N$\geq60$, $R>20,000$). Old stars are too
faint to be studied in this way beyond the nearest dwarf spheroidal
members of the Local Group ($\leq100$~kpc).  Even in the nearest
gas--rich irregulars of the Local Group, abundance measurements are
limited to a few very luminous supergiants (\textit{e.g.}\ Venn
\textit{et~al.}\ 2004), and results can only be obtained for a few
elements.  Moreover, these bright supergiants are very massive (young)
and so only identify the recent gas composition.  The history of a
galaxy can only be learned from long--lived, low--mass stars, which
are too faint for high--resolution spectroscopy beyond the nearest
dwarf spheroidal galaxies (\textit{e.g.}\ Shetrone \textit{et~al.}\
2003).

These limitations drive us to target Globular clusters (GCs).  Unlike
single stars, high--resolution spectra \emph{can} be obtained of
unresolved GCs out to $\sim4$~Mpc with \emph{current} telescopes.  GCs
are bright enough ($-10<M_v<-6$ mag) and have low enough velocity
dispersions ($2-20$ km/s) that even weak lines ($\sim15$ m\AA) could
be detected in spectra of their integrated light.  Moreover, several
lines of evidence from photometry and low--resolution spectroscopy
already suggest that GCs trace the star formation and the global
formation history of their parent galaxy: the number of GCs per galaxy
is roughly constant relative to the total light and mass of a galaxy;
young GCs are found in regions of active star formation
(\textit{e.g.}\ Schweizer \& Seitzer 1993, Barth
\textit{et~al.}\ 1995); and the properties of both metal--rich (red)
and now also metal--poor (blue) GC systems seem to correlate with the
properties of their parent galaxies (\textit{e.g.}\ Zepf \& Ashman
1993, Geisler, Lee \& Kim 1996; Gebhardt \& Kissler--Patig 1999; Kundu
\& Whitmore 2001; Larsen \textit{et~al.}\ 2001; Strader, Brodie, \&
Forbes 2004).  All of which implies that GCs are bright, observable
markers of the chemical enrichment record of normal galaxies.

Of course, significant progress has also been made in understanding
the chemical enrichment histories of elliptical galaxies and bulges
from analysis of low--resolution, integrated--light (IL) spectra using
line indexes, such as the Lick system (Burstein \textit{et~al.}\ 1984,
Faber \textit{et~al.}\ 1985, Worthey \textit{et~al.}\ 1994, Worthey \&
Ottaviani 1997, Trager \textit{et~al.}\ 1998).  These indexes provide
estimates of age and composite ``metallicity'' (``Z'') based on lines
from multiple elements such as Fe, Mg, and Ca (see Rose 1984, Worthey
\textit{et~al.}\ 1994, and Worthey \& Ottaviani 1997).  GCs are ideal
targets for these index systems as they are nearly ideal, single--age
stellar populations.

A great deal of progress has recently been made with line index
systems by using various kinds of principle component analysis to
obtain better leverage on Z and also to obtain an ``E'' parameter,
which includes all of the elements (O, Mg, C, etc) typically found to
be enhanced in elliptical and bulge (rapidly forming) populations
(Proctor \& Sansom 2002; Proctor, Forbes, Beasley 2004; Strader \&
Brodie 2004).  However, even for GCs, information from line indexes
are limited by two difficulties.  First, E and Z do not have the
interpretive power of detailed abundances because they include
elements which form in multiple sites and include elements which are
known to be affected by GC self--enrichment (see discussions in
Gratton, Sneden \& Carretta 2004). These indexes do not entirely
isolate the [$\alpha$/Fe] ratios which constrain formation timescales.
Second, line index results are very sensitive to the calibration of
the system, and it can be hard to obtain stellar libraries of the
appropriate age and enrichment or to know what calibration to use.
This complication has been widely discussed in the literature
(\textit{e.g.}\ Trippico \& Bell 1995; Trager \textit{et~al.}\ 2000ab;
Thomas, Maraston, \& Bender 2003; Proctor \textit{et~al.}\ 2004;
Tantalo \& Chiosi 2004).  Adjustments to the calibration system have
been published recently based on stellar models which include a range
of enrichments. However, not all of these agree with each other and
the models themselves are not well tested against observed
populations.  This concern is particularly relevant given recent work
on M31, which shows that the Milky Way GC system may not be generally
representative of GCs in even spiral galaxies; M31 appears to contain
a disk GC system (Morrison \textit{et~al.}\ 2004) which is at least partly
composed of young GCs (0.1--0.8 Gyrs; Barmby \textit{et~al.}\ 2000,
Beasley \textit{et~al.}\ 2004) with different abundance patterns than
are seen in the Milky Way (Burstein \textit{et~al.}\ 2004).

Fortunately, it \emph{is} possible to obtain high--resolution spectra
of unresolved GC.  A typical GC has velocity dispersions of
$2<\sigma_v<20$ km/s (Pryor \& Meylan 1993, see Figure 1), so that
spectra of their integrated light can have line widths in the range
$60,000>\lambda/\Delta\lambda>6,500$.  For comparison, an elliptical
has $\sigma_v\sim 200$ km/s or greater and limiting line widths of
$\lambda/\Delta\lambda \sim 1500$.  Absorption features as weak as a
$\sim15$ m\AA\ can therefore be seen in high--resolution, {\it integrated
 light} spectra of all but the most massive GCs  with S/N$\approx
60-90$.

Detailed abundances have never been obtained for unresolved GCs
because current methods of abundance analysis only work for individual
stars.  We describe here progress that we have made in adapting the
basic techniques of single star abundance analysis to spectra of
integrated--light.

\section{A Method for Light--Weighted Abundance Analysis}
\label{sec:method}

Abundances are not directly observable quantities; the strength of any
given absorption line is a function not only of the abundance of the
element, but also of the physical properties (e.g. mass and
temperature) of the star.  The method we are developing for analyzing
IL spectra is based on the standard techniques used to analyze
individual stars.

The standard technique for abundance analysis of single stars involves
comparing the observed equivalent widths (EWs) of lines for each
species with those predicted by model stellar atmospheres and spectral
line synthesis.  The model atmospheres have the following, observationally
constrained parameters: effective temperature ($T_{\rm eff}$, from
$B\!-\!V$), specific gravity ($\log g$, from luminosity), and
microturbulence ($\xi$).  In principle, the only free parameter is
abundance, [A/H].  So far, we have used Kurucz models in which Fe alone
is varied, so that [A/H] is roughly [Fe/H].  The Kurucz model
atmosphere grids (available from R. Kurucz at {\tt
http://cfaku5.harvard.edu}) give optical depth, temperature, pressure,
and electron density in each of 64 layers.  The EW of any particular
spectral line is then found by radiative transfer through these
layers. Given the model atmosphere, the EW of a line is uniquely
determined by the wavelength ($\lambda$), excitation potential (EP),
and {\it gf} value of the line, and by the abundance of the element in
question, $[X/H]$.  Of these, only [X/H] is adjusted.  In standard
analysis, the best model atmosphere is found by adjusting [Fe/H] in
the line synthesis and model atmospheres until a unique,
self--consistent solution is found from all Fe I and Fe II lines
($\sim 100$ lines).  A strong constraint on the correct model
atmosphere comes from requiring that Fe I lines over the full range of
$\lambda$, EP, and EW give a stable Fe solution.  This requirement is
used to adjust $T_{\rm eff}$, $\log g$, and $\xi$  to find the
appropriate stellar atmosphere model.  Lines from any other element
are then uniquely determined by the input abundance of that element
([X/H]) in the line synthesis step.  The right abundance is found by
adjusting this value until predicted EWs match those observed.

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=250pt]{6397cmd.ps}
\end{center}
\caption{A color--magnitude diagram for NGC~6397 (Kaluzny 1997).  The
scale on the left shows the cumulative $V$--band flux coming from
bright to faint stars.  The boxes sub--divide the stellar population
into groups containing $\sim5$\% of the $V$--band light, as indicated.
An analysis based on this CMD, for example, would
include one stellar atmosphere model per box.}
\label{fig:6397cmd}
\end{figure}

To interpret an IL spectrum, we have developed an original method of
detailed chemical abundance analysis which involves producing a
light--weighted EW, building on those techniques used for individual
stars. For a {\it resolved} Galactic GC, the observed CMD tells us the
exact fraction of light coming from every stellar type in the
cluster. For NGC~6397, for instance, Figure 2 shows the CMD of
NGC~6397 with boxes indicating groups of stars which each contain
roughly 5\% of the cluster's light.  Based on this observed CMD, we
can compute a model atmosphere for each box using the Kurucz model
atmospheres (ATLAS9) for a range in abundance, [A/H].  We have
developed code which then synthesizes the individual Fe lines by
iteratively calling the spectrum synthesis program MOOG (Sneden 1974)
and varying the Fe in the line synthesis only until consistent
solutions are obtained.  With this Fe abundance, we can then compute a
final stellar atmosphere for each box. We then identify the abundances
of other elements by adjusting [X/H] (in the line synthesis only) to
match the observed EWs for a species.  In this analysis, a
microturbulence scaling law is adopted ($\xi \propto \log gf$).

\begin{figure}[!t]
\begin{center}
\includegraphics[width=130pt]{ep9.00_004.eps}\hspace{-0.13in}
\includegraphics[width=130pt]{ep9.30_004.eps}\hspace{-0.13in}
\includegraphics[width=130pt]{ep10.00_004.eps}
\end{center}
\caption{Excitation potential vs [Fe/H] for all iron lines measured in
the IL of 47 Tuc.  Each plot shows the best--fit [Fe/H] solution for
one isochrone.  A systematic negative trend is clear in the 1 Gyr
isochrone, but weak in older isochrones.  Scatter is largely due to
uncertain \emph{gf} values.  It should be possible to obtain greater
precision by using \emph{gf} values tuned to the abundance solution of
Arcturus, a nearby metallicity reference star.  In RGB star analysis,
such techniques reduce the scatter in EP vs.\ [Fe/H] to $\pm0.05$ {\it
rms}.}
\label{fig:fevsep}
\end{figure}

As a first test of this basic strategy, we have analyzed the
integrated light spectrum of NGC~6397 using its observed CMD (see
Figure 2) to define the contributing stellar atmosphere models and
their weights.  From this analysis, we successfully obtain abundances
for Fe--peak, $\alpha$--elements, and light elements which are in good
agreement with the results from other groups for single stars in
NGC~6397 (see Bernstein \& McWilliam 2001 and Bernstein \& McWilliam
2005).  However, when we target extragalactic GCs, we will not be
able to use resolved CMDs to empirically identify the right mix of
stellar populations.  We must therefore learn to use theoretical CMDs
--- isochrones based on stellar evolution tracks --- to analyze the
IL spectra. To do this, and to understand the limitations and impact
of stellar evolution models themselves, we have observed a set of
GCs in the Milky Way and LMC to serve as a ``training set''.

\section{A Training Set of Galactic and LMC GCs}
\label{sec:tset}

Our ``training set'' consists of seven Milky Way and eight LMC
GCs. These span the range of velocity dispersions, abundances,
ages, and HB morphologies available in the Milky Way and LMC systems
(see Figure \ref{fig:mvsigma} and Table \ref{table:tset}).  For each
of the training set GCs, we have obtained integrated light spectra by
uniformly scanning the central $32\times32$ arcsec$^2$ and
$12\times12$ arcsec$^2$ in the Milky Way and LMC GCs, respectively.
Crucially, these clusters are all spatially resolved, so detailed
color magnitude diagrams (CMDs) can provide age estimates and \emph{a
priori} knowledge of the member stellar populations, including HB
morphology.  We can use this information to refine our methods and
determine how such variables affect our results.  Moreover, we can
also identify individual red giant branch (RGB) stars from these CMDs
and use spectra of these individual stars to obtain ``fiducial''
abundances by standard methods which are consistent with our own
atmospheric modeling and synthesis for the integrated light.  Below
we describe the results we have obtained for one Milky
Way GC: NGC~104 (47~Tuc).

To analyze 47 Tuc as we would an unresolved GC, we use isochrones as
the template for its stellar population.  The analysis here uses
isochrones based on the stellar evolution tracks of the Padova group
(Girardi \textit{et~al.}\ 2000) and a Kroupa IMF (Kroupa \& Boily
2002).  We complete the same analysis for isochrones with a range of
age and abundance (1--16 Gyrs, $-2.3<$[Fe/H]$<0.2$).  The first
question one might ask is whether the resulting $\sim15$ Gyr,
[Fe/H]$\sim -0.7$ isochrone \emph{looks} like the observed CMD of
47~Tuc.  Two differences are evident, although not unexpected
(\textit{e.g.}\ Bergbusch \& VandenBerg 2001): the isochrones include
lower mass stars than are observed in the clusters and do not
reproduce the observed AGB bump. The same problems exist with the
BASTI (Pietrinferni \textit{et~al.}\ 2004) stellar evolution models.
To begin, we use the Padova isochrones ``as is,'' with no adjustments
for mass segregation or the AGB bump.

\begin{figure}[!t]
\begin{center}
  \includegraphics[width=375pt]{FeIvsFeII.cps}
\end{center}
\caption{ The inferred [Fe/H] abundance from Fe I and Fe II lines for
Padova isochrones of age 1--16 Gyrs and abundances of [A/H]$=-2.35$ to
0.2.  These isochrones do not include $\alpha$ enhancement.  The Fe I
lines give stable solutions with $\lambda$ and EP for isochrones with
ages older than about 3 Gyrs (see Figure \ref{fig:fevsep}).  The Fe II
and Fe I lines give the same solution (the lines connecting Fe I and
Fe II solutions cross) at a narrow range in abundances, namely at
[Fe/H] = $-0.63 \pm 0.05$.}
\label{fig:FeIFeII}
\end{figure}

For each isochrone, we split the stellar population into roughly 20
groups each containing roughly 5\% of the light, similar to the boxes
shown in Figure \ref{fig:6397cmd}.  We then produce a stellar
atmosphere for each group and synthesize a light--weighted EW for
about 70 Fe I and 10 Fe II lines by combining the synthesized
EWs from MOOG for each group, as described in \S\ref{sec:method}
Note that with the isochrones, the parameters of the atmosphere model
are completely dictated by the parameters of the isochrone.  We adjust
only the abundance of Fe in the line synthesis step to obtain a good
match between observed and predicted EWs for all Fe I and Fe II
line. Figure \ref{fig:fevsep} shows the best--match Fe value for
Fe I lines as a function of excitation potential for three
different isochrones.  The average of these solutions gives us the
inferred [Fe/H] value (and a statistical uncertainty) for each
isochrone.  Preferred solutions have small {\it rms} scatter and no slope
with excitation potential, wavelength, and EW.  We then do the same
test for Fe II lines.  A plot showing the inferred value of [Fe/H] from
Fe I and Fe II lines for each isochrone is shown in Figure
\ref{fig:FeIFeII}.  Surprisingly, we find that best--fit [Fe/H]
solution from Fe I lines is nearly independent of the input isochrone
{\it metallicity}, while Fe II lines are nearly independent of the
input {\it age}.  Together, the Fe I and Fe II lines provide an
excellent constraint on the abundance of the cluster and a way to
break the age--metallicity degeneracy that will never be achievable at
low resolution.

From the best--fit isochrone, we get [Fe/H]=$-0.63\pm0.03$ ($1\sigma$
statistical uncertainty).  We then adjust the luminosity function of
one isochrone to remove low--mass stars and add an AGB bump.  With
this adjusted isochrone, we find [Fe/H]$=-0.73$.  Note that these
results differ by only 0.1 dex, and that both results are roughly
within the range of values obtained recently by different groups from
individual stars (Kraft \& Ivans 2003, Carretta \textit{et~al.}\
2004).  Nevertheless, the change in our [Fe/H] solution with these
adjustments emphasizes the importance of the luminosity function to
our results.  We will use the rest of the training set to determine the
best approach to dealing with these differences between isochrones
and real CMDs.

We use the adjusted isochrone to then obtain abundances for a broad
range of key elements (see Table \ref{table:47tuc}).  Our final
abundances are in good agreement with those in the literature for
individual 47 Tuc stars.  Moreover, we find the expected abundances
patterns $\alpha$-- (enriched relative to solar), Fe--peak (roughly
solar), and r--process (enriched) elements that one would expect for
an old globular cluster.  We also detect enrichment of Na and Al,
which is a consistent with self--enrichment through proton--burning in
AGB stars and is found for some individual stars (with star--to--star
variations) in 47 Tuc (Carretta \textit{et~al.}\ 2004). This suggests
that we will be able to measure such abundance variations in
unresolved clusters even if present in only a fraction of the stars in
a given cluster!  Finally, we note that the age of 47 Tuc is only
constrained from our analysis to be greater than $\sim3$ Gyrs. We may
be able to improve this by calibrating the $\log gf$ values of
individual Fe lines to the abundance solution for Arcturus; doing
so should reduce the scatter in plots like those shown in Figure
\ref{fig:fevsep} to $\pm 0.05$, as it does in the analysis of
individual stars (see Bernstein \& McWilliam 2005).  However, it
is also clear from evolutionary tracks and observed CMDs of Milky Way
GCs that stellar populations themselves simply do not vary much with
age after a few Gyrs (\textit{e.g.}\ Yi \textit{et~al.}\ 2003).

\section{Further Spectroscopic Constraints on Integrated Light Abundances}

An important issue to explore is the influence of Horizontal Branch
(HB) morphology on our abundance analysis because it is not uniquely
correlated with age or metallicity and therefore not reliably
characterized by the isochrones.  Balmer line EWs \emph{and profiles}
will be helpful in this regard, because they are very sensitive to the
light fraction contributed by hot stars (\textit{e.g.}\ Bernstein \&
McWilliam 2002; Schiavon \textit{et~al.}\ 2004; Bernstein \& McWilliam
2005;).  To this end, we have compared observed Balmer profiles for
several of our training set clusters with synthesized, light--weighted
Balmer profiles based on the isochrones. Our preliminary work suggests
that we will need to also synthesize blended lines in the Balmer line
wings in order to get accurate Balmer EWs and profiles, particularly
in the metal rich clusters. This may be a particular problem for
low--resolution spectra. Similar issues have already been noted in the
literature to the extent that Balmer line indexes in low--resolution
spectra seem to be sensitive to the calibration of the line system
(\textit{e.g.}\ Proctor, Forbes \& Beasley 2004).  Additional
constraints may also come from CaII lines, which have also been shown
to track HB color in low--resolution spectra (e.g. Burstein
\textit{et~al.}\ 2004; Proctor, Forbes \& Beasley 2004).  Independent
models for horizontal branch morphology can also be combined with the
isochrone models directly to explore the impact on inferred abundances

We are also planning to explore a broader range of parameters in the
isochrones themselves as more models become publicly available.  For
example, stellar evolution tracks with different abundance ratios
(e.g. enriched in $\alpha$--elements or with different helium
fractions) may yield systematically different abundances and should
also be tested against the training set.

\section{Applications}

It is already possible to obtain spectra of extragalactic GCs that can
be analyzed with this technique.  In the local group, abundances of
the more distant LMC GCs could be obtained in a few hours with 4--m
class telescopes taking IL spectra, instead of requiring many hours of
exposure time on 8--m class telescopes observing individual cluster
stars.  Beyond the LMC, an handful of galaxies within 4--5~Mpc can be
observed using using current 6.5--10~m telescopes.  We have already
obtained a few spectra ($R\sim 20,000$, S/N$=60-900$ at $H\alpha$) of
confirmed GCs in galaxies NGC~1313 and NGC~5128 (both $\sim4$~Mpc
away) with the MIKE echelle spectrograph (Bernstein
\textit{et~al.}\ 2003) on the Magellan Telescopes.  In the future,
hundreds of galaxies in the local universe could be observed with the
next generation of ground--based 20--30~m telescopes.  In addition,
the abundances and integrated light spectra we already have in hand
can be used to better understand and calibrate the low--resolution
data and line--index systems.  This is important because line indexes
will always be able to reach galaxies at greater distances than
the high--resolution technique described here.

\acknowledgements We gratefully acknowledge the help of Las Campanas
technical staff and particularly telescope operators Fernando Peralta
and Herman Olivares. R.~A.~B.\ acknowledges partial support from NASA
through Hubble Fellowship grant HF--01088.01--97A awarded by Space
Telescope Science Institute, which is operated for NASA by the
Association of Universities for Research in Astronomy, Inc.\ under
contract NAS 5--2655.  A.~M.\ acknowledges support from NSF grants
AST--96--18623 and AST--00--98612.  We also thank the organizers for a
very enjoyable and stimulating conference.

\setcounter{table}{0}

\begin{landscape}
\begin{table*}[!h]
\renewcommand{\arraystretch}{1.0}
\begin{center}
\caption[]{\footnotesize
`Training set'' clusters in the Milky Way and LMC.}
\begin{tabular}{l c c c c c c c c c c }
\scriptsize
\hfil \\
\hline
\hline
GC      & RA    &  Dec  & $\mu_V$  & $V_{\rm tot}$ & $M_{V{\rm tot}}$ & HBR$^a$   & $ R_\odot$ & ${r_{\rm core}}^b$ & $\sigma_v{\small\rm core}$  &  [Fe/H]\\
        &{(J2000)}&(J2000) & ${{\rm mag}/{\rm asec}^{2}}$ & mag & mag      & & kpc        & pc              & km/s        &      \\
\hline
\multicolumn {6}{l}{\underline{Milky Way clusters}} &  & \\
NGC 104	&  00 24 05.2  	&  --72 04 51  	&  14.42  & 3.95  & --9.42   & --0.99    &   4.5    & 0.40    &  9.8 & --0.65 \\
NGC 362	&  01 03 14.3  	&  --70 50 54 	&  14.79  & 6.40  & --8.40   & --0.87    &   8.5    & 0.17c   &  6.3 & --1.16 \\
NGC 2808&  09 12 02.6   & --64 51 47     &  15.17  & 6.20  & --9.39   & --0.49    &   9.6    & 0.26    & 14.2 & --1.15 \\
NGC 6093&  16 17 02.5	&  --22 58 30	&  15.38  & 7.33  & --8.23   & 0.93     &  10.0	  & 0.15    & 12.5 & --1.75 \\
NGC 6388&  17 36 17.0  	&  --44 44 06  	&  14.50  & 6.72  & --9.82   &          &  11.5    & 0.12    & 18.9 & --0.60 \\
NGC 6397&  17 40 41.3  	&  --53 40 25  	&  16.66  & 5.73  & --6.63   & 0.98     &   2.3    & 0.05c   &  3.3 & --1.95 \\
NGC 6752&  19 10 51.8 	& --59 58 55 	&  15.20  & 5.4   & --7.73   & 1.00     &   4.0    & 0.17c   &  4.5 & --1.83 \\
\hfil \\
\multicolumn {7}{l}{LMC clusters} & log(Age) & \\
\hline
\multicolumn{7}{l}{{\underline{Old ($>5$Gyrs)}}}  \\
NGC 1916 	& 5 18 39.00   & --69 24 24     & 15.8   & 9.88  & --8.96  & & $>$9	& $<1.5$& $8.2\pm1.2$ 	&--2.08   \\
NGC 2019 	& 5 31 56.7    & --60 09 33     & 15.5  	& 10.7 	& --7.94  & & $>$9	& 0.9	& $7.5 \pm 1.3$	&--1.3   \\
NGC 2005 	& 5 30 10.36   & --69 45 09     & 15.8 	& 11.57 & --7.48	 & & $>$9	& $1.3$ & $8.1 \pm 1.3$	&--1.92   \\
\multicolumn{9}{l} {{\underline{Intermediate age (0.1--1.5G yrs)}}} \\
NGC 1866 	&  5 13 39     & --65 27 54     & 18.5   & 9.76  & --8.74	 & & 8.1	&11.7   & ...		& --0.51 \\%\pm0.08^c$ \\
NGC 1978 	&  5 28 45.00  & --66 14 12     &        & 10.7	& --7.7	 & & 9.3	& ...   & $3.0\pm 0.5$ 	& $-0.91$\\%\pm0.15^c$  \\
\multicolumn{9}{l} {{\underline{Young ($\sim10$ Myr)}}} \\
NGC 1711  	& 4 50 37     & --69 59 06     & 17.0    & 10.11 & --8.3	 & & 7.4	& 5.7	&  ...	        & $-0.6$  \\
NGC 2002	& 5 30 21     & --66 53 00     & 16.0	& 10.1	& --8.3 	 & & 7.2 	&3.5	& ... 	        &  ...    \\
NGC 2100 	& 5 42 08     & --69 12 42     & 17.0    & 9.6   & --8.8   & & 7.2	&6.7	& ...	        & ...     \\
\hline
\hline
\label{table:tset}
\end{tabular}
\end{center}
\noindent{\footnotesize Notes: (a) HBR, the horizontal--branch ratio,
  is defined as $(B-R)/(B+V+R)$.  (b) A ``c'' indicates that the
  cluster is core--collapsed.  For the Milky Way GCs, values for 1--D
  velocity dispersion, $\sigma_v$, are from Pryor \& Meylan (1993);
  all other values are from Harris (1996).  Data for the LMC clusters
  comes from various sources, including Johnson \textit{et~al.}\ 2004, Hill
  (2004), Olsen \textit{et~al.}\ 1998, Olszewski \textit{et~al.}\ 1996
  and references therein.}
\end{table*}
\end{landscape}

\begin{table}[!ht]
\renewcommand{\arraystretch}{0.7}
\begin{center}
\caption[]{\footnotesize Element Abundances for 47 Tuc from
Integrated Light Spectra}
\begin{tabular}{l c c c c c| c }
\scriptsize
\hfil\\
\tableline
\tableline
& & & & & &  \\
{ Species, X} &
{ $\epsilon(X)\,^{a}$} &
{ $\sigma$} &
{ $N_{\rm lines}$} &
{ $(\sigma/\sqrt{N})^b$} &
{ [$X$/Fe]$\,^c$ }&
{ [$X$/Fe]$\,^d$ } \\
& & & & & &  \\
\tableline
& & & & & &  \\
\multicolumn{6}{l} {\underline{Fe--peak}}\\
Sc II	& 2.53	& ...  	& 1	& ...  & +0.13	& +0.13   \\
V  I	& 3.40	& 0.31 	& 5	& 0.16 & +0.11	& +0.05\\
Cr I 	& 4.87	& 0.14	& 3 	& 0.10 & --0.19	& +0.11\\
Mn I	& 4.49	& 0.30	& 4	& 0.17 & --0.31	& --0.29\\
Fe I	& 6.78 	& 0.24	& 71	& 0.03 & --0.73	& --0.67, --0.63$^f$\\
Fe II	& 6.81	& 0.15	& 8	& 0.06 & --0.70	& --0.56, --0.58$^f$\\
Ni I	& 5.47  & 0.18	& 12	& 0.05 & --0.05	& +0.06\\
\multicolumn{6}{l} {\underline{$\alpha$--elements}}\\
Si I	& 7.16 	& 0.21 	& 6     & 0.09 	& +0.33 & +0.30 \\
Ca I	& 5.81	& 0.24	& 12	& 0.07	& +0.19 & +0.20 \\
Ti I$^g$& 4.55	& 0.24 	& 13	& 0.07	& +0.34 & +0.26 \\
Ti II$^g$& 4.68 & 0.16	& 3 	& 0.11 	& +0.44 & +0.38 \\
& & & & & &  \\
\multicolumn{6}{l} {\underline{Light elements with possible intra--cluster variations}}\\
Mg I	& 7.02	& ... 	& 1	& ... 	& +0.17 & +0.40 \\
Al I	& 6.19  & 0.02  & 2	& 0.02	& +0.43 & NA    \\
& & & & & &  \\
\multicolumn{6}{l} {\underline{Neutron--capture elements  (s--, r--process)}}\\
Y II	& 1.34 	& ...  	& 1 	& ...	& --0.19  & +0.49$^e$\\
Y I:	& 1.29  & ...  	& 1 	& ...	& --0.24  &  NA\\
Zr I	& 1.80  & 0.21	& 2 	& 0.21  & --0.08	 & --0.22$^e$\\
Ba II	& 1.41  & 0.09	& 3	& 0.05  & --0.11  & NA	\\
La II	& 0.57  & 0.38	& 2	& 0.26  & +0.05  & NA\\
Nd II	& 0.80	& ...  	& 1	& ...	& +0.01  & NA\\
Eu II	& --0.12 & ...  	& 1 	& ...	& +0.04  & +0.36$^e$\\
\tableline
\tableline
\label{table:47tuc}
\end{tabular}
\end{center}
\noindent{\footnotesize{Notes: (a) The abundance of
element $X$ is defined as  $\epsilon(X) = \log_{10}[{\rm N(X)/N(H)}] +12$.
(b) Realistic uncertainties probably lie between the {\it rms} error,
$\sigma$,  and the error in the mean, $\sigma/\sqrt{N}$.
(c) $[X/{\rm Fe}] = \epsilon(X/{\rm Fe}) -\epsilon(X/{\rm Fe})_\odot $.
For Fe I  and Fe II this column indicates [Fe/H].
(d) From Carretta \textit{et~al.}\ (2004).
(e) From Brown \& Wallerstein (1989).
(f) From  Kraft \& Ivans (2003) using the Kurucz models;
also, values of --0.70 and --0.65 were obtained from Fe I and Fe II,
respectively, using the Gustafsson \textit{et~al.}\ (1975) models.
(g) Ti can also be considered an Fe--peak element and does
show different abundance ratios than Si in Galactic GCs
(Lee \& Carney 2002).}}
\end{table}