text:
# Untitled Document

Results and limits of the soliton theory of DNA transcription

Paper presented at the conference Nonlinear Phenomena in Biology, Pushchino, Russia, 22-27 June 1998
Giuseppe Gaeta

Dipartimento di Fisica, Universitá di Roma

00185 Roma (Italy)

gaeta@roma1.infn.it

Abstract.  It has been suggested by several authors that nonlinear excitations, in particular solitary waves, could play a fundamental functional role in the process of DNA transcription, providing the opening of the double chain needed for RNA Polymerase to be able to copy the genetic code. Some models have been proposed to model the relevant DNA dynamics in terms of a reduced number of effective degrees of freedom. Here I discuss advantages and disadvantages of such an approach, and discuss in more detail one of the models, i.e. the one proposed by Yakushevich.

1. Introduction

The first step in the replication of DNA  is its transcription, from the original contained in the cell to a copy – the RNA messenger – which will then be used as a “master copy” for producing actual copies of the genetic information. The evolutionary advantage of such a messenger is obvious: in this way, the original DNA is opened – and thus less protected – for as small a time as possible.

The transcription process is carried out by a specialized enzyme, the RNA Polymerase (RNAP); here we are mostly interested in the dynamical aspect of the transcription, which is roughly speaking as follows. The RNAP opens a ”transcription bubble” of a size of about 15-20 base pairs, and then travels along the DNA chain keeping the size of the open region more and less constant, i.e. providing at the same time to open the chain in front of it and closing back the one behind.

The process undergoes several (even quite long in time) stops, but in the active phase the RNAP proceeds along the DNA chain at a speed of 500 – 1000 base pairs per second. Since each base pair is linked by two or three hydrogen bonds, the energy involved in such a process, even considering only the one to open (and close) the DNA chain, is of the order of thousands H bonds per second.

Obviously, as biological systems live at a temperature of the order of 300 Kelvin, thermal energy is widely available, but the problem is how is this energy focused at the right place and with the right timing to operate the process.

Needless to say, the problem of energy transport and focusing is of more general relevance in Biology and biological systems – and one indeed which has to be studied by Biological Physics. Thus, the study of this problem in the dynamics of DNA transcription is not only relevant per se and for the fundamental role of DNA transcription and then replication in life, but can also have a more general relevance for the study of energy transport in living systems.

It has been suggested by several authors that Nonlinear Dynamics can have a word to say concerning this, and that energy transport could happen by means of nonlinear excitations travelling along biological chains, such as long proteins or indeed DNA.

The most widely known theory in this context is probably that of ”Davydov’s soliton” , which is of essentially quantum nature. Theories concerning nonlinear excitations in DNA chains and their functional role have been proposed by several authors, following the seminal work of Englander, Kallenbach, Heeger, Krumhansl and Litwin  (and sometimes it is difficult to understand priorities, so I will not assign any); these present differences, but many more common aspects.

Excellent reviews exist on Nonlinear Waves in biological systems ; I would also like to mention, as a source of a number of inspiring material on the relations between nonlinear structures and biology (or more specifically biological molecules) the volume collecting the proceedings of the Les Houches 1994 School .

My goal here is not to review and/or compare the different theories proposed for these would-be excitations of DNA connected with transcription, but rather discuss the merit and limitations of the general idea of what I will call, with a deliberately ambiguous expression, ”soliton-helped transcription”.

Roughly speaking, the idea is that there are nonlinear excitations travelling along the DNA chains, causing a local opening of the double chain. The RNAP could then travel along with these, and use the opening of the chain to read the DNA sequence without having to focus the energy needed to open the double chain.

Another advantage of this theory is connected with the closing of the double chain after the RNAP has passed: if this is done in a non-coordinated way, it will generate a substantial quantity of random motion and thus of thermal energy; if instead this correspond to the passing of a nonlinear excitation (maybe in this context it is more appropriate to call our soliton a “nonlinear coherent structure”), no thermal energy is generated when the system recovers its local fundamental state.

I will of course also discuss quantitative results and predictions, and not just qualitative ideas; in doing this I will refer to a specific model, the one suggested by L.V. Yakushevich  in the late eighties, on which I have done some work and which is well suited to analytical investigation. Needless to say, the very fact that it is amenable to analytic results means that it is an oversimplified model, and it is not difficult to find in the literature – or to build by oneself – more detailed models. However, as I will advocate in more detail in the sequel, I believe the study of such simple models can teach us a lot.

Before proceeding, there is a point that should be made clear for the sake of the more mathematical reader: the nonlinear excitations we wish to consider are not necessarily solitons in the proper mathematical sense , but rather solitary waves; however, it is by now common – and more attractive – to call them solitons, and I will conform to this usage.

Also, the word “soliton” is used in two senses, a topological one  together with the dynamical one . The solitons to be considered here will be both solitary waves in dynamical sense and $`S^1`$-topological solitons.

2. What should a model contain ?

The DNA molecule is an awfully (or wonderfully) complicated one; each segment corresponding to a base pair has about 100 degrees of freedom, and a DNA molecule can contain $`10^{10}`$ (or even $`10^{13}`$ as for salamanders) such segments, thus having $`10^{12}`$ degrees of freedom.

As if this was not enough, DNA is very sensitive on the ”boundary conditions”, i.e. to the characteristics (physical and chemical) of its environment, and no chemist or biochemist would think of DNA “in abstract”, i.e. without specifying in which environment it operates.

With such numbers (even apart from the problem of modelling the interactions of the DNA molecule with its environment) it is clear that there is no hope of detailed mathematical modelling. Actually, these are so large that we are quite justified in applying the thermodynamical limit (i.e. consider an infinite chain), and resorting to statistical mechanics ; indeed, several approaches have been attempted in this direction, and with some success.

However, the dynamical models which have been proposed to study the dynamics of DNA transcription, only consider an extremely limited number of degrees of freedom (one per base in the Yaushevich model); indeed, it seems hopeless to study a model with more than a very few degrees of freedom per segment and investigate if it has solitons, not to say determine the detailed dynamics of these.

Thus the first question we should ask is: does it make any sense to study such simplified models of such an enormously complicated system as DNA ?

Substituting simple mathematical models to complicate real system is, of course, what has always been done in theoretical – and not only theoretical – Physics, and it has been successful many times; but this is of course not a good reason to say that it can be done here as well, not to say that in a biological system the complexity is somehow much more fundamental, inherent to the system, than in Physics, where we can e.g. study gravitation without the disturbances of air by dropping masses in a vacuum tube: in Biology, “simplifying the system” most probably means killing it, and biologists and biochemists are well justified in their diffidence towards “Physics-style” simplified models.

I would like at this point to open a short parenthesis to mention the problem of dissociation (or melting) of DNA: experimental data shows there is a very sharp transition as temperature is raised, even in relatively short specimen of DNA, leading to dissociation. For some time this has been a serious obstacle in the attempt to model DNA as a one-dimensional system, since one would expect that there is no phase transition for one-dimensional chains .

However, the theorem on the absence of phase transitions in one dimension only apply to Ising-like (spin) systems, while the situation is far more delicate for systems the state of whose components is described by a continuous local order parameter. Basically, in the former case one describes the system by means of a transfer matrix  (which is a Markov one, so relevant restriction on its eigenvalues exist – in particular, the fundamental state cannot be degenerate, and so phase transitions are not possible), while in the latter one one has to use a transfer operator ; it is by making use of this formalism, following work done by J.A. Krumhansl and J.R. Schrieffer in 1975, that A. Bishop, Th. Dauxois, and M. Peyrard proved the existence of a ”dissociation” phase transition in DNA considered as a one dimensional system ; see also .

I do not want to discuss their work here – among other reasons, because I could only provide a bad version of their excellent papers – but I want to stress that their success was not only of mathematical or qualitative nature: their theory compares successfully to experimental data on the detailed (spatiotemporal) dynamics of DNA melting; I want to stress again that they are able to predict not only average quantities, as it should anyway be the case with a Statistical Mechanics approach, but a spatiotemporal pattern: this is a much more severe test of a model, and the model of Peyrard and Bishop  for DNA melting which was used (with some minor modifications) in this study did brilliantly pass the test.

Now, what are the characteristics of this model ? Well, it is pretty much similar to the models we consider for transcription, and in particular to the Yakushevich one. Indeed it models DNA as a one-dimensional chain, and by singling out one degree of freedom per base – corresponding to ”radial” displacements along the axis joining the two bases of a pair – that is, the degree of freedom thought to be the most relevant for the process under study.

Thus the success of Bishop, Peyrard and Dauxois in the study of DNA melting is a witness to the fact that these – apparently, hopelessly oversimplified – models can indeed provide a good understanding, qualitative but also quantitative, of the behaviour of DNA in biologically relevant processes. Simple models can be relevant even in very complicated systems !

I would take this to the extreme, and say that the simplest the model, the greater its value: not only because a simple model can be studied more throughly, but first and foremost because if a simple model works, it means it has focused on the right point, i.e. that we have reached an understanding of what are the most relevant aspects (degrees of freedom) of our system for the process under study, and thus it adds more to our knowledge.

Now, conforted by the success of Bishop, Peyrard and Dauxois, we can look back at the simple counting of degrees of freedom in DNA with a more optimistic attitude: it is true that DNA is terribly complicated, but it is also true that it performs a huge variety of essential tasks. While the amazement for how seemingly inconciliable characteristics needed to perform these different tasks can convive in a single molecule should remain, it is possible to think of a process-oriented modelling of DNA, i.e. to obtain a model of DNA which does not have any ambition – or hope – to describe it in general, but only limitedly to a specific (although relevant) process we want to study.

This of course is dependent upon our ability to correctly identify the specific aspects (or degrees of freedom) of DNA which are essential for the process we study. This is the problem whenever we try to formulate a model, and the Peyrard-Bishop model shows that it is possible to do this even in DNA.

A less extreme reduction – keeping a few degrees of freedom per base rather than only one – has been recently considered by Zhang and Collins , again obtaining quite encouraging results. In particular, they compared results obtained by a simplified models with those obtained by an all-atoms simulations, showing that modelling DNA dynamics in terms of a few degrees of freedom is by all means satisfactory, at least for the analysis of specific processes.

Thus, let us come back to discuss what a model of DNA nonlinear excitations which could play a role in transcription should contain.

A first obvious but relevant observation is that the DNA molecule, like any other one, obeys Quantum Mechanics (it should be recalled that both the Peyrard-Bishop and the Yakushevich model are classical); moreover, we are interested in its behaviour in living condition, i.e. in particular at a quite precise temperature (around 300 Kelvin, or a bit more for humans).

Now it happens that at this temperature most of the degrees of freedom of the atoms constituting the DNA momlecule will be essentially frozen due to Quantum Mechanics (their excitation energy being much higher), and thus many bonds between atoms can be considered as rigid; on the other side, there are degrees of freedom whose excitation energy is much lower than the thermal scale, so that they can be considered to behave classically.

The consequence of this on modelling is clear: we can safely consider the degrees of freedom of the first kind as non-existent, and concentrate on the remaining – effective – degrees of freedom. These should still be treated quantistically, but for some of them – those of the second kind – we can also use a classical modellization.

It should be clear that the (quite general) reduction to “effective” degrees of freedom is by no means sufficient if we want to consider only one – or even a few – degree of freedom per base: we should have some physical understanding of the process we want to model in order to understand which are the relevant degrees of freedom, which cannot always be just a few.

In the case of transcription the relevant process undergone by the conformation of DNA is its “unwinding”  so that the bases, which usually are stored in the interior of the molecule, as protected as possible from external agents, come to be accessible for reading by the RNAP.

Among the degrees of freedom which are relevant to such a conformational change, there is one which is “softer” than any other, corresponding grosso modo to rotations of the bases in a plane orthogonal to the axis of the double helix. Actually such a rotation is, obviously, more complicate than this, but we will investigate the hypothesis that the comparatively small movements around this main rotation are not of fundamental importance (see ).

It is also reasonable, looking at the geometry of the molecule, to model the conformational change only in terms of this; obviously a lot of other degrees of freedom are involved in such a conformational change, but the idea is that they will play the role of “slave modes”, i.e. there will be small adjustements involving these but which essentially follow the relevant changes undergone by the “master modes”, those of the fundamental degree of freedom.

Such a way of thinking, in terms of relevant variables and slave ones, has been successfull in a wealth of situations, and in Statistical Mechanics  or in Hydrodynamics  one has a very precise way of defining them; however, there is no a priori argument showing that this should be the right way of looking at this specific problem of DNA transcription: we can only justify it a posteriori if the model works properly.

On the other side, if the model works, this means not only that we have the correct values for coupling constants and so on, but also something more fundamental: that we have individuated the most relevant degree of freedom for the process we are studying.

3. The model and its prediction.

I want now to briefly present the model which has been proposed by L.V. Yakushevich (with some small modification), and its predictions concerning the existence of solitons and their characteristics.

The discussion of the limitations, or of how acceptable are a number of approximations it embodies, will be postponed to the following section.

As I mentioned in the Introduction, several other models have been proposed; I will mention here  (and apologize to other authors I have inadvertentdly forgotten); in many cases these give predictions similar to those of the Yakushevich model, in particular for the successful ones, but my goal here is not to provide a review of the different models available, and I prefer to restrict my discussion to the model I am familiar with. This is also justified by the fact that this model has been quite successfull; and I do not think one can obtain much more without introducing more detailed models and designing specific experiments to test their (presumably, only numerical) predictions. If we look at the Yakushevich “classification” of DNA models , I would say that her model provides what is reasonable to ask to a model of its level, and that one should now pass to a deeper level to really test the theory and discriminate among different versions.

In the Yakushevich model, we consider rotations of the bases in a plane orthogonal to the double helix axis (which of course is not completely true), and any other movement is not considered. Notwithstanding the fact that the forces involved in hydrogen bonds are highly directional, the forces among the bases are considered to be linear in the distance $`d`$ between the extremities interacting to form the H bonds. This will produce a nonlinear dependence on the rotation angles.

An even more severe approximation is introduced concerning the characteristics of the bases and their interactions: all bases are considered as identical, and so the same also holds for base-pair interactions.

Let us consider the interaction among bases on the two chains at site $`n`$. If $`r`$ is the distance from the center of rotation to the site responsible for the intra-pair interaction, $`\theta _1(n)`$ and $`\theta _2(n)`$ the angles (both measured in, say, counterclockwise direction ) of rotation, and $`L`$ the distance at rest, we would then have for the distance $`d`$,

$$d^2=[2L+2rr(\mathrm{cos}\theta _1+\mathrm{cos}\theta _2)]^2+[r(\mathrm{sin}\theta _1+\mathrm{sin}\theta _2)]^2.$$
$`(1)`$
However, in the Yakushevich model one makes the approximation $`L=0`$ (see next section).

Thus the potential “rotation” energy will be $`V_r=(K_r/2)d^2`$, with $`K_r`$ a real constant and $`d^2`$ the simplified expression corresponding to taking $`L=0`$ in (1). Notice that the force this produces is nonlinear in $`\theta _1,\theta _2`$.

Together with this, we will have an elastic interaction with nearest neighbour bases along the two chains, with interaction “torsional” energy $`V_t=(K_t/2)(\theta _1(n)\theta _1(n1))^2`$ and the like for the other chain and for the interaction with bases at site $`(n+1)`$.

There is another kind of interaction which should be taken into account: bases which are half-pitch of the helix apart are actually near enough in ambient space, and they interact through the formation of water filament; the energy associated to this “helicoidal” interaction is $`V_h=(K_h/2)[\theta _1(n+p)\theta _2(n)]^2`$ and the like for the other chain, interchaninging indices 1 and 2, and for interaction with bases at site $`(np)`$. Here $`p`$ is the half-pitch of the helix measured in units of base sites, which for DNA means $`p=5`$. This interaction, although weaker than the other ones, produces a qualitative change in the dispersion relations and in the predictions of the model.

Thus, introducing also the kinetic energy $`K=I/2\dot{\theta }_1(n)+I/2\dot{\theta }_2(n)`$, we can finally write the total energy for the (Yakushevich model of the) DNA double chain:

$$\begin{array}{cc}\hfill H=& \underset{n}{}\frac{1}{2}I\left[\dot{\theta }_1^2(n)+\dot{\theta }_2^2(n)\right]+\hfill \\ & +\underset{n}{}\frac{1}{2}K_t\left[\left(\theta _1(n)\theta _1(n1)\right)^2+\left(\theta _2(n)\theta _2(n1)\right)^2\right]+\hfill \\ & +\underset{n}{}\frac{1}{2}K_r\left[6r^2+2r^2\mathrm{cos}(\theta _1(n)\theta _2(n))4r^2(\mathrm{cos}(\theta _1(n))+\mathrm{cos}(\theta _2(n)))\right]+\hfill \\ & +\underset{n}{}\frac{1}{2}K_t\left[\left(\theta _1(n)\theta _2(np)\right)^2+\left(\theta _2(n)\theta _1(np)\right)^2\right].\hfill \end{array}$$
$`(2)`$
¿From this we can easily derive the equations of motion for the system; these are more usefully written in terms of the variables

$$\psi _n:=\frac{\theta _1(n)+\theta _2(n)}{2},\chi :=\frac{\theta _1(n)\theta _2(n)}{2}$$
$`(3)`$
and are

$$\begin{array}{cc}\hfill \ddot{\psi }_n=& \alpha \mathrm{sin}\psi _n\mathrm{cos}\chi _n+\beta \left(\psi _{n+1}2\psi _n+\psi _{n1}\right)+\gamma \left(\psi _{n+p}2\psi _n+\psi _{np}\right)\hfill \\ \hfill \ddot{\chi }_n=& \alpha \mathrm{sin}\chi _n(\mathrm{cos}\psi _n\mathrm{cos}\chi _n)+\beta \left(\chi _{n+1}2\chi _n+\chi _{n1}\right)\gamma \left(\chi _{n+p}+2\chi _n+\chi _{np}\right)\hfill \end{array}$$
$`(4)`$
(notice the sign differences in the $`\gamma `$ terms), where $`\alpha =K_r/I`$, $`\beta =K_t/I`$ and $`\gamma =K_h/I`$.

In these equations we find the different elastic constants $`K_r,K_t,K_h`$ or equivalently $`\alpha ,\beta ,\gamma `$; these could be seen as free parameters to be fitted, but it is also possible to give an estimate on their values based on first principles, i.e. on an estimate of the strength of the interactions they should model. The values given in  are as follows:

$$K_r=0.13\mathrm{eV}/\mathrm{rad}^2,K_t=0.025\mathrm{eV}/\mathrm{rad}^2,K_h=0.009\mathrm{eV}/\mathrm{rad}^2,I=310^{37}\mathrm{cm}^2\mathrm{g}.$$
$`(5)`$
Now, in the analysis of the Yakushevich model it has been usual to pass to the continuum approximation, i.e. having fields $`\psi (x,t),\chi (x,t)`$ instead of the sequences of variables $`\{\psi _n(t)\}`$ and $`\{\chi _n(t)\}`$; in this way the chains of ODEs (4) are replaced by coupled PDEs

$$\begin{array}{cc}\hfill \psi _{tt}=& \alpha \mathrm{sin}\psi \mathrm{cos}\chi +\beta \delta ^2\psi _{xx}+\gamma 𝒲_+(\psi )\hfill \\ \hfill \chi _{tt}=& \alpha \mathrm{sin}\chi (\mathrm{cos}\psi \mathrm{cos}\chi )+\beta \delta ^2\chi _{xx}\gamma 𝒲_{}(\chi )\hfill \end{array},$$
$`(6)`$
where $`\delta `$ represents the distance between consecutive bases along the double chain axis ($`\delta 3.4`$ Angstroms), and the nonlinear terms $`𝒲_\pm `$ are defined by

$$𝒲_\pm (\varphi ):=\varphi (x+p\delta ,t)2\varphi (x,t)+\varphi (xp\delta ,t).$$
$`(7)`$
This passage to the continuum is questionable, and indeed in their study of the melting transition Peyrard and Dauxois have pointed out interesting effects due precisely to the discrete nature of the model (there is a Peierls barrer restraining energy flow); we will discuss this point in next section as well.

We can then analyze (6); first of all – considering the linear terms – we obtain the dispersion relations for the $`\psi `$ and the $`\chi `$ branch; these are respectively

$$\begin{array}{cc}\hfill \omega _\psi ^2=& \alpha +\beta \delta ^2q^2+4\gamma \mathrm{sin}^2(p\delta q/2)\hfill \\ \hfill \omega _\chi ^2=& \beta \delta ^2q^2+4\gamma \mathrm{cos}^2(p\delta q/2).\hfill \end{array}$$
$`(8)`$
Notice that now (due indeed to the introduction of the “helicoidal” interactions) the minimum of the dispersion curve, and thus the threshold for excitation of linear modes, corresponds in terms of the natural wavenumber unit $`\xi =p\delta q/2`$, to $`\xi =1.4`$; the corresponding wavelength will be the first to be excited and should act as a seed for the formation of nonlinear structure. Thus, it gives a rough estimate for the size of excitations, which is $`\lambda 11\delta `$.

Similarly, the corresponding $`\omega =\omega (\xi _0)`$ gives a rough estimate for the timescale of small oscillations, which is of 57 ps, and similarly one obtains an estimate of the amplitude of “linear” oscillations, which is $`0.77\pi `$. This seems sufficient to ignite the unwinding and the nonlinear regime.

We stress that the minima of the dispersion relations are essentially dependent on the helicoidal interaction term (without this we actually would have phonons in the model). Data for linearized dynamics are in a way even more significant when we analyze DNA melting rather than transcription, but this lies outside our scope here; see  for a discussion. Also, the presence of helicoidal interactions lead to crossings between the two branches of the dispersion relations, which could play a role in energy transfers .

It is remarkable that the estimate of the size of collective motions is of the same order of magnitude as experimental data on the size of transcription bubbles. Again, for melting models the agreement is significant also for other relevant quantities, again providing the right order of magnitude.

Let us now look more precisely to solitons, or solitary waves, in (6). These will be solution of the form $`\varphi (x,t)=\varphi (xvt)\varphi (z)`$, where $`\varphi =\psi ,\chi `$. Moreover they should satisfy

$$\underset{z\pm \mathrm{}}{lim}\varphi _z=\mathrm{\hspace{0.17em}0},\underset{z\pm \mathrm{}}{lim}\varphi =\mathrm{\hspace{0.17em}2}m_\pm \pi $$
$`(9)`$
in order to have a finite energy. They can thus be characterized in terms of their winding numbers $`\mu =(m_+m_{})`$ for $`\psi `$ and $`\chi `$ ; these corresponds to the phase shift of the solitons.

The lowest energy solitons will be those with $`\mu =(1,0)`$ and $`\mu =(0,1)`$, and we concentrate on these.

If we set the helicoidal interaction to zero, it is possible to determine exact solitary wave solutions. These are given by

$$\psi (z)=\mathrm{arccos}\left(\frac{\mathrm{sinh}^2(az)1}{\mathrm{sinh}^2(az)+1}\right)$$
$`(10)`$
and $`\chi (z)=0`$ in the case of the $`(1,0)`$ soliton, while the $`(0,1)`$ soliton is given by $`\psi (z)=0`$ and

$$\chi (z)=\mathrm{arccos}\left(\frac{a^2z^21}{a^2z^2+1}\right).$$
$`(11)`$
Here $`a=[2\alpha /(\beta \delta ^2v)]^{1/2}`$, where $`v`$ is the soliton’s speed.

A numerical integration  shows that the introduction of the (weak) helicoidal terms does not really alter the fully nonlinear structures.

The solitonic solution can exist, in principle, with any speed $`v<\beta \delta ^2`$. Moreover, the dependence of their energy on the speed  is given by a relation of the type $`E=c_1/(c_2v^2)`$, as it should for waves, so that up to high speed the energy is very little sensitive to speed; see  for details. This means that in this model there is no selection of soliton speed. A refinement of the model  seems to be able to provide such a selection based on some conditions at interface of A-T and G-C regions, but in this model the soliton can not actually exist in proper terms, i.e. as solitary waves travelling along the chain with no dispersion.

As for the size of the soliton, this will correspond to the size $`s`$ of the region in which $`\psi `$ differs from its asymptotic value by more than a given quantity $`\epsilon `$, and thus will depend on this. We call $`z_+,z_{}`$ the values of $`z`$ at which $`\psi (z_{})=\epsilon `$, $`\psi (z_+)=2\pi \epsilon `$. For $`az>>1`$, the argument of the arccos in (10) is $`f(z)12\mathrm{exp}[2az]`$, and thus if $`\mathrm{cos}(\epsilon )=1\lambda (\epsilon )`$, we have to solve $`\lambda (\epsilon )=2\mathrm{exp}[2az_\pm ]`$; obviously $`z_{}=z_+`$, and so $`s=|a^1\mathrm{ln}[\lambda (\epsilon )/2]|`$.

As for $`a`$, when $`v=0`$ we get $`a=\sqrt{2\alpha /(\beta \delta ^2)}:=a_0\delta ^1`$, and the soliton size measured in base sites units is thus

$$s=\left|\frac{1}{a_0}\mathrm{ln}[\lambda (\epsilon )/2]\right|\delta .$$
$`(12)`$
If we choose $`\epsilon =5^o0.044`$, we get $`\lambda (\epsilon )2.910^7`$ and $`\mathrm{ln}[\lambda (\epsilon )/2]15.75`$; with the values given in (5) we get $`a_03.22`$ and thus it results $`s5\delta `$.

This is smaller than the observed size of the transcription bubble, which is of order 15 - 20 bases, but again we have got the right order of magnitude.

4. Limitations of the model.

In the previous section we have seen that the Yakushevich model, with the introduction of “helicoidal” terms, is able to give a prediction of the order of magnitude of several physical quantities for coherent structures and excitations in the DNA chain. This is specially significative when we consider the crude approximations which have been introduced in the formulation of the model, so that we are encouraged to think that the model does indeed correctly identify the essential degrees of freedom for the dynamics related to the transcription process.

In this section we want instead to emphasize the limitations and the shortcomings of the model; in most cases this will actually be not specific to the Yakushevich model, but common to the whole set of simple models we have advocated above.

We will actually first of all mention a shortcoming of teh Yakushevich model itself, which was noticed by Gonzalez and Martin Landrove : the seemingly unoffensive approximation $`L=0`$, see eq.(1), is actually selecting a special case: any nonzero value of $`L`$ would cause qualitative (and not only quantitative) differences in some basic behaviour of the model; in particular the pairing interaction would be through nonlinear terms only, with no linear ones. As this situation is unphysical, we see that the apparently special choice $`L=0`$ does actually produce a “generic” behaviour; this has also suggested to consider the Yakushevich model as an effective one.

Having mentioned this fact, we can pass to aspects which are not really specific of the Yakushevich model, although we will refer to it for the sake of definiteness.

A first obvious limitation is in the fact that bases are considered as identical. Needless to say, this is not the case. If one could be tempted to think that distinguishing between Purines and Pyrimidines only could be acceptable, as a first approximation, in view of the modellization adopted (the physical quantities entering in the model are homogeneous enough for the two purines and for the two pyrimidines, see ), considering all bases as equivalent is highly questionable: e.g., masses range from 110 to 150 atomic mass units, and inertia moment from 1500 to 2500 m.u. Angstrom.

Also, the energies involved in the base-base interactions are quite different for A-T and C-G pairs: the H-bond pairing energy is of 7.0 kcal/mole in the first case and 16.8 kcal/mole in the second one; the opening energy is respectively of 4.0 and 7.5 kcal/mole, and the stacking energies in the case of a homogeneous sequence is of 5.4 kcal/mole for a sequence of A-T (or T-A) pairs, and of 8.3 kcal/mole for G-C (or C-G) pairs.

This calls immediately for a more realistic modellization of the DNA chain, taking into account the different characteristics of different bases. It is quite clear, however, that renouncing to the spatial homogeneities of the chain (and to any kind of spatial periodicity as well if we want to model realistic sequences) means that such more realistic models will not have proper solitary waves, i.e. solutions of the form $`\varphi (x,t)=\varphi (xvt)`$ or the discrete analogue.

Naturally, from the point of view of biological significance of solutions we do not need that these are exactly solitary waves propagating with no dispersion and/or change in shape: it is sufficient that they are sufficiently stable over sufficiently long space and time intervals.

Obviously an analytical investigation of a realistic model with a true base sequence is outside reach, so this should be studied numerically. As far as I know, such a numerical study has never been conducted – opposite to the detailed studies conducted by Dauxois and Peyrard for the related modelling of the dissociation process by a realistic model based on the Peyrard-Bishop one – and so we cannot draw any conclusion in this respect (I do not know of any apriori obstacle to such a numerical study; in particular, it seems that for the nonlinear structure the weak helicoidal interaction, which was essential at the level of dispersion relations, has very little influence, so that one could study a simple “planar” model ).

Partial informations can be extracted for special kind of base sequences (at least in the purine-pyrimidine approximation); thus, for example, alternating regions such as “TATA boxes” would select a zero soliton speed . This is an interesting result, as TATA boxes are indeed know to code the pausing sites for the transcription process.

However, for more general sequences, one has reasons to think that a more detailed Yakushevich-like model could actually be worse than the original simple one from the point of view of correspondence to experimental observations; this could be due to the fact that the structures we are interested in involve several base pairs, so that there could be a mechanism of “effective self-averaging” in the dynamics .

Another relevant limitation of Yakushevich-like models is that there seem to be no selection of the soliton speed (see however below): solitons can exist at all speed below the maximal one $`v_{\mathrm{max}}=\beta \delta ^2`$, and energy considerations will not lead to a real selection among low speeds. This is not really encouraging, as transcription speed is not enormously fluctuating: we would expect there is some mechanism to regulate the soliton speed. This would also be essential to the functional role which the soliton is supposed to have, as it is necessary in this picture that the soliton speed is compatible with the speed of the transcription by the RNAP.

One should also recall that, although most of the investigation of the Yakushevich model have been based on the continuum approximation, such an approximation is not justified: actually, as mentioned above, in the case of DNA denaturation Dauxois and Peyrard have pointed out the essential role played by discreteness of the chain.

Again, it seems that no thorough investigation of the Yakushevich model on a chain (even an homogeneous one) have been conducted, so that properly speaking we cannot be sure about the existence and characteristic of soliton-like excitations in this case.

Another relevant point to be recalled is that a soliton in a continuous homogeneous medium would move with zero energy barrer. Essentially, this can be seen realizing that the field Lagrangian will involve $`\varphi _t^2`$ and $`\varphi _x^2`$, and when $`\varphi =\varphi (xvt)`$ these two terms become $`\varphi _z^2`$ and $`v^2\varphi _z^2`$: thus a soliton field configuration with zero speed will have an energy $`E_0`$, and one with a speed $`v=\epsilon `$, an energy $`E^{}=(1+\epsilon ^2)E_0`$.

However, the discrete structure of the chain lead to an effective barrer (known as Peierls barrer) for the motion. The result is that excitations of sufficiently low energy can be pinned on-site precisely by the presence of this barrer and thus by the discrete structure. This suggests that in a discrete Yakushevich model only solitons of sufficient high energy and thus speed could actually move along the chain, passing the Peierls barrer.

It should be noticed that this could provide the required mechanism for soliton speed, imposing $`v>v_0`$; the significance of this would of course depend on the actual value of $`v_0`$ and how it compares with the observed speed of the transcription bubble.

Anotrher major limitation of the Yakushevich – and similar – models is precisely the fact of focusing on the rotational motion connected with the reading of the DNA sequence by RNAP: indeed, the rigidity (or softness) of the concerned degree of freedom is quite comparable to the one of “radial” motion as considered in the Peyrard-Bishop model; thus, it would be highly desirable to consider these two degrees of freedom together, and combine the Peyrard-Bishop and the Yakushevich models into a single one.

Finally, the models – such as Yakushevich’s one – studied in connection with the problem of DNA transcription are all purely deterministic and only concern the DNA chain itself, i.e. with no interaction with its environment. One hopes that in a large part the chemico-physical parameters describing the environment could somehow be embodied in effective values of the constants defining the model (or one could limit to consider “standard” environment), but however the DNA has to work at finite, actually quite different from zero, temperature: thus unavoidably one has to consider thermal energy of the environment and energy exchanges involving this.

In considering this problem, one could resort to Statistical Mechanics considerations (which lie outside the scope of the present short report); or more in the spirit of Yakushevich-like models one should introduce random forces modelling thermal effects. It is not clear, nor – as far as I am aware – it has been studied if the solitons would survive the introduction of such random forces.

On the other side, the role of thermal energy is not necessarily negative from the point of view of existence of solitons: it is well possible that these same random forces would help the initial formation of nonlinear sructures by pumping energy into the more easily excitable modes until they ignite the formation of soliton-like collective structures.

It is maybe worth mentioning that statistical approaches to inhomogeneous models have been also considered recently; I will just mention the work of Hisakado and Wadati , where the distribution of masses and physical constants of the bases is a gaussian white noise, and the more specific work of Homma , which sets the basis for a statistical mechanics treatment of double chain models with sine-like nonlinearities, i.e. of the kind we consider here.

In this respect, it should be recalled that the whole process of transcription is strongly dependent on temperature (as is also that of denaturation), and in in vitro experiments where temperature is varied the transcription bubble forms sharply in a very small temperature range, see  and references therein: thus we are compelled to introduce dependence upon temperature in any modellization of the process.

This problem is considered e.g. in ; see also  for consideration of the environment in the somehow related problem of a single polyethilene chain.

5. Summary and conclusions.

I have given a short discussion of the idea that nonlinear excitations could play a role in the process of DNA transcription, i.e. that the transcription bubble could correspond to a solitary wave travelling along the chain, which the RNAP could then “surf” in order to access the base sequence with no energy to provide for opening the double helix.

I have discussed at some length the general idea of providing a simple model for a specific DNA process, and argued that despite the tremendous complexity of the DNA model this approach is not bound to fail. I have then recalled the main features of the model proposed by Yakushevich, mentioning some encouraging achievements and several limitations.

These limitations, however, more than being inherent to the model are limitation of the studies conducted so far. It is clear that the model is too simple to be valid as it is, and that it is needed to go “one step further” in the Yakushevich classification of DNA models , but only a more thorough analysis can focus on the detailed refinements which are needed.

In particular, I have pointed out several directions in which I believe it is necessary to generalize the model and to investigate its behaviour, such as considering real base sequences and thermal effects.

I hope this conference can, among other things, also stimulate a reflection on this theme, and maybe new works to analyze more realistic models of the nonlinear dynamics involved in DNA transcription.

References.

1. F. Calogero and A. Degasperis, Spectral transforms and solitons, North Holland (Amsterdam), 1982

2. P.L. Christiansen and V. Muto, Physica D 83 (1993), 93

3. Th. Dauxois, Phys. Lett. A 159 (1991), 390

4. Th. Dauxois and M. Peyrard, Phys. Rev. Lett. 70 (1993), 3935

5. Th. Dauxois, M. Peyrard and A.R. Bishop, Phys. Rev. E 47 (1993), R44 and 684

6. Th. Dauxois, M. Peyrard and C.R. Willis, Physica D 57 (1992), 267

7. A.S. Davydov, J. Theor. Biol. 38 (1973), 559

8. A.S. Davydov, Biology and Quantum Mechanics, Pergamon (London), 1982

9. P.G. Drazin and R.S. Johnson, Solitons: an introduction, Cambridge University Press (Cambridge), 1989, 1993

10. A. Dubrovin, S.P. Novikov and A. Fomenko, Modern Geometry, Springer (Berlin), 1984

11. S.W. Englander, N.R. Kallenbach, A.J. Heeger, J.A. Krumhansl and S. Litwin, Proc. Natl. Acad. Sci. U.S.A. 777 (1980), 7222

12. G. Gaeta, Phys. Lett. A 143 (1990), 227

13. G. Gaeta, Phys. Lett. A 168 (1992), 383

14. G. Gaeta, Phys. Lett. A 172 (1993), 365

15. G. Gaeta, Phys. Lett. A 179 (1993), 167

16. G. Gaeta, Phys. Lett. A 190 (1994), 301

17. G. Gaeta, C. Reiss, M. Peyrard and Th. Dauxois, Rivista del Nuovo Cimento 17 (1994), n.4

18. J.A. Gonzalez and M. Martin Landrove, Phys. Lett. A 191 (1994), 409

19. M. Hisakado and M. Wadati, J. Phys. Soc. Jap. 64 (1995), 1098

20. S. Homma, Physica D 113 (1998), 202

21. J.A. Krumhansl and J.R. Schrieffer, Phys. Rev. B 11 (1975), 3535

22. L.D. Landau and E.M. Lifsitz, Statistical Physics, Pergamon (London), 1955

23. V. Lisy and V.K. Fedyanin, J. Biol. Phys. 18 (1991), 127

24. S.K. Ma, Statistical Mechanics, World Scientific (Singapore), 1985

25. P. Manneville, Dissipative structures and weak turbulence, Academic Press (N.Y.), 1990

26. V. Muto, P.S. Lomdahl and P.L. Christiansen, Phys. Lett. A 136 (1989), 33; Phys. Rev. A 42 (1990), 7452

27. V. Muto, A.C. Scott and P.L. Christiansen, Physica D 44 (1990), 75

28. M. Peyrard ed., Nonlinear excitations in biomolecules, Springer (Berlin), 1995

29. M. Peyrard and A.R. Bishop, Phys. Rev. Lett. 62 (1989), 2755

30. M. Peyrard, Th. Dauxois, H. Hoyet and C.R. Willis, Physica D 68 (1993), 104

31. E.W. Prohofsky, Phys. Rev. A 38 (1988), 1538 and 5332

32. M. Remoissenet, Waves called solitons, Springer (Berlin), 1994

33. W. Saenger, Principles of nucleid acid structure, Springer (Berlin) 1984, 1988

34. A.C. Scott, Phys. Rep. 217 (1992), 1

35. A.C. Scott, Solitary waves in Biology, in

36. E.B. Starikov, Mol. Biol. 24 (1990), 1194

37. S. Takeno and S. Homma, Progr. Theor. Phys. 70 (1983), 308; 72 (1984), 679; 77 (1987), 548

38. M. Techera, L.L. Daemen and E.W. Prohofsky, Phys. Rev. A 40 (1989), 6636; 42 (1990), 1008

39. L.V. Yakushevich, Studia Biophys. 121 (1987), 201

40. L.V. Yakushevich, Phys. Lett. 136 (1989), 413

41. L.V. Yakushevich, Quart. Rev. Biophys. 26 (1993), 201

42. L.V. Yakushevich, Physica D 79 (1994), 77

43. S. Yomosa, J. Phys. Soc. Jap. 51 (1982), 3318; Phys. Rev. A 27 (1983), 2120; 30 (1984), 474

44. L.L. van Zandt, Phys. Rev. A 40 (1989), 6134; 42 (1990), 5036

45. F. Zhang and M.A. Collins, Phys. Rev. E 49 (1994), 5804

46. F. Zhang and M.A. Collins, Phys. Rev. E 52 (1995), 4217


format: text

text:
# The Θ⁺ (1540) as a heptaquark with the overlap of a pion, a kaon and a nucleon

## Abstract

We study the very recently discovered $`\mathrm{\Theta }^+`$ (1540) at SPring-8, at ITEP and at CLAS-Thomas Jefferson Lab. We apply the same RGM techniques that already explained with success the repulsive hard core of nucleon-nucleon, kaon-nucleon exotic scattering, and the attractive hard core present in pion-nucleon and pion-pion non-exotic scattering. We find that the $`KN`$ repulsion excludes the $`\mathrm{\Theta }^+`$ as a $`KN`$ s-wave pentaquark. We explore the $`\mathrm{\Theta }^+`$ as heptaquark, equivalent to a $`N+\pi +K`$ borromean boundstate, with positive parity and total isospin $`I=0`$. We find that the kaon-nucleon repulsion is cancelled by the attraction existing both in the pion-nucleon and pion-kaon channels. Although we are not yet able to bind the total three body system, we find that the $`\mathrm{\Theta }^+`$ may still be a heptaquark state.

In this talk we study the exotic hadron $`\mathrm{\Theta }^+`$ update  (narrow hadron resonance of 1540 MeV decaying into a $`nK^+`$) very recently discovered at Spring-8 Spring-8 , and confirmed by ITEP ITEP  and by CLAS at the TJNL CEBAF . This is an extremely exciting state because it may be the first exotic hadron to be discovered, with quantum numbers that cannot be interpreted as a quark and an anti-quark meson or as a three quark baryon. Very recent studies propose that the $`\mathrm{\Theta }^+`$ is a pentaquark state recent . Exotic multiquarks are expected since the early works of Jaffe Jaffe ; Strottman , and some years ago Diakonov, Petrov and Polyakov Diakonov  applied skyrmions to a precise prediction of the $`\mathrm{\Theta }^+`$. The nature of this particle, and its isospin, parity Hyodo  and angular momentum are yet to be determined.

We start in this talk by reviewing the Quark Model (QM) and the Resonating Group Method (RGM) Wheeler , which are adequate to study states where several quarks overlap. We show that the QM, together with chiral symmetry, produces hard core hadron-hadron potentials, which can be either repulsive or attractive. First we apply the RGM to show Bicudo ; Bender  that the exotic $`NK`$ hard core s-wave interaction is repulsive. This is consistent with the experimental data Barnes , see Fig. 1. We think that this excludes the $`\mathrm{\Theta }^+`$ as a bare s-wave pentaquark $`uddu\overline{s}`$ state or a tightly bound s-wave $`NK`$ narrow resonance. The observed mass of the $`\mathrm{\Theta }^+`$ is larger than the sum of the $`K`$ and $`N`$ masses by $`1540939493=118MeV`$, and this does not suggest a simple $`KN`$ binding. However a $`\pi `$ could also be present in this system, in which case the binding energy would be of the order of $`20MeV`$. Moreover this state of seven quarks would have a positive parity, and would have to decay to a p-wave $`NK`$ system, which is suppressed by angular momentum, thus explaining the narrow width of the $`\mathrm{\Theta }^+`$. With this natural description in mind we then apply the RGM to show that the $`\pi N`$ and $`\pi K`$ hard core interactions are attractive. Finally we put together the $`\pi N`$, $`\pi K`$ and $`NK`$ interactions to show that the $`\mathrm{\Theta }^+`$ is possibly a borromean borromean  three body s-wave bound state of a $`\pi `$, a $`N`$ and a $`K`$, with positive parity and total isospin $`I=0`$.

The RGM was used by Ribeiro Ribeiro  to show that in exotic hadron-hadron scattering, the quark-quark potential together with Pauli repulsion of quarks produces a repulsive short range interaction. For instance this explains the $`NN`$ hard core repulsion, preventing nuclear collapse. Deus and Ribeiro Deus  used the same RGM to show that, in non-exotic channels, the quark-antiquark annihilation could produce a short core attraction. The RGM computes the effective hadron-hadron interaction using the matrix elements of the microscopic quark-quark interactions. The wave functions of quarks are arranged in anti-symmetrized overlaps of simple colour singlet hadrons.

When spontaneous chiral symmetry breaking is included in the quark model Bicudo3 ; Bicudo4  the annihilation potential becomes crucial Bicudo1 ; Bicudo2  to understand the low mass of the $`\pi `$. The annihilation potential $`A`$ is also present in the $`\pi `$ Salpeter equation where it cancels most of the kinetic energy and confining potential $`2T+V`$. From the quark model with chiral symmetry breaking we get the matrix element Bicudo3 ,

$`A_{S=0}`$ $`=`$ $`{\displaystyle \frac{2}{3}}(2M_NM_\mathrm{\Delta }).`$ (1)
While the standard quark-quark potential provides a repulsive (positive) overlap, it is clear in eq. (1) that the annihilation potential provides an attractive (negative) overlap. This confirms that the hard core can be attractive for non-exotic channels where annihilation occurs, while it is repulsive for exotic channels.

We summarize Bicudo ; Bicudo1 ; Bicudo2  the effective potentials for the different channels,

$`V_{KN}`$ $`=`$ $`{\displaystyle \frac{2\frac{4}{3}\stackrel{}{\tau }_A\stackrel{}{\tau }_B}{\frac{5}{4}+\frac{1}{3}\stackrel{}{\tau }_A\stackrel{}{\tau }_B}}{\displaystyle \frac{(M_\mathrm{\Delta }M_N)}{3}}\left({\displaystyle \frac{2\sqrt{\pi }}{\alpha }}\right)^3e^{\frac{p_{\lambda }^{}{}_{}{}^{2}}{2\beta ^2}}{\displaystyle \frac{d^3p_\lambda ^{}}{(2\pi )^3}e^{\frac{p_{\lambda }^{}{}_{}{}^{2}}{2\beta ^2}}}`$
$`V_{\pi N}`$ $`=`$ $`{\displaystyle \frac{2}{9}}(2M_NM_\mathrm{\Delta })\stackrel{}{\tau }_A\stackrel{}{\tau }_B𝒩_{\alpha }^{}{}_{}{}^{2},`$
$`V_{\pi K}`$ $`=`$ $`{\displaystyle \frac{8}{27}}(2M_NM_\mathrm{\Delta })\stackrel{}{\tau }_A\stackrel{}{\tau }_B𝒩_{\alpha }^{}{}_{}{}^{2},`$ (2)
where $`\stackrel{}{\tau }`$ are the isospin matrices. This parametrization in a separable potential enables us to use standard techniques Bicudo  to exactly compute the scattering $`T`$ matrix. The scattering length $`d\delta _0/dk_{c.m.}`$ is,

$$a=\frac{\sqrt{\pi }}{\alpha }\frac{4\mu v}{\alpha ^2+\frac{\beta }{\alpha }4\mu v}.$$
(3)
The parameters and results for the relevant channels are summarized in Table 1. We have fitted $`\alpha `$ with the $`I=1`$ $`KN`$ scattering. We use $`\beta =\alpha `$ in the repulsive channels and fit it with the appropriate scattering lengths in the attractive $`I=1/2`$ pionic channels.

Let us first apply our method to the $`KN`$ exotic system, where the anti-quark $`\overline{s}`$ is present. In this case the $`I=0`$ channel is less repulsive than the $`I=1`$ channel. With our method we reproduce the $`KN`$ exotic s-wave phase shifts – see Fig. 1 – where indeed there is no evidence for the $`\mathrm{\Theta }^+`$ state. In what concerns the $`\pi N`$ system and the $`\pi K`$ systems, the corresponding parameters in Table 1 are almost identical. There we find repulsion for $`I=3/2`$ and attraction for $`I=1/2`$. The repulsion in the $`I=3/2`$ channel prevents a bound state in this channel. In what concerns the $`I=1/2`$ channel, the attraction is not sufficient to provide for a bound state, because the $`\pi `$ is quite light and the attractive potential is narrow. From eq. (3) we conclude that we have binding if the reduced mass is $`\mu \frac{\alpha ^3}{4\beta v}`$. With the present parameters $`\mu `$ should be of the order of 1 Fm-1. This is larger than the $`\pi `$ mass, therefore it is not possible to bind the $`\pi `$ to the $`K`$ or to the $`N`$. All that we can get is a very broad resonance. For instance in the $`\pi K`$ channel this is the kappa resonance Rupp1 , which has been recently confirmed by the scientific community. However, with a doubling of the interaction, produced by a $`K`$ and a $`N`$, we expect the $`\pi `$ to bind.

We now investigate the borromean borromean  binding of the exotic $`\mathrm{\Theta }^+`$ constituted by a $`N`$, $`K`$ and $`\pi `$ triplet. In what concerns isospin, we need the $`\pi `$ to couple to both $`N`$ and $`K`$ in $`I=1/2`$ states for attraction, and the only candidate for binding is the total $`I=0`$ state, see Fig. 2.

Since the $`\pi `$ is much lighter that the $`NK`$ system, we study the borromean binding adiabatically. As a first step, we start by assuming that the $`K`$ and $`N`$ are essentially stopped and separated by $`\stackrel{}{r}_N\stackrel{}{r}_K=2a\widehat{e}_z`$. This will be improved later. We also take advantage of the similarities in Table 1 to assume that the two heavier partners of the $`\pi `$ have a similar mass of 3.64 Fm-1 and interact with the $`\pi `$ with the same separable potential. Then we solve the bound state equation for a $`\pi `$ in the potential $`V_{\pi N}+V_{\pi K}`$, where this potential is wider the direction of the $`KN`$ $`z`$ axis. The $`\pi `$ energy is depicted if Fig. 3, and it is negative as expected. Once the $`\pi `$ binding energy is determined, we include it in the potential energy of the $`KN`$ system, which becomes the sum of the repulsive $`KN`$ potential and the $`\pi `$ energy. We find that for short distances the total potential is indeed attractive. Finally, using this $`KN`$ potential energy we solve the schrödinger equation for the system, thus including the previously neglected $`K`$ and $`N`$ kinetic energies. However here we are not able to bind the $`KN`$ system, because the total effective $`KN`$ potential is not sufficiently attractive to cancel the positive $`KN`$ kinetic energy. Nevertheless the $`K`$ is heavier than the $`\pi `$, thus a small enhancement of the attraction would suffice to bind the heptaquark. We remark that existing examples of narrow resonances with a trapped $`K`$ are the $`f_0(980)`$ Rupp1 , the $`D_s(2320)`$ Rupp2  very recently discovered at BABAR Babar , and possibly the $`\mathrm{\Lambda }(1405)`$. Moreover we expect that meson exchange interactions and the irreducible three-body overlap of the three hadrons, that we did not include here, would further increase the attractive potential. Therefore it is plausible that a complete computation will eventually bind the $`K\pi N`$ system.

We conclude that the $`\mathrm{\Theta }^+`$ hadron very recently discovered cannot be an s-wave pentaquark. We also find that it may be a heptaquark state, composed by the overlap of a $`\pi `$, a $`K`$ and a nucleon. This scenario has many interesting features. The $`\mathrm{\Theta }^+`$ would be, so far, the only hadron with a trapped $`\pi `$. Moreover the $`\pi `$ would be trapped by a rare three body borromean effect. And the decay rate to a $`K`$ and a $`N`$ would be suppressed since the $`\pi `$ needs to be absorbed with a derivative coupling, while the involved hadrons have a very low momentum in this state. Because the $`\mathrm{\Theta }^+`$ would be composed by a $`N`$ and two pseudoscalar mesons, its parity would be positive, $`J^P=\frac{1}{2}^+`$, in agreement with the prediction of Diakonov, Petrov and Polyakov. The isospin of the $`\mathrm{\Theta }^+`$ would be $`I=0`$ in order to ensure the attraction of the $`\pi `$ both by the $`N`$ and the $`K`$.

We are very grateful to George Rupp for pointing our attention to the pentaquark state. The work of G. Marques is supported by Fundação para a Ciência e a Tecnologia under the grant SFRH/BD/984/2000.


text:
# Constructing Ensembles of Pseudo-Experiments

## I Introduction

Many statistical analyses in physics are based on a frequency interpretation of probability. For example, the result of measuring a physical constant $`\theta `$ can be reported in the form of a $`1\alpha `$ confidence interval $`[X_1,X_2]`$, with the understanding that if the measurement is replicated a large number of times, one will have $`X_1\theta X_2`$ in a fraction $`1\alpha `$ of the replications. This type of interpretation therefore requires the definition of a reference set of similar measurements:

> The reference set of a measurement is the ensemble of experiments in which the actually performed experiment is considered to be embedded for the purpose of interpreting its results in a frequentist framework.

A major appeal of frequentism among physicists is its empirical definition of probability. By the strong law of large numbers, probabilities can be approximated in finite ensembles, and such approximations converge to the true value as the ensemble size increases. In other words, frequentist confidence statements are experimentally verifiable.

Physicists use Monte Carlo generated ensembles in various applications: to check a fitting algorithm for the presence of bias, non-Gaussian pulls, or other pathologies; to calculate the coverage of confidence intervals or upper limits; to average out statistical fluctuations in order to isolate systematic effects; to calculate goodness-of-fit measures and significances; to design experiments; etc. When constructing ensembles to address these questions, one needs to pay attention to a number of subtle issues that arise in a frequentist framework: what is the correct stopping rule?; is it appropriate to condition, and if so, on what statistic?; how should nuisance parameters be handled?

The aim of this paper is to draw attention to these issues and to propose some recommendations where possible. We start by discussing basic frequentist principles in section II and illustrate them with an example of conditioning in section III. The importance of stopping rules is argued in section IV. Finally, some purely frequentist methods to handle nuisance parameters are described in section V.

## II Frequentist Principles

In order to deserve the label frequentist, a statistical procedure and its associated ensemble must satisfy two core principles, which we examine in the next two subsections.

### II.1 The Frequentist Guarantee

The first principle states the aims of frequentism:

> Frequentist Guarantee baya2003 :
> In repeated use of a statistical procedure, the long-run average actual accuracy should not be less than (and ideally should equal) the long-run average reported accuracy.

To clarify this principle, we return to the $`1\alpha `$ confidence interval procedure mentioned in the Introduction. Let $``$ be an ensemble of intervals obtained by applying this procedure many times on different, independent data. The actual accuracy of an interval in $``$ is $`1`$ or $`0`$: either the interval covers the true value of the parameter of interest, or it does not. The average actual accuracy is therefore simply the fraction of intervals in $``$ that cover. On the other hand, the average reported accuracy is $`1\alpha `$. The reported accuracy is often the same for all intervals in $``$, but in some settings it is possible to report a different, data-dependent accuracy for each interval. Thus, averaging the reported accuracy is not necessarily a trivial operation. A procedure that satisfies the Frequentist Guarantee is said to have coverage.

In a sense, the Frequentist Guarantee is only weakly constraining, because it does not require a procedure to have coverage when applied to repeated measurements of the same quantity. To see how this is relevant, consider the construction of a 68% confidence interval for the mean $`\mu `$ of a Poisson distribution. One procedure is to take all $`\mu `$ values satisfying $`(n\mu )^2/\mu 1`$, where $`n`$ is the observed number of events. The resulting interval actually undercovers for many values of $`\mu `$ and overcovers for other values, so that the Frequentist Guarantee appears to be satisfied on average. To make this statement more precise we need a weighting function with which to carry out the average over $`\mu `$. A simple proposal is to perform local smoothing of the coverage function, resulting in local average coverage baya2003 .

Physicists may object to this notion of local average coverage on the grounds that they sometimes repeatedly measure a given constant of nature and are then interested in the coverage obtained for that particular constant, not in an average coverage over “nearby” constants. A possible answer is that one rarely measures the quantity of interest directly. Rather, one measures a combination of the quantity of interest with calibration constants, efficiencies, sample sizes, etc., all of which vary from one measurement to the next, so that an effective averaging does take place.

Finally, it could be argued that even Bayesians should subscribe to some form of the Frequentist Guarantee. If, over repeated use, a 95% credible Bayesian interval fails to cover the true value more than 30% of the time (say), then there must be something seriously wrong with that interval.

### II.2 The Anticipation Criterion

Although the Frequentist Guarantee specifies how a statistical procedure should behave under many repetitions of a measurement, it does not indicate what constitutes a valid repetition, and hence a valid ensemble. To the extent that this question involves the notion of randomness, it is well beyond the scope of this paper. From a practical standpoint however, one would like to stipulate that all effects susceptible to interfere with that randomness must be recognized as such and included in the construction of the ensemble, i.e. “anticipated”brow1977 . Hence the second principle:

> Anticipation Criterion:
> Ensembles must anticipate all elements of chance and all elements of choice of the actual experiments they serve to interpret.

To clarify, “elements of chance” refers to statistical fluctuations of course, but also to systematic uncertainties when the latter come from nuisance parameters that are determined by auxiliary measurements. On the other hand, “elements of choice” refers to actions by experimenters, in particular how they decide to stop the experiment, and what decisions they make after stopping.

One can identify several levels of anticipation. At the highest level, the data collection and analysis methods, as well as the reference ensemble used to interpret results, are fully specified at the outset. They do not change once the data is observed. The reference ensemble is called “unconditional”.

At the second highest level, the data collection and analysis methods are fully specified at the outset, but the reference ensemble is not. The latter will be fully determined once the data is observed, and is therefore “conditional”. Although a conditional ensemble is not known before observing the data, it is a subset in a known partition of a known unconditional ensemble.

The lowest level of anticipation is occupied by Bayesian methods, which fully condition on the observed data. The reference ensemble collapses to a point and can therefore no longer be used as a reference.

As the level of anticipation decreases, the reference ensemble becomes smaller. A remarkable result is that within the second level of anticipation one can refine the conditioning partition to the point where it is possible to give a Bayesian interpretation to frequentist conclusions, and vice-versa jber1997 .

## III Conditioning

To illustrate the interplay between anticipation and conditioning, we present here a famous example originally due to Cox cox1958 . Suppose we make one observation of a rare particle and wish to estimate its mass $`\mu `$ from the momenta of its decay products. For the sake of simplicity, assume that the estimator $`X`$ of $`\mu `$ is normal with mean $`\mu `$ and variance $`\sigma ^2`$. There is a 50% chance that the particle decays hadronically, in which case $`\sigma =10`$; otherwise the particle decays leptonically and $`\sigma =1`$. Consider the following 68% confidence interval procedures:

1. Unconditional
If the particle decayed hadronically, report $`X\pm \delta _h`$, otherwise report $`X\pm \delta _{\mathrm{}}`$, where $`\delta _h`$ and $`\delta _{\mathrm{}}`$ are chosen so as to minimize the expected length $`\delta =\delta _h+\delta _{\mathrm{}}`$ subject to the constraint of 68% coverage. This yields $`\delta _{\mathrm{}}=2.20`$ and $`\delta _h=5.06`$. The expected length is 7.26.
2. Conditional
If we condition on the decay mode, then the best interval is $`X\pm 10`$ if the particle decayed hadronically, and $`X\pm 1`$ otherwise. So the expected length is 11.0 in this case.

Note that in both cases we used all the information available: the measurement $`X`$ as well as the decay mode. Both procedures are valid; the only difference between them is the reference frame. The unconditional ensemble includes both decay modes, whereas the conditional one only includes the observed decay mode.

The expected length is shorter for unconditional intervals than for conditional ones. Does this mean we should quote the former? If our aim is to report what we learned from the data we observed, then clearly we should report the conditional interval. Suppose indeed that we observed a hadronic decay. The unconditional interval width is then 10.12, compared to 20.0 for the conditional one. The reason the unconditional interval is shorter is that, if we could repeat the experiment, we might observe the particle decaying into the leptonic mode. However, this is irrelevant to the interpretation of the observation we actually made. This example illustrates a general feature of conditioning, that it usually increases expected length, and reduces power in test settings.

Another aspect of the previous example is that the conditioning statistic (the decay mode) is ancillary: its distribution does not depend on the parameter of interest (the particle mass). This is not always the case. Suppose for example that we are given a sample from a normal distribution with unit variance and unknown mean $`\theta `$, and that we wish to test $`H_0:\theta =1`$ versus $`H_1:\theta =+1`$. The standard symmetric Neyman-Pearson test based on the sample mean $`\overline{X}`$ as test statistic rejects $`H_0`$ if $`\overline{X}>0`$. It makes no distinction between $`\overline{X}=0.5`$ and $`\overline{X}=5`$, even though in the latter case we certainly feel more confident in our rejection of $`H_0`$. Although $`\overline{X}`$ is not ancillary, it is possible to use it to calculate a conditional “measure of confidence” to help characterize one’s decision regarding $`H_0`$ kief1977 . Unfortunately, a general theory for choosing such conditioning statistics does not exist.

## IV Stopping Rules

Stopping rules specify how an experiment is to be terminated. High-energy physics experiments are often sequential, so it is important to properly incorporate stopping rules in the construction of ensembles.

As a first example, consider the measurement of the branching fraction $`\theta `$ for the decay of a rare particle $`A`$ into a particle $`B`$. Suppose we observe a total of $`n=12`$ decays, $`x=9`$ of which are $`AB`$ transitions, and the rest, $`r=3`$, are $`A\to ̸B`$ transitions. We wish to test $`H_0:\theta =1/2`$ versus $`H_1:\theta >1/2`$.

A possible stopping rule is to stop the experiment after observing a total number of decays $`n`$. The probability mass function (pmf) is then binomial:

$$f(x;\theta )=\left(\begin{array}{c}n\\ x\end{array}\right)\theta ^x(1\theta )^{nx},$$
(1)
and the $`p`$ value for testing $`H_0`$ is:

$$p_b=\underset{i=9}{\overset{12}{}}\left(\begin{array}{c}12\\ i\end{array}\right)\theta ^i(1\theta )^{12i}=\mathrm{\hspace{0.33em}0.075}.$$
(2)
An equally valid stopping rule is to stop the experiment after observing a number $`r`$ of $`A\to ̸B`$ decays. Now the pmf is negative binomial:

$$f(x;\theta )=\left(\begin{array}{c}r+x1\\ x\end{array}\right)\theta ^x(1\theta )^r,$$
(3)
and the $`p`$ value is:

$$p_{nb}=\underset{i=9}{\overset{\mathrm{}}{}}\left(\begin{array}{c}2+i\\ i\end{array}\right)\theta ^i(1\theta )^3=\mathrm{\hspace{0.33em}0.0325}.$$
(4)
If we adopt a 5% threshold for accepting or rejecting $`H_0`$, we see that the binomial model leads to acceptance, whereas the negative binomial model leads to rejection.

Here is a more intriguing example jber1988 . Imagine a physicist working at some famous particle accelerator and developping a procedure to select collision events that contain a Higgs boson. Assume that the expected rate of background events accepted by this procedure is known very accurately. Applying his technique to a given dataset, the physicist observes $`68`$ events and expects a background of $`50`$. The (Poisson) probability for 50 to fluctuate up to 68 or more is 0.89%, and the physicist concludes that there is significant evidence against $`H_0`$, the background-only hypothesis, at the 1% level.

Is this conclusion correct? Perhaps the physicist just decided to take a single sample. But what would he have done if this sample had not yielded a significant result? Perhaps he would have taken another sample! So the real procedure the physicist was considering is actually of the form:

* Take a data sample, count the number $`n_1`$ of Higgs candidates, and calculate the expected background $`b`$;
* If IP$`(Nn_1|b)\alpha `$ then stop and reject $`H_0`$;
* Otherwise, take a second sample with the same expected background, count the number $`n_2`$ of Higgs candidates and reject $`H_0`$ if IP$`(Nn_1+n_2|\mathrm{\hspace{0.17em}2}b)\alpha `$.

For this test procedure to have a level of 1%, $`\alpha `$ must be set at 0.67%. Since the actual data had a $`p`$ value of 0.89%, the physicist should not have rejected $`H_0`$.

So now the physicist finds himself forced to take another sample. There are two interesting cases:

1. The second sample yields 57 candidate events, for a total of 125. The probability for the expected background (100 events now) to fluctuate up to 125 or more is $`0.88\%>0.67\%`$, so the result is not significant. However, the result would have been significant if the physicist had not stopped halfway through data taking to calculate the $`p`$ value!
2. The second sample yields 59 candidate events, for a total of 127. The $`p`$ value is now 0.52% and significance has been obtained, unless of course the physicist was planning to take a third sample in the event of no significance.

Bayesian methods are generally independent of the stopping rule. It is therefore somewhat ironic that frequentists, who start from an objective definition of probability, should end up with results that depend on the thought processes of the experimenter.

## V Nuisance Parameters

Most problems of inference involve nuisance parameters, i.e. uninteresting parameters that are incompletely known and therefore add to the overall uncertainty on the parameters of interest. To fix ideas, assume that we have a sample $`\{x_1,\mathrm{},x_n\}`$ whose probability density function (pdf) $`f(\stackrel{}{x};\mu ,\nu )`$ depends on a parameter of interest $`\mu `$ and a nuisance parameter $`\nu `$, and that the latter can be determined from a separate sample $`\{y_1,\mathrm{},y_m\}`$ with pdf $`g(\stackrel{}{y};\nu )`$. Correct inference about $`\mu `$ must then be derived from the joint pdf

$$h(\stackrel{}{x},\stackrel{}{y};\mu ,\nu )f(\stackrel{}{x};\mu ,\nu )g(\stackrel{}{y};\nu ).$$
(5)
What is often done in practive however, is to first obtain a distribution $`\pi (\nu )`$ for $`\nu `$, usually by combining measurement results with a sensible guess for the form of $`\pi (\nu )`$. Inference about $`\mu `$ is then based on:

$$h^{}(\stackrel{}{x};\mu )f(\stackrel{}{x};\mu ,\nu )\pi (\nu )𝑑\nu .$$
(6)
Although this technique borrows elements from both Bayesian and frequentist methodologies, it really belongs to neither and is more properly referred to as a hybrid non-frequentist/non-Bayesian approach.

We illustrate the handling of nuisance parameters with a simple $`p`$ value calculation. Suppose that a search for a new particle ends with a sample of $`n_0=12`$ candidates over a separately measured background of $`\nu _0=5.7\pm 0.47`$, where we ignore the uncertainty on the standard error $`0.47`$. Let $`\mu `$ be the unknown expected number of new particles among the 12 candidates. We wish to test $`H_0:\mu =0`$ versus $`H_1:\mu >0`$.

A typical model for this problem consists of a Poisson density for the number of observed candidates and a Gaussian for the background measurement. Using equation (6) with a simple Monte Carlo integration routine, one obtains a $`p`$ value of $`1.6\%`$. For reference, when there is no uncertainty on $`\nu _0`$ the $`p`$ value is $`1.4\%`$.

While there are many purely frequentist approaches to the elimination of nuisance parameters, few of these have general applicability. Concentrating on the latter, we discuss the likelihood ratio and confidence interval methods in the next two subsections.

### V.1 Likelihood Ratio Method

The likelihood ratio statistic $`\lambda `$ is defined by:

$$\lambda =\frac{\begin{array}{c}sup(\mu ,\nu |n_0,\nu _0)\hfill \\ \mu =0\hfill \\ \nu 0\hfill \end{array}}{\begin{array}{c}sup(\mu ,\nu |n_0,\nu _0)\hfill \\ \mu 0\hfill \\ \nu 0\hfill \end{array}},$$
(7)
where, for $`\nu _0\mathrm{\Delta }\nu `$:

$$(\mu ,\nu |n_0,\nu _0)\frac{(\mu +\nu )^{n_0}}{n_0!}e^{\mu \nu }e^{\frac{1}{2}\left(\frac{\nu \nu _0}{\mathrm{\Delta }\nu }\right)^2}.$$

Simple calculus leads to:

$$\begin{array}{cccc}2\mathrm{ln}\lambda & =& 2\left(n_0\mathrm{ln}\frac{n_0}{\widehat{\nu }}+\widehat{\nu }n_0\right)+\left(\frac{\widehat{\nu }\nu _0}{\mathrm{\Delta }\nu }\right)^2\hfill & \text{if}n_0>\nu _0,\\ & =& 0\hfill & \text{if}n_0\nu _0,\\ \multicolumn{4}{c}{\text{with:}\widehat{\nu }=\frac{\nu _0\mathrm{\Delta }\nu ^2}{2}+\sqrt{\left(\frac{\nu _0\mathrm{\Delta }\nu ^2}{2}\right)^2+n_0\mathrm{\Delta }\nu ^2}.}\end{array}$$

Since $`\lambda `$ depends on $`n_0`$ and $`\nu _0`$, its distribution under $`H_0`$ depends on the true expected background $`\nu _t`$. A natural simplification is to examine the limit $`\nu _t\mathrm{}`$. Application of theorems describing the asymptotic behavior of $`2\mathrm{ln}\lambda `$ must take into account that for $`n_0<\nu _0`$ the analytical maximum of the likelihood lies outside the physical region $`\mu 0`$. The correct asymptotic result is that, under $`H_0`$, half a unit of probability is carried by the singleton $`\{2\mathrm{ln}\lambda =0\}`$, and the other half is distributed as a chisquared with one degree of freedom over $`0<2\mathrm{ln}\lambda <+\mathrm{}`$.

For our example the expected background is only 5.7 particles however, so one may wonder how close this is to the asymptotic limit. Here is an algorithm to check this. Choose a true number of background events $`\nu _t`$ and repeat the following three steps a large number of times:

1. Generate a Gaussian variate $`\nu _0`$ with mean $`\nu _t`$ and width $`\mathrm{\Delta }\nu `$;
2. Generate a Poisson variate $`n_0`$ with mean $`\nu _t`$;
3. Calculate $`\lambda `$ from the generated $`\nu _0`$ and $`n_0`$.

The $`p`$ value is then equal to the fraction of pseudo-experiments that yield a likelihood ratio $`\lambda `$ smaller than the $`\lambda _0`$ obtained from the observed data.

Note that this algorithm does not “smear” the true value of any parameter, in contrast with equation (6). The price for this is that the result depends on the choice of $`\nu _t`$. For $`\nu _t`$ varying from $`0.5`$ to $`50`$, the $`p`$ value ranges from $`0.48`$ to $`1.2\%`$. A general prescription for dealing with a $`p`$ value dependence on nuisance parameters is to use the so-called supremum $`p`$ value:

$$p_{sup}=\underset{\nu }{sup}\text{IP}(2\mathrm{ln}\lambda 2\mathrm{ln}\lambda _0|\mu ,\nu )|_{\mu =0}$$

From a frequentist point of view, the supremum $`p`$ value is valid, in the sense that:

$$\text{IP}(p_{sup}\alpha )\alpha ,\text{for each}\alpha [0,1],$$
(8)
regardless of the true value of the nuisance parameter. Although it is often difficult to calculate a supremum, in this case it turns out to equal the asymptotic limit to a good approximation. In our example $`2\mathrm{ln}\lambda _0=5.02`$ and corresponds to $`p_{sup}p_{\mathrm{}}=1.25\%`$.

As the attentive reader will have noticed, the $`p`$ value is smaller for $`\mathrm{\Delta }\nu =0.47`$ than for $`\mathrm{\Delta }\nu =0`$. This is a consequence of the discreteness of Poisson statistics; it does not violate inequality (8) because $`p_{sup}`$ actually overcovers a little when $`\mathrm{\Delta }\nu =0`$. To avoid the bias resulting from this overcoverage, the use of mid-$`p`$ values is sometimes advocated for the purpose of comparing or combining $`p`$ values berr1995 .

### V.2 Confidence Interval Method

The supremum $`p`$ value introduced in the previous section can be defined for any test statistic, although it will not always give useful results. If for example in our new particle search we take the total number $`n_0`$ of observed candidates as test statistic, the $`p`$ value will be 100% since the background $`\nu `$ is unbounded from above. A more satisfactory method proceeds as follows berg1994 ; silv1996 . First, construct a $`1\beta `$ confidence interval $`C_\beta `$ for the nuisance parameter $`\nu `$, then maximize the $`p`$ value over that interval, and finally correct the result for the fact that $`\beta 0`$:

$$p_\beta =\underset{\nu C_\beta }{sup}\text{IP}(Nn_0|\mu ,\nu )|_{\mu =0}+\beta .$$

It can be shown that this is also a valid $`p`$ value.

For the sake of illustration with our example, we consider three choices of $`\beta `$ and construct the corresponding $`1\beta `$ confidence intervals for $`\nu _t`$:

$$\begin{array}{cccc}1\beta =\hfill & 99.5\%:\hfill & C_{0.005}\hfill & =[4.38,\mathrm{\hspace{0.17em}7.02}]\hfill \\ 1\beta =\hfill & 99.9\%:\hfill & C_{0.001}\hfill & =[4.15,\mathrm{\hspace{0.17em}7.25}]\hfill \\ 1\beta =\hfill & 99.99\%:\hfill & C_{0.0001}\hfill & =[3.87,\mathrm{\hspace{0.17em}7.53}]\hfill \end{array}$$

To calculate the $`p`$ value, a good choice of statistic is the maximum likelihood estimator of the signal, i.e. $`\widehat{s}n_0\nu _0`$. Under $`H_0`$, the survivor function of $`\widehat{s}`$ is given by:

$$\text{IP}(S\widehat{s})=\underset{k=0}{\overset{\mathrm{}}{}}\frac{1+\text{erf}\left(\frac{k\nu _t\widehat{s}}{\sqrt{2}\mathrm{\Delta }\nu }\right)}{1+\text{erf}\left(\frac{\nu _t}{\sqrt{2}\mathrm{\Delta }\nu }\right)}\frac{\nu _t^k}{k!}e^{\nu _t}$$

We then find:

$$\begin{array}{ccccc}1\beta =\hfill & 99.5\%:\hfill & p_\beta =\mathrm{\hspace{0.33em}1.6}\%\hfill & +0.5\%\hfill & =\mathrm{\hspace{0.33em}2.1}\%\hfill \\ 1\beta =\hfill & 99.9\%:\hfill & p_\beta =\mathrm{\hspace{0.33em}1.7}\%\hfill & +0.1\%\hfill & =\mathrm{\hspace{0.33em}1.8}\%\hfill \\ 1\beta =\hfill & 99.99\%:\hfill & p_\beta =\mathrm{\hspace{0.33em}1.88}\%\hfill & +0.01\%\hfill & =\mathrm{\hspace{0.33em}1.89}\%\hfill \end{array}$$

An important point about the confidence interval method is that, in order to satisfy the Anticipation Criterion, the value of $`\beta `$ and the confidence set $`C_\beta `$ must be specified before looking at the data. Since $`p_\beta `$ is never smaller than $`\beta `$, the latter should be small. In particular, if $`p_\beta `$ is used in a level-$`\alpha `$ test, then $`\beta `$ must be smaller than $`\alpha `$ for the test to be useful.

## VI Summary

From the practical point of view of someone analyzing data, the most critical property of frequentist ensembles is their “anticipatoriness.” This requires that all the structural elements of an analysis (i.e. test sizes, interval procedures, bin boundaries, stopping rules, etc.) be in place before looking at the data. The only exception to this requirement occurs in situations where conditioning is both possible and appropriate. Even in that case, the conditioning partition itself must be specified beforehand.