text:
# MAVID: Constrained ancestral alignment of multiple sequences

## 1 Introduction

The multiple alignment problem is difficult for many reasons (see, for example, \[Notredame 2002\]), and thus remains unsolved despite much progress over the past two decades. To appreciate the complexity of the problem, it is instructive to observe that five DNA sequences of length five have approximately $`1.0510^{18}`$ alignments! \[Slowinski 1998\] Even three sequences of length five have over 14 billion alignments, and thus it is clear that an alignment of the human, mouse and rat genomes must be based on a relatively simple optimization criteria, and even then must involve heuristics to reduce the complexity of the problem.

Despite the overwhelming complexity of the problem, it is important to observe that the multiple alignment problem for genomic sequences is not equivalent to the mathematical problem of producing an optimal alignment maximizing some score function \[Gusfield 1997\]. The biological problem consists of correctly aligning homologous bases to each other, thus correctly identifying conserved non-coding regions in introns and intergenic regions, exons in orthologous genes, and groups of orthologous genes that form larger blocks of homology.

In this paper, we propose a method capable of rapidly aligning multiple large genomic regions by incorporating biologically meaningful heursitics, with theoretically sound alignment strategies. The core of our approach is a probabilistic ancestral alignment scheme \[Feng 1987, Gonnet 1996, Holmes and Bruno 2003, Hein 2001, Löytynoja 2003\]. This involves the progressive alignment of ancestor sequences (inferred using maximum-likelihood estimation within a probabilistic evolutionary model \[Felsenstein 1981\]) along a phylogenetic guide tree. Although a comprehensive review of progressive alignment is beyond the scope of this paper, it is important to point out that probabilistic approaches have been proposed and implemented (e.g. \[Holmes and Bruno 2003, Hein 2001\]), although existing methods are not scalable to very large problems.

In order to incorporate biological information into the alignment procedure, the progressive alignment is constrained by gene-based anchors. These anchors are pre-computed based on ab-initio gene predictions and their protein alignments, and form part of the input to the program. In addition, non-trivial positional constraints \[Myers 1997, Hardison 1993\] are pre-computed and ensure that the progressive alignment steps respect a pre-computed homology map for the sequences.

The alignment of the ancestor sequences is based on the AVID \[Bray 2003B\] alignment method, thus allowing for the rapid alignment of very large genomic sequences, even in between gene anchors that may be far apart. This fast alignment, along with the speedup obtained by using constraints, allows for an iterative alignment approach alternating between the progressive alignment step and phylogenetic tree construction (based on the alignment). In fact, as we show, it is possible to start with a random initial tree and converge to the correct guide tree thus eliminating the need for an expensive pairwise alignment step (quadratic in the number of sequences) at the beginning of the progressive alignment \[Thompson 1994\].

We have combined all these ideas into a new program called MAVID, which we used to align the human, mouse and rat genomes. We also show that MAVID is suitable for aligning very large numbers of sequences, and is therefore practical for the alignment of multiple HIV genomes \[Korber 2001\] or hundreds of mitochondrial sequences \[Hernnstadt 2002\]. Finally, we demonstrate the accuracy of MAVID on the benchmark cystic fibrosis (CFTR) gene region \[Green 2003\] and show that it compares favorably to existing alignment methods.

## 2 Method

Our method consists of a core progressive ancestral alignment step, which can incorporate pre-processed constraints (see Figure 1).

In order to clarify the presentation of the method, we begin by defining some terminology: A match is any identified similar region between two sequences. A match does not have to be an exact at the sequence level; for example, we could declare a match between two orthologous gene regions even if the sequence does not match exactly. A maximal exact match is a match that is exact at the sequence level, and is maximal (i.e. cannot be extended on either side without creating a mismatch). An anchor is a match that is used in the alignment. A constraint $`a_ib_j`$ in a multiple alignment means that position $`i`$ in sequence $`a`$ must appear before position $`j`$ in sequence $`b`$ in the multiple alignment.

### 2.1 Progressive Alignment

Let $`T`$ be a phylogenetic tree, i.e. a binary tree with branch lengths from an evolutionary model. In order to build progressive alignments we need a root for $`T`$. Sometimes the root is known because of the availability of an outlier. In the case in which the tree $`T`$ is unrooted, we use the midpoint method: select the root to be halfway in between the two leaves of $`T`$ which are farthest apart. This method produces the correct root if the molecular clock assumption holds, and a good approximation otherwise.

We associate multiple alignments to the vertices of $`T`$ recursively, starting from the leaves. For a vertex $`x`$ in $`T`$, let $`T_x`$ be the subtree consisting of $`x`$ and all the vertices beneath $`x`$. If $`x`$ is a leaf then $`T_x`$ is a “trivial tree”(i.e. a tree consisting of only one vertex) and we will label $`x`$ with the sequence corresponding to $`x`$ in the phylogenetic tree. This sequence can be considered a “trivial multiple alignment”(i.e. an alignment of only one sequence). If $`x`$ is an internal node then it has two children, $`u`$ and $`v`$, which are labeled with the multiple alignments $`A_u`$ and $`A_v`$, respectively. We then construct an alignment of all the sequences in $`T_x`$ by aligning the two alignments $`A_u`$ and $`A_v`$. This procedure is applied recursively so the program works its way from the leaves of the phylogenetic tree to the root, at which point it will have constructed an alignment for all the sequences in $`T`$.

The key difference between our progressive alignment schema, and more standard methods, is that instead of aligning $`A_u`$ and $`A_v`$ directly, we first infer ancestral sequences $`s_u`$ and $`s_v`$ using standard phylogenetic models for inference of the common ancestor \[Felsenstein 1981\]. We used the general reversible model, with the rate matrices obtained by \[Yap and Pachter\].

The discussion above ignores the issue of gaps. Gaps can be modeled as a fifth symbol which is equivalent to assigning a linear gap penalty. We have implemented the procedure in this form, but affine gap penalties are preferable. Furthermore, it is desirable to infer the deletion or insertion of bases in the ancestor, and models for this already exist (e.g. TKF \[Thorne 1992\]). For the human-mouse-rat alignment the issue of properly scoring gaps while inferring the ancestral sequence was not critical, and so we did not score with a sophisticated model, however future work will build on probabilistic insertion/deletion models developed in \[Thorne 1991, Thorne 1992\] and which have already been used to develop multiple alignment algorithms \[Holmes and Bruno 2001\]. It is important to note that in the MAVID alignment scheme gaps also play a role in Smith-Waterman alignments of the ancestral sequences (see below).

After the ancestral sequence calculation, $`s_u`$ and $`s_v`$ are aligned with AVID \[Bray 2003B\]. AVID is a hierarchical global pairwise alignment program that iteratively anchors exact matches and wobble matches (i.e. matches which are exact except for possible mismatches every third base) until a final Smith-Waterman alignment step of remaining regions. In the Smith-Waterman phase of AVID, the match and mismatch scores are again assigned according to a substitution matrix corresponding to the branch length between $`u`$ and $`v`$ (using the same rate matrix as for the ancestral inference). Gap scores were assigned using the AVID protocol, however both the gap open and gap extension scores were scaled according to the evolutionary distance.

The alignment of the ancestral sequences is then used to ”glue” together $`A_u`$ and $`A_v`$ to produce a new multiple alignment which is assigned to the vertex $`x`$. In particular, if position $`i`$ in $`s_u`$ matches position $`j`$ in $`s_v`$, then column $`i`$ in the multiple alignment assigned to $`u`$ is aligned with column $`j`$ in the multiple alignment assigned to $`v`$. Gaps in the ancestral sequence alignments lead to gaps in the multiple alignment in the obvious way. The procedure terminates with a final pairwise alignment at the root node of the tree.

### 2.2 Exon Anchoring and Constraints

Our gene matches and constraints are based on a homology map for the sequences; this is a map that identifies the order and orientation of matching gene runs between the sequences (Colin Dewey, manuscript in preparation). First, pairwise gene matches are computed between all the sequences. Gene predictions are generated using GENSCAN \[Burge and Karlin 1997\], and every pair of predicted genes from every pair of sequences is aligned using the translated BLAT tool \[Kent 2002\]. GENSCAN was selected because it is sensitive, and BLAT provides a fast way of obtaining protein alignments between large numbers of sequences. It should be noted that these programs can easily be replaced for different types of organisms, for example a viral gene finding program is more suitable for virus alignments.

Genes are considered to match if they form a reciprocal best hit. The genes matches are assembled into runs which then form the basis of the homology map. Genome sequence coordinates for exon matches are inferred from the protein alignments, thus producing a set of pairwise matches between all the sequences. These matches are used in the obvious way: when aligning $`A_u`$ and $`A_v`$ at node $`x`$, all matches that are between a sequence in $`u`$ and a sequence in $`v`$ are collected. Every such match can be converted into a match between the ancestral sequences $`s_u`$ and $`s_v`$ which is then used in the AVID alignment.

In addition to anchoring the alignment of the ancestral sequences, the exon matches can be used in more subtle ways to shape the final multiple alignment. It is illustrative to consider one base-pair anchors, i.e. single matches between the sequences. Suppose we have sequences $`a`$, $`b`$, and $`c`$, and that $`a_i`$ is anchored to $`c_x`$ and $`b_j`$ is anchored to $`c_y`$. If we are aligning sequences $`a`$ and $`b`$ then the given anchors to $`c`$ do not allow us to anchor the alignment, but they do allow us to constrain it. If $`x`$ is less than (resp. greater than) $`y`$, we must have that $`a_i`$ comes before (resp. after) $`b_j`$ in the alignment of $`a`$ and $`b`$ if we are to produce an alignment of $`a`$, $`b`$, and $`c`$ which is consistent with both of the anchors. In the language of Miller and Myers \[Myers 1997\], the two anchors provide explicit constraints on the alignment (namely that $`a_ic_x,a_ic_x,b_jc_y,`$ and $`b_jc_y`$) but they also provide implicit constraints which are implied by transitivity: if $`xy`$ then we have $`a_ic_xc_yb_j`$, and so $`a_ib_j`$. This information can be used in the alignment of the ancestral sequences by requiring potential anchors between the sequences to satisfy the constraints.

Thus when constructing the multiple alignment at node $`x`$, every triplet of sequences $`(a,b,c)`$ with $`a`$ in $`u`$, $`b`$ in $`v`$ and $`c`$ not in $`x`$ provides a potential constraint for the alignment. This can lead to a combinatorial explosion of constraints. If there are $`n`$ sequences in the alignment then there are $`O(n^3)`$ such triplets, each of which may imply many constraints. Fortunately we do not need to find the set of all possible constraints, many of which will be redundant. Instead we wish to find a set of prime constraints (i.e. a set such that no constraint is implied by the others) which is equivalent to, but potentially much smaller than, the set of all constraints implied by the gene matches. Such a set can be inferred from the homology map. If there are $`m`$ sets of orthologous exons (not all of which will be in every sequence) then at node $`x`$ there can be at most $`O(m)`$ prime constraints and a prime set which is equivalent to all possible constraints can easily be found in $`O(mk)`$ time where $`k`$ is the number of leaves below $`x`$. Thus all the set of all prime constraints can be found in $`O(mk^2)`$ time with a small constant factor. Matches between the ancestral sequences which are inconsistent with this set of constraints can then be filtered out in time $`O(Nlogm)`$ where $`N`$ is the total number of matches. For typical values of $`m`$ and $`k`$, the time taken computing and utilizing the constraints is neglible.

Figure 2 shows an example of a constraint, and how it is enforced in the AVID alignment of the ancestral sequences.

The pre-processing step of finding all exon matches is quadratic in the number of sequences, however since the protein alignments are gene based, they are typically computed on less than $`5\%`$ of the sequence. Thus, the gene matching is actually significantly faster than translated match finding which requires searching the entire sequence in all three frames and on both strands. Furthermore, by comparing only the proteins produced by a gene prediction program, the program implicitly takes into account splice sites and other gene features in building gene anchors. It is also important to note that this approach is completely ab initio even though a gene finding step is necessary, no information beyond the sequences is used. For this paper we performed alignments using this strategy in order to demonstrate the performance of an ab initio approach. However it is possible to make use of mRNA and EST data, thus incorporating known biological annotation about the sequences into the alignment.

### 2.3 Tree Building

Most multiple alignment programs require pairwise alignments of all the sequences to build an initial guide tree. This step requires a quadratic number of sequence alignments and is infeasible for large numbers of sequences. We utilize an iterative method to obtain a guide tree using only a linear number of alignments.

The initial guide tree is selected randomly from the set of complete binary trees (or almost complete binary trees in the case where the number of sequences is not a power of 2). For a given number of nodes, these are the binary trees with minimal depth and thus initial errors in pairwise alignments have less opportunity to propagate through the tree. The sequences are aligned using this random tree and then a phylogenetic tree is inferred from the resulting multiple alignment. The likelihood of the tree given the alignment can be used as a quantitative measure of the quality of the tree and the process is iterated until the alignment and tree are satisfactory.

For small numbers of sequences, the inference of the tree from the multiple alignment can be done using maximum-likelihood methods and accounts for only a small percentage of the running time. However as the number of sequences increases, we have found that ML reconstruction becomes impractical and neighbor-joining must be used. Because pairwise alignments are easy to infer from a multiple alignment, we can perform neighbor-joining reconstruction rapidly even with large numbers of sequences. We have tested MAVID with the fastDNAml \[Olsen 1994\] program for smaller datasets and the CLUSTALW implementation of neighbor joining for larger problems.

Instead of computing all pairwise alignments, only $`O(nk)`$ alignments are necessary to perform $`n`$ iterations with $`k`$ sequences. We found that for typical alignment problems only a small number of iterations were necessary (see results on the HIV sequences below). It is important to note that our iterative method (multiple alignment alternating with neighbor joining) is considerably less sophisticated than ML methods such as SEMPHY \[Friedman 2002\] or MCMC sampling methods that search through combined alignment/tree space \[Hein 2000\]. However our approach is scalable to large problems and, as we have pointed out, appears to converge quickly in practice.

## 3 Results

### 3.1 A human, mouse and rat whole genome multiple alignment

We aligned the human (April 2003), mouse (February 2003), and rat (June 2003) genomes using MAVID. A homology map for the genomes was built by Colin Dewey (manuscript in preparation), and was used to generate gene anchors and constraints. Figure 3 summarizes the exon coverage of the alignment on chromosome 20: it shows how many of the RefSeq genes were covered by anchors (and therefore automatically aligned correctly), and how many were subsequently aligned by MAVID. Chromosome 20 was chosen because it aligns almost completely with mouse chromosome 2, and therefore the quoted numbers should be useful for comparing MAVID to other alignment approaches which do not explicitly separate out orthologous from paralogous alignments.

The MAVID alignments have been used to estimate evolutionary rates for the genomes, and to identify evolutionary hotspots where one of the rodent genomes has been evolving much more slowly than the other (these results are reported in another companion paper \[Yap and Pachter\]). They are also used to support the K-BROWSER \[Chakrabarti and Pachter\], which is a new browser especially designed to view multiple genomes, their associated annotations and alignments.

### 3.2 CFTR region: 21 organisms

We aligned 1.8Mb of human sequence together with the homologous regions from 20 other organisms (baboon, cat, chicken, chimp, cow, dog, dunnart, fugu, hedgehog, horse, lemur, macaque, mouse, opposum, pig, platypus, rabbit, rat, tetraodon and zebrafish) for a total of 23Mb of sequence. This sequence has been generated by NISC as part of a comprehensive project to sequence a number of regions in the genome in multiple organisms for evolutionary and functional studies. However it is important to note that some of the sequences remain incomplete, contributing to the difficulty of the multiple alignment problem. A subset of this data-set (13 organisms) has recently been used \[Brudno 2003\] as a benchmark to compare pairwise and multiple alignment programs.

The map-building step takes approximately 15 minutes on a 2.6GhZ processor, with peak memory usage of roughly 700Mb for GENSCAN. The subsequent MAVID alignment takes another 24 minutes, for a total of about 40 minutes. The tree reconstruction step takes less than a minute using neighbor joining. Thus, an iterative approach to building the tree is feasible, and a stable tree was constructed after only two rounds of alignment.

It is difficult to assess the overall quality of the alignment, but one feature which can be verified is the alignment of exons. In order to do so, we projected the alignment onto the human sequence in order to produce pairwise alignments between human and each of the other 20 sequences. This analysis was complicated by the fact that the sequencing is not complete and so not every exon has been sequenced in every organism. In order to address this shortcoming, we calculated the fraction of human exons which were aligned with each of the sequences. An exon was considered to be aligned if at least 70% of it was covered by alignment, and at least 50% of the bases were matching.

The MAVID alignments were compared to MLAGAN, version 1.1 \[Brudno 2003\]. MLAGAN is the only other program we know of that is able to align the 21 sequences in a reasonable period of time (the running time of MLAGAN on the 21 sequences is roughly 6 hours). DIALIGN \[Morgenstern 1998\], also designed for large genomic regions, was too slow for processing the sequences: even with the new CHAOS/DIALIGN program, aligning only 4 of the sequences took 14 hours. The MGA program \[Hohl 2002\] is designed for very similar sequences only, and so was not suitable for the diverse fish, bird, marsupial mammal alignments.

Table1: Comparison of MLAGAN and MAVID multiple alignments on the 21 organism CFTR alignment. Each row in the table shows the % coverage of the alignable exons in an organism as calculated by extracting the pairwise alignment with human from the multiple alignment. The different columns in the table correspond to the number of sequences in the multiple alignnment. Thus, the first column corresponds to the pairwise alignment of human-zebrafish and the last column to the multiple alignment of all 21 sequences). Organisms were added according to the k-MST, so that the alignment problems are as difficult as possible (mutually most distant organisms).

In order to better understand how alignment accuracy varies with the number of sequences being aligned, we compared the MAVID and MLAGAN alignments on 21 different datasets, beginning with a pairwise alignment of the most distant organisms, and adding in one (mutually most distant) organism at a time all the way up to a comparison on the 21 sequences. In order to do this we computed the human clamped k-MST trees \[Boffelli 2003\], i.e. the subtrees on $`k`$ leaves of maximum weight with the human sequence as one of the leaves. Thus, alignments were computed first for human and zebrafish alone, then for human, zebrafish and dunnart, and eventually all of the 21 sequences. The order in which they were added was: human, zebrafish, dunnart, platypus, hedgehog, chicken, rat, fugu, cow, rabbit, dog, oppossum, tetraodon, lemur, horse, pig, cat, mouse, baboon, macaque, chimp. Exon coverage was calculated by first running TBLASTX to identify human exon homologs in the other species using the same criteria as in \[Brudno 2003\], and then computing coverage with respect to the identified exon sets.

The results of the alignments show that both programs correctly aligned mammalian sequences. The alignment of distant organisms shows much greater variability with respect to the sequences included in the alignment problem. For example, adding fugu to an alignment of human, zebrafish, and dunnart may improve the alignment, but as Table 1 demonstrates, adding platypus can degrade it. MAVID shows significant improvement over MLAGAN in this respect.

### 3.3 HIV1/SIV: complete genomes from 242 individuals

The HIV databases maintained at LANL contains a collection of HIV-1, HIV-2, and SIV sequences, carefully linked with individuals and their histories. We extracted the complete genomic sequences of HIV1 and SIV from this database (currently totaling 242) and aligned them with MAVID. The alignment of the sequences takes 2.5 minutes. A phylogenetic tree was constructed with neighbor joining taking an additional 30 seconds. Again, it is difficult to assess the quality of the alignment, but it is accurate enough that the different strains cluster in the inferred tree (see Figure 4). In order to understand the stability of the tree building/alignment iteration, we examined 100 alignment runs on a reduced sequence set with the recombinant strains removed (since recombination is not addressed by standard phylogenetic models). We found that the correct strains were grouped together within one round of alignment (starting with a random tree) in all the 100 runs. To validate the iterative alignment/tree building procedure, we also examined the number of sub-trees that were fixed with each successive round of alignment. The tree inferred after the second round of alignment has $`45.87`$ of its subtrees in agreement with the tree after the first run (on average), and the number fixed between the second and third rounds is $`54.08`$ on average. Starting with different random trees and aligning for three rounds, we find that $`54.16`$ of the sub-trees agree on average. Our criteria for comparing trees (exact agreement of sub-trees) is rather strict, and these numbers are very encouraging (one would expect about one non-trivial sub-tree to agree for two random trees). Although a MAVID multiple alignment combined with a neighbor joining tree may not be as accurate as hand edited alignments followed by maximum likelihood tree building, it can serve to provide very fast results that can then form the basis for further refinement. An alignment problem of this size is not practical on a standard desktop computer if all of the pairwise alignments are computed in order to build an initial guide tree (as is done in many programs, e.g. CLUSTALW).

The alignments and phylogenetic trees for all the above sequences are downloadable at
http://baboon.math.berkeley.edu/mavid/data

## 4 Conclusion

As we have outlined in the introduction, we view the genomic multiple sequence alignment problem as a biological alignment problem, rather than a purely mathematical one. That is, the incorporation of biologically relevant information (in our case ab-initio gene predictions) is critical to building accurate alignments and correctly identifying homologous relationships. Our method of incorporating constraint information into the alignments helps address one of the primary objections to progressive alignment strategies, namely that progressive alignment is “local” with the alignment at each node containing only information about the sequences below it. The application of constraint information can be thought of as a ”look ahead” step which helps to fix potential problems.

Our approach is also consistent with a number of other ideas which we have not yet implemented, but which could be easily integrated into MAVID and will improve results. Iterative refinement, the process of re-aligning across an edge in the tree, fits in naturally with our framework \[Gotoh 1993, Gotoh 1996\]. Similarly, the homology map which MAVID uses can indicate information about inversions and duplications, and this can be used to correctly align regions containing rearrangements.

MAVID compares favorably with existing programs. As we have pointed out, it is significantly more accurate than MLAGAN on the alignment of the CFTR benchmark region. MLAGAN is the only other program we know of that can even align such a large dataset. We also know of no other programs that can quickly align hundreds of viral or mitochondrial genomes.

A MAVID web server has been operational for over six months and processes over 700 requests a month \[Bray 2003A\]. Alignment requests have ranged from large genomic regions in mammals, fish, flies and plants, to alignments of viruses, mitochondria, and other bacterial genomes. MAVID can be accessed at http://baboon.math.berkeley.edu/mavid/. The program is freely available for academic and non-profit use.

## 5 Acknowledgments

We thank Colin Dewey for providing access to his homology map making programs. We thank Von Bing Yap for helping with the evolutionary models used in MAVID. Thanks to Ingileif Bryndís Hallgrímsdóttir for her help throughout the project and for her comments on the final manuscript. The data used in the multiple alignment of the CFTR region was generated by the NIH Intramural Sequencing Center (www.nisc.nih.gov), and was used subject to their six month hold policy. The HIV sequences were downloaded from the HIV database (hiv-web.lanl.gov). Thanks also to the rat sequencing consortium, both for providing the rat sequence to align, and for facilitating helpful collaborations and discussions. Finally, we thank the anonymous reviewers for their insightful comments and suggestions. This work was partially supported by funding from the NIH (grant R01-HG02362-01) and the Berkeley PGA grant from the NHLBI.

## Website References

http://www.nisc.nih.gov/ – NIH Intramural Sequencing Center
http://hiv-web.lanl.gov/ – LANL HIV Databases
http://baboon.math.berkeley.edu/mavid/ – The MAVID Webserver
http://baboon.math.berkeley.edu/mavid/data/ – Supplemental Data
http://hanuman.math.berkeley.edu/kbrowser/ – KBROWSER


format: text

text:
# Threshold Perspectives on Meson Production

## 1 Introduction

High precision data from the present generation of cooler rings, IUCF, CELSIUS, and COSY, have contributed significantly over the last decade to our present knowledge and understanding of threshold meson production (for a recent review see Moskal et al. (2002)).

Due to the high momentum transfers required to create a meson or mesonic system in production experiments close to threshold the short range part of the interaction is probed. In nucleon–nucleon scattering, for mesons in the mass range up to $`1\text{GeV}/\text{c}^2`$ distances from $`0.53\text{fm}`$ ($`\pi ^0`$) down to less than $`0.2\text{fm}`$ ($`\varphi `$) are involved. At such short distances it is a priori not clear, whether the relevant degrees of freedom are still baryons and mesons, or rather quarks and gluons. As there is no well defined boundary, one goal of the threshold production approach is to explore the limits in momentum transfer for a consistent description using hadronic meson exchange models. Within this framework, questions concerning both the underlying meson exchange contributions and especially the role of intermediate baryon resonances have to be answered.

Another aspect which enriches the field of study arises from the low relative centre–of–mass velocities of the ejectiles: Effects of low energy scattering are inherent to the observables due to strong final state interactions (FSI) within the baryon–baryon, baryon–meson, and meson–meson subsystems. In case of short–lived particles, low energy scattering potentials are otherwise difficult or impossible to study directly.

## 2 Dynamics of the Two Pion System

In $`\gamma `$ and $`\pi `$ induced double pion production on the nucleon the excitation of the $`N^{}(1440)`$ $`\text{P}_{11}`$ resonance followed by its decay to the $`N\sigma `$ channel, i.e. $`N^{}(1440)p(\pi \pi )_{I=l=\mathrm{\hspace{0.17em}0}}`$, is found to contribute non–negligibly close to threshold Oset and Vicente-Vacas (1985); Bernard et al. (1995); Gómez Tejedor and Oset (1996). Nucleon–nucleon scattering should provide complementary information, eventually on the $`\pi \pi `$ decay mode of the $`N^{}(1440)`$, which plays an important part in understanding the basic structure of the second excited state of the nucleon Morsch and Zupranski (2000); Krehl et al. (2000); Hernández et al. (2002).

Exclusive CELSIUS data from the PROMICE/WASA setup on the reactions $`pppp\pi ^+\pi ^{}`$, $`pppp\pi ^0\pi ^0`$ and $`pppn\pi ^+\pi ^0`$ Brodowski et al. (2002); Johanson et al. (2002); Pätzold et al. (2003) are well described by model calculations Alvarez-Ruso et al. (1998): For the $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$ channels, the reaction preferentially proceeds close to threshold via heavy meson exchange and excitation of the $`N^{}(1440)`$ Roper resonance, with a subsequent pure s–wave decay to the $`N\sigma `$ channel111For the $`pn\pi ^+\pi ^0`$ final state, this reaction mechanism is trivially forbidden by isospin conservation. An underestimation of the total cross section data Johanson et al. (2002) by the model predictions Alvarez-Ruso et al. (1998) might be explained by the neglect of effects from the $`pn`$ final state interaction in the calculation Alvarez-Ruso (2002).. While nonresonant contributions are expected to be small, resonant processes with Roper excitation and decay via an intermediate $`\mathrm{\Delta }`$ ($`pppN^{}p\mathrm{\Delta }\pi pp\pi \pi `$) and $`\mathrm{\Delta }\mathrm{\Delta }`$ excitation ($`pp\mathrm{\Delta }\mathrm{\Delta }p\pi p\pi `$) are strongly momentum dependent and vanish directly at threshold. Double $`\mathrm{\Delta }`$ excitation, which is expected to dominate at higher excess energies beyond $`\text{Q}=250\text{MeV}`$ Alvarez-Ruso et al. (1998) involves higher angular momenta and consequently strongly anisotropic proton and pion angular distributions. On the other hand, the Roper decay amplitude via an intermediate $`\mathrm{\Delta }`$ depends predominantly on a term symmetric in the pion momenta (eq.(1)), leading to the $`p(\pi ^+\pi ^{})_{I=l=\mathrm{\hspace{0.17em}0}}`$ channel and an interference with the direct $`N\sigma `$ decay.

Experimentally, for the reaction $`pppp\pi ^+\pi ^{}`$ at excess energies of $`\text{Q}=64.4\text{MeV}`$ and $`\text{Q}=75\text{MeV}`$ angular distributions give evidence for only s–waves in the final state, in line with a dominating $`pppN^{}pp(\pi ^+\pi ^{})_{I=l=\mathrm{\hspace{0.17em}0}}`$ process, with the initial inelastic $`pp`$ collision governed by heavy meson ($`\sigma ,\rho `$) exchange.

Roper excitation is disclosed in the $`p\pi ^+\pi ^{}`$ invariant mass distribution (Fig.1a), where the data are shifted towards higher invariant masses compared to phase space in agreement with resonance excitation in the low energy tail of the $`N^{}(1440)`$. Compared with Monte Carlo simulations including both heavy meson exchange for $`N^{}`$ excitation, and $`pp`$ S–wave final state interaction, but only the direct decay $`N^{}p(\pi ^+\pi ^{})_{I=l=\mathrm{\hspace{0.17em}0}}`$ (dotted lines), the production process involves additional dynamics, which is apparent from discrepancies especially in observables depending on the $`\pi `$ momentum correlation $`\stackrel{}{\text{k}}_{\mathrm{\hspace{0.17em}1}}\stackrel{}{\text{k}}_{\mathrm{\hspace{0.17em}2}}`$, i.e. $`\pi ^+\pi ^{}`$ invariant mass $`\text{M}_{\pi \pi }`$ (Fig.1b) and opening angle $`\delta _{\pi \pi }=\mathrm{}(\stackrel{}{\text{k}}_{\mathrm{\hspace{0.17em}1}}\stackrel{}{\text{k}}_{\mathrm{\hspace{0.17em}2}})`$ (Fig.1c). A good description of the experimental data is achieved including the $`N^{}(1440)`$ decay via an intermediate $`\mathrm{\Delta }`$ and its destructive interference with the direct decay branch to the $`N\sigma `$ channel (solid lines) in the ansatz for the Roper decay amplitude Alvarez-Ruso et al. (1998):

$$𝒜1+c\stackrel{}{\text{k}}_{\mathrm{\hspace{0.17em}1}}\stackrel{}{\text{k}}_{\mathrm{\hspace{0.17em}2}}\left(3\text{D}_{\mathrm{\Delta }^{++}}+\text{D}_{\mathrm{\Delta }^0}\right),$$
(1)
where the first term describes the direct decay, the parameter $`c`$ adjusts the relative strengths of the two decay routes, and $`\text{D}_\mathrm{\Delta }`$ denote the $`\mathrm{\Delta }`$ propagators. A fit to the data allows to determine the ratio of partial decay widths $`(\text{M}_N^{})=\mathrm{\Gamma }_{N^{}\mathrm{\Delta }\pi N\pi \pi }/\mathrm{\Gamma }_{N^{}N\sigma }`$ at average masses $`<\text{M}_N^{}>`$ corresponding to excess energies $`\text{Q}=64.4\text{MeV}`$ and $`\text{Q}=75\text{MeV}`$ relative to the $`\pi ^+\pi ^{}`$ threshold. The numerical results, $`(1264)=0.034\pm 0.004`$ and $`(1272)=0.054\pm 0.006`$, exhibit the clear dominance of the direct decay to the $`N\sigma `$ channel in the low energy region of the Roper resonance. On the other hand they indicate the strong energy dependence of the ratio from the momentum dependence in the decay branch via an intermediate $`\mathrm{\Delta }`$, which will surpass the direct decay at higher energies Pätzold et al. (2003). A model dependent extrapolation based on the validity of ansatz (1) leads to $`(1440)=3.9\pm 0.3`$ at the nominal resonance pole in good agreement with the PDG value of $`4\pm 2`$ Hagiwara et al. (2002).

Within the experimental programme to determine the energy dependence of the $`N^{}N\pi \pi `$ decay exclusive data (for details see Bashkanov et al. (2003)) have been taken simultaneously at the CELSIUS/WASA facility on both the $`pp\pi ^+\pi ^{}`$ and $`pp\pi ^0\pi ^0`$ final states. In case of the $`\pi ^+\pi ^{}`$ system the preliminary results at an excess energy of $`\text{Q}=75\text{MeV}`$ are in good agreement with the relative strength of the decay routes adjusted to an extrapolated ratio $`(1440)=3`$. However, at slightly higher excess energy ($`\text{Q}=127\text{MeV}`$) the data might be equally well described by a value $`(1440)=1`$, which is noticeably favoured at both excess energies by the data on $`\pi ^0\pi ^0`$ production, indicating distinct underlying dynamics in $`\pi ^0\pi ^0`$ and $`\pi ^+\pi ^{}`$ production. One difference becomes obvious from the isospin decomposition of the total cross section Johanson et al. (2002): An isospin $`\text{I}=1`$ amplitude in the $`\pi \pi `$ system, and accordingly a p–wave admixture, is forbidden by symmetry to contribute to the neutral pion system in contrast to the charged complement. A p–wave component was neglected so far in the analysis, since the unpolarized angular distributions show no deviation from isotropy. However, there is evidence for small, but non–negligible analysing powers from a first exclusive measurement of $`\pi ^+\pi ^{}`$ production with a polarized beam at the COSY–TOF facility Bashkanov et al. (2003); Clement et al. (2003), suggesting higher partial waves especially in the $`\pi \pi `$ system.

At higher energies, i.e. $`\text{Q}=208\text{MeV}`$ and $`\text{Q}=286\text{MeV}`$ with respect to the $`\pi ^+\pi ^{}`$ threshold, preliminary data for both $`\pi ^+\pi ^{}`$ and $`\pi ^0\pi ^0`$ from CELSIUS/WASA rather follow phase space than expectations based on a dominating $`pppN^{}pp\sigma `$ reaction mechanismBashkanov et al. (2003). At these energies, the $`\mathrm{\Delta }\mathrm{\Delta }`$ excitation process should influence observables significantly, and, thus, a phase space behaviour becomes even more surprising, unless the $`\mathrm{\Delta }\mathrm{\Delta }`$ system is excited in a correlated way.

## 3 The Proton–Proton–Eta Final State

As a general trait in meson production in nucleon–nucleon scattering, the primary production amplitude, i. e. the underlying dynamics can be regarded as energy independent in the vicinity of threshold Moalem et al. (1996); Bernard et al. (1999); Gedalin et al. (1999). Consequently, for s–wave production processes, the energy dependence of the total cross section is essentially given by a phase space behaviour modified by the influence of final state interactions. In Fig.2 total cross section data obtained in proton–proton scattering are shown for the pseudoscalar isosinglet mesons $`\eta `$ and $`\eta ^{}`$ Moskal et al. (2003a).

In both cases, the energy dependence of the total cross section deviates significantly from phase space expectations. Including the on–shell $`^1\text{S}_0`$ proton–proton FSI enhances the cross section close to threshold by more than an order of magnitude, in good agreement with data in case of $`\eta ^{}`$. As expected from kinematical considerations Moskal et al. (2002) the cross section for $`\eta `$ production deviates from phase space including the $`pp`$ FSI at excess energies $`\text{Q}40\text{MeV}`$, where the $`^1\text{S}_0`$ final state is no longer dominant compared to higher partial waves. Deviations at low excess energies seem to be well accounted for by an attractive proton–$`\eta `$ FSI (dashed line), treated phenomenologically as an incoherent pairwise interaction Bernard et al. (1999); Moskal et al. (2002); Schuberth (1995). In comparison to the proton–$`\eta ^{}`$ (Fig.2) and proton–$`\pi ^0`$ systems only the $`p\eta `$ interaction is strong enough to become apparent in the energy dependence of the total cross section Moskal et al. (2000b). In differential observables, effects should be more pronounced in the phase space region of low proton–$`\eta `$ invariant masses. However, to discern effects of proton–$`\eta `$ scattering from the influence of proton–proton FSI, which is stronger by two orders of magnitude, requires high statistics measurements, which have only become available recently Calén et al. (1999); Abdel-Bary et al. (2003); Moskal et al. (2003a):

Close to threshold, the distribution of the invariant mass of the proton–proton subsystem is characteristically shifted towards low invariant masses compared to phase space (dotted line in Fig.3). This low–energy enhancement is well reproduced by modifying phase space with the $`^1\text{S}_0`$ $`pp`$ on–shell interaction. A second enhancement at higher $`pp`$ invariant masses, i.e. low energy in the $`p\eta `$ system, is not accounted for even when including additionally the proton–$`\eta `$ interaction incoherently (dashed line). However, including a P–wave admixture in the $`pp`$ system by considering a $`^1\text{S}_0^3\text{P}_0\text{s}`$ transition in addition to the $`^3\text{P}_0^1\text{S}_0\text{s}`$ threshold amplitude, excellent agreement with the experimental invariant mass distribution is obtained (dashed–dotted line Nakayama et al. (2003)). In return, with the P–wave strength adjusted to fit the invariant mass data, the approach fails to reproduce the energy dependence of the total cross section (Fig.2) below excess energies of $`\text{Q}=40\text{MeV}`$. Preliminary calculations considering only s–waves in the final state but using a rigorous three–body treatment of the $`pp\eta `$ final state actually decrease the cross section at large values of the $`pp`$ invariant mass (dashed–double–dotted lines Fix and Arenhövel (2002)) compared to an incoherent two–body calculation within the same framework. However, close to threshold the energy dependence of the total cross section is enhanced compared to the phenomenological incoherent treatment and the data (Fig.2). Although part of this enhancement has to be attributed to the neglect of Coulomb repulsion in the $`pp`$ system, consequently overestimating the $`pp`$ invariant mass at low values, qualitatively the full three–body treatment has opposite effects compared to a P–wave admixture in the proton–proton system in view of both the total cross section as well as the $`pp`$ invariant mass distribution. In the approximate description of the total cross section by the phenomenological s–wave approach with an incoherent FSI treatment these two effects seem to cancel casually.

Close to threshold, resonance excitation of the $`\text{S}_{11}(1535)`$ and subsequent decay to the $`p\eta `$ final state is generally222The $`pnd\eta `$ excitation function has been interpreted to provide direct experimental evidence for $`\text{S}_{11}(1535)`$ excitation Calén et al. (1997). It should be noted, however, that in Peña et al. (2001) for $`\eta `$ production in proton–proton scattering short range nucleonic currents, i.e. $`\sigma `$ and $`\omega `$ exchange are found to be much stronger compared to the contribution from resonance currents. believed to be the dominant $`\eta `$ production mechanism Batinić et al. (1997); Santra and Jain (1998); Gedalin et al. (1998); Bernard et al. (1999); Fäldt and Wilkin (2001); Nakayama et al. (2002); Baru et al. (2003); Nakayama et al. (2003). In this context, the issue of the actual excitation mechanism of the $`\text{S}_{11}(1535)`$ remains to be addressed. The $`\eta `$ angular distribution is sensitive to the underlying dynamics: A dominant $`\rho `$ exchange favoured in Fäldt and Wilkin (2001) results in an inverted curvature of the $`\eta `$ angular distribution compared to $`\pi `$ and $`\eta `$ exchanges which are inferred to give the largest contribution to resonance excitation in Nakayama et al. (2002). In the latter approach the interference of the pseudoscalar exchanges in the resonance current with non–resonant nucleonic and mesonic exchange currents turns the curvature to the same angular dependence as expected for $`\rho `$ exchange. Presently, due to the statistical errors of the available unpolarised data at an excess energy of $`\text{Q}40\text{MeV}`$ Calén et al. (1999); Abdel-Bary et al. (2003) it is not possible to differentiate between a dominant $`\rho `$ or $`\pi `$, $`\eta `$ exchange, as discussed in Abdel-Bary et al. (2003). Data recently taken at the CELSIUS/WASA facility with statistics increased by an order of magnitude compared to the available data might provide an answer in the near future Zlomanczuk (2003).

Spin observables, like the $`\eta `$ analyzing power, should even disentangle a dominant $`\rho `$ meson exchange and the interference of $`\pi `$ and $`\eta `$ exchanges in resonance excitation with small nucleonic and mesonic currents Nakayama et al. (2002), which result in identical predictions for the unpolarised $`\eta `$ angular distribution. First data Winter et al. (2002) seem to favour the vector dominance model, but final conclusions both on the underlying reaction dynamics and the admixture of higher partial waves Nakayama et al. (2003) have to await the analysis of data taken with higher statistics for the energy dependence of the $`\eta `$ analysing power Czyżykiewicz et al. (2003).

## 4 Associated Strangeness Production

In elementary hadronic interactions with no strange valence quark in the initial state the associated strangeness production provides a powerful tool to study reaction dynamics by introducing a “tracer” to hadronic matter. Thus, quark model concepts might eventually be related to mesonic or baryonic degrees of freedom, with the onset of quark degrees of freedom expected for kinematical situations with large enough transverse momentum transfer.

First exclusive close–to–threshold data on $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^0`$ production Balewski et al. (1998); Sewerin et al. (1999) obtained at the COSY–11 facility showed at equal excess energies below $`\text{Q}=13\text{MeV}`$ a cross section ratio of

$$_{\mathrm{\Lambda }/\mathrm{\Sigma }^0}\left(\text{Q}13\text{MeV}\right)=\frac{\sigma \left(pppK^+\mathrm{\Lambda }\right)}{\sigma \left(pppK^+\mathrm{\Sigma }^0\right)}=28_9^{+6}$$
(2)
exceeding the high energy value ($`\text{Q}300\text{MeV}`$) of 2.5 Baldini et al. (1988) by an order of magnitude.

In the meson exchange framework, estimates for $`\pi `$ and $`K`$ exchange contributions based on the elementary scattering processes do not reproduce the experimental value (2) Sewerin et al. (1999); Gasparian et al. (2000). However, inclusive $`K^+`$ production data in $`pp`$ scattering at an excess energy of $`\text{Q}=252\text{MeV}`$ with respect to the $`pK^+\mathrm{\Lambda }`$ threshold show enhancements at the $`\mathrm{\Lambda }p`$ and $`\mathrm{\Sigma }N`$ thresholds of similar magnitude Siebert et al. (1994). Qualitatively, a strong $`\mathrm{\Sigma }^0N\mathrm{\Lambda }p`$ final state conversion might account for both the inclusive SATURNE results as well as the $`\mathrm{\Sigma }^0`$ depletion in the COSY–11 data. Evidence for such conversion effects is known e. g. from fully constrained kaon absorption on deuterium via $`K^{}d\pi ^{}\mathrm{\Lambda }p`$ Tan (1969).

In exploratory calculations performed within the framework of the Jülich meson exchange model Gasparian et al. (2000), taking into account both $`\pi `$ and $`K`$ exchange diagrams in a coupled channel approach, a final state conversion is rather excluded as origin of the experimentally observed ratio: While $`\mathrm{\Lambda }`$ production is found to be dominated by kaon exchange, both $`\pi `$ and $`K`$ exchange turn out to contribute to the $`\mathrm{\Sigma }^0`$ channel with similar strength. Qualitatively, this result is experimentally confirmed at higher excess energies between $`\text{Q}=200\text{MeV}`$ and $`\text{Q}=430\text{MeV}`$ from polarization transfer measurements from the DISTO experiment Balestra et al. (1999); Maggiora (2002); Laget (1991). It is concluded in Gasparian et al. (2000), that only a destructive interference of $`\pi `$ and $`K`$ exchange might explain the experimental value (2). $`\mathrm{\Sigma }`$ production in different isospin configurations should provide a crucial test for this interpretation, since for the reaction $`ppnK^+\mathrm{\Sigma }^+`$ the interference pattern is found to be opposite compared to the $`pK^+\mathrm{\Sigma }^0`$ channel. Data close to threshold have recently been taken at the COSY–11 facility Rożek and Grzonka (2002).

However, within an effective Lagrangian approach Shyam et al. (2001) both $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^0`$ production channels are concluded to be dominated by $`\pi `$ exchange and excitation of the $`\text{S}_{11}(1650)`$ close to threshold, while at excess energies above $`\text{Q}=300\text{MeV}`$ the $`N^{}(1710)`$ governs strangeness production333For further complementary theoretical approaches see references in Moskal et al. (2002); Kowina et al. (2003a); Wagner et al. (2003).. In this energy range the influence of resonances becomes evident from recent data on invariant mass distribution determined at COSY–TOF Wagner et al. (2003).

To study the transition region between the low–energy enhancement (2) and the high energy value measurements have been extended up to excess energies of $`\text{Q}=60\text{MeV}`$ Kowina et al. (2003a, b): In order to describe the energy dependence of the total cross section for $`\mathrm{\Lambda }`$ production, in addition to phase space the $`p\mathrm{\Lambda }`$ final state interaction has to be taken into account. In contrast, $`\mathrm{\Sigma }^0`$ production is satisfactorily well described by phase space behaviour only Kowina et al. (2003a). This qualitatively different behaviour might be explained by the $`\mathrm{\Sigma }^0p`$ FSI being much weaker compared to the $`\mathrm{\Lambda }p`$ system. However, the interpretation implies dominant S–wave production and reaction dynamics that can be regarded as energy independent. Within the present level of statistics, contributions from higher partial waves can be neither ruled out nor confirmed at higher excess energies for $`\mathrm{\Sigma }^0`$ production.

The energy dependence of the production ratio $`_{\mathrm{\Lambda }/\mathrm{\Sigma }^0}`$ is shown in Fig.4 in comparison with theoretical calculations obtained within the approach of Gasparian et al. (2000) assuming a destructive interference of $`\pi `$ and $`K`$ exchange and employing different choices of the microscopic hyperon nucleon model to describe the interaction in the final state Gasparyan (2002).

The result crucially depends on the details — especially the off–shell properties — of the hyperon–nucleon interaction employed. At the present stage both the good agreement found in Gasparian et al. (2000) with the threshold enhancement (2) and for the Nijmegen model (dashed line in Fig.4) with the energy dependence of the cross section ratio should rather be regarded as accidental444In the latter case an SU(2) breaking in the $`^3\text{S}_1`$ $`\mathrm{\Sigma }N`$ channel had to be introduced Maessen et al. (1989) resulting in an ambiguity for the $`\mathrm{\Sigma }^0p`$ amplitude.. Calculations using the new Jülich model (solid line in Fig.4) do not reproduce the tendency of the experimental data. It is suggested in Gasparyan (2002) that neglecting the energy dependence of the elementary amplitudes and higher partial waves might no longer be justified beyond excess energies of $`\text{Q}=20\text{MeV}`$. However, once the reaction mechanism for close–to–threshold hyperon production is understood, exclusive data should provide a strong constraint on the details of hyperon–nucleon interaction models.

## 5 Present and Future

Intermediate baryon resonances emerge as a common feature in the dynamics of the exemplary cases for threshold meson production in nucleon–nucleon scattering discussed in this article. However, this does not hold in general for meson production in the $`1\text{GeV}/\text{c}^2`$ mass range (for a discussion on $`\eta ^{}`$ production see Moskal et al. (2003b)). Moreover, the extent to which resonances are evident in the observables or actually govern the reaction mechanism depends on the specific channels, which differ in view of the level of present experimental and theoretical understanding.

The $`N^{}(1440)`$ resonance dominates $`\pi ^+\pi ^{}`$ production at threshold, and exclusive data allow to extract resonance decay properties in the low–energy tail of the Roper. Dynamical differences between the different isospin configurations of the $`\pi \pi `$ system and the behaviour at higher energies remains to be understood with first experimental clues appearing.

With three strongly interacting particles in the final state, a consistent description of $`\eta `$ production close to threshold requires an accurate three–body approach taking into account the possible influence of higher partial waves. High statistics differential cross sections and polarization observables coming up should straighten out both the excitation mechanism of the $`N^{}(1535)`$ and the admixture of higher partial waves.

At present, the available experimental data on the elementary strangeness production channels give evidence for both an important role of resonances coupling to the hyperon–kaon channels and on a dominant non–resonant kaon exchange mechanism. Experiments on different isospin configurations, high statistics and spin transfer measurements close to threshold should disentangle the situation in future.

From the cornerstone of total cross section measurements, it is apparent from the above examples to what extent our knowledge is presently enlarged by differential observables and what will be the impact of polarization experiments in future to get new perspectives in threshold meson production.

The author gratefully acknowledges the pleasure to work with the CELSIUS/WASA and COSY–11 collaborations, and, in particular, thanks M. Bashkanov, H. Clement, R. Meier, P. Moskal and W. Oelert for helpful discussions . This work has been supported by The Swedish Foundation for International Cooperation in Research and Higher Education (STINT Kontrakt Dnr 02/192).