Josep Sardanyés
email: josep.sardanes@upf.edu

Santiago F Elena

Error threshold in rna quasispecies models with complementation

Keywords: complementation, defective interfering particles, error threshold, quasispecies theory, RNA viruses, systems biology

A general assumption of quasispecies models of replicons dynamics is that the fitness of a genotype is entirely determined by its sequence. However, a more biologically plausible situation is that fitness depends on the proteins that catalyze metabolic reactions, including replication. In a stirred population of replicons, such as viruses replicating and accumulating within the same cell, the association between a given genome and the proteins it encodes is not tight as it can be replicated by proteins translated from other genomes. We have investigated how this complementation phenomenon affects the error threshold in simple quasispecies mean field models. We first studied a model in which the master and the mutant genomes code for wild-type and mutant replicases, respectively. We assume that the mutant replicase has a reduced activity and that that the wild-type replicase does not have increased affinity for the master genome. The whole pool of replicases can bind and replicate both genomes. We then analyze a different model considering a more extreme case of mutant genomes, the defective interfering particles (DIPs) described in many cases of viral infection. DIPs, with a higher replication rate owed to their shorter genomes, do not code for replicase, but they are able of using the replicase translated from the master genome. Our models allow to study how the probability of interaction between the genomes and the whole pool of replicases affects the error threshold. In both systems we characterize the scenario of coexistence between master and mutant genomes, providing the critical values of mutation rate, μ c, and the critical interaction rate between master genomes and replicases, γ c, at which the quasispecies enters into error catastrophe, a situation in which the mutant genomes dominate the population. In both cases, we showed that the error-threshold transition is given by transcritical-like bifurcations, suggesting a continuous phase transition. We have also found that the region in the parameter space (μ, γ) in which the master sequence survives is reduced when DIPs are introduced into the system.

I. INTRODUCTION

RNA viruses are the most important pathogens infecting plants and animals, and they are also a continuous source of emerging infectious diseases. The extremely short generation time, large population size and high mutation rates make RNA viruses excellent tools for experimental evolution and for testing basic principles of Evolutionary Theory [START_REF] Elena | Virus evolution: insights from an experimental approach[END_REF]. RNA viruses are characterized as highly polymorphic populations which are usually assimilated to quasispecies by many virologists [START_REF] Nowak | What is a quasispecies?[END_REF][START_REF] Domingo | RNA virus mutations and fitness for survival[END_REF][START_REF] Domingo | Quasispecies and RNA virus evolution: Principles and Consequences[END_REF]. The quasispecies population structure arises as a consequence of fast replication coupled with high mutation rates [START_REF] Domingo | RNA genetics[END_REF][START_REF] Domingo | Mutations and rapid evolution of RNA viruses[END_REF][START_REF] Domingo | Viruses at the edge of adaptation[END_REF][START_REF] Domingo | Nucleotide sequence heterogeneity of an RNA phage population[END_REF]. The mathematical quasispecies theory [START_REF] Eigen | Self-organization of matter and the evolution of biological macromolecules[END_REF][START_REF] Eigen | Molecular quasi-species[END_REF][START_REF] Eigen | The Molecular Quasispecies[END_REF][START_REF] Bull | Quasispecies Made Simple[END_REF] describes populations of replicons as a collection of a master genomes and a diverse cloud of mutants, which all together contribute to the phenotypic properties of the entire population. Under this view, the entire quasispecies, rather than individuals, is the target of natural selection [START_REF] Eigen | Viral quasispecies[END_REF][START_REF] Biebricher | The errorthreshold[END_REF]. A remarkable result of the quasispecies theory is the prediction of the existence of a phase transition at a critical mutation rate, μ c , the error-threshold, beyond which the master sequence is lost and the quasispecies is dominated by the mutant genomes. Roughly speaking, the critical mutation rate can be obtained as μ c ≈ ν -1 , being ν the sequence length. The theory predicts that replicons with mutation rates μ > μ c , might enter into the so-called error catastrophe regime, involving the out-competition of the master sequence by the pool of mutant genomes. Hence, genomic information is lost as the population enters into a drift phase [START_REF] Eigen | Self-organization of matter and the evolution of biological macromolecules[END_REF][START_REF] Schuster | How do RNA molecules and viruses explore their worlds?[END_REF].

The standard quasispecies mathematical model assumes that the fitness of any particular genotype entirely depends on its genome, without specifically distinguishing between effects at the genotypic and protein (phenotypic) levels. However, this assumption is highly unrealistic, because the fitness of a genome would depend on the functionality of the proteins it encodes. Furthermore, as master and mutant genomes coexist and both may contribute to the pool of proteins, it is conceivable that wild-type proteins may, for instance, act over mutant genomes, replicating and/or encapsidating them and FIG. 1: Schematic representation of the flow reactor approach used as abstract model to study quasispecies dynamics with complementation. The master (white) and the mutant (red) genomes encode for wild-type and mutant replicases, respectively. In Section II we consider the case where the mutant genomes code for deficient replicase with a reduced polymerization activity. In Section III we study the dynamics considering the mutant genomes as defective interfering particles (DIPs), which do not encode for a replicase (i.e., the process framed with the dotted line is removed). In both cases, the mutant genomes take advantage of the wild-type replicase to make more copies of themselves. Both systems are equivalent to a two-member hypercycle formed by a replicator and a parasite, which share the replicase (solid and dashed arrows indicate, respectively, replication and translation processes).

vice versa, in a clear case of functional trans complementation [START_REF] Wilke | Phenotypic mixing and hiding may contribute to memory in viral quasispecies[END_REF]. Henceforth, the phenotype of a genome may not reflect its genotype, a situation which in the case of viruses is known as phenotypic mixing and hiding [START_REF] Novick | Virus strains of identi-cal phenotype but different genotype[END_REF][START_REF] Brenner | Genetic control and phenotypic mixing of the adsorption cofactor requirement in bacteriophages T2 and T4[END_REF][START_REF] Huang | Pseudotype formation between enveloped RNA and DNA viruses[END_REF][START_REF] Holland | Virus mutation frequencies can be greatly underestimated by monoclonal antibody neutralization of virions[END_REF]. Moreover, a substantial amount of evidences, gathered with viruses as different as the murine cytomegalovirus [START_REF] Cicin-Sain | Frequent coinfection of cells explains functional in vivo complementation between cytomegalovirus variants in the multiply infected hosts[END_REF], foot-and-mouth disease virus [START_REF] García-Arriaza | Evolutionary transition towards defective RNAs that are infectious by complementation[END_REF] or tobacco mosaic virus [START_REF] Fraile | A quantitative analysis of complementation of deleterious mutants in plant virus populations[END_REF], give support to the fact that deletion mutants can be replicated in cells coinfected with the full-genome helper viruses.

Several mathematical models have been proposed that take into consideration trans interactions between different viral genomes with different fitness properties [START_REF] Gao | Complementation and epistasis in viral coinfection dynamics[END_REF].

For instance, studying the dynamics of defective interfering particles (DIPs) as an extreme case of complementation [START_REF] Bangham | Defective interfering particles: effects in modulating virus growth and persistence[END_REF][START_REF] Frank | Within-host spatial dynamics of viruses and defective interfering partiles[END_REF][START_REF] Kirkwood | Cycles, chaos, and evolution in virus cultures: A model of defective interfering particles[END_REF][START_REF] Szathmáry | Natural Selection and Dynamical Coexistence of Defective and Complementing Virus Segments[END_REF][START_REF] Szathmáry | Co-operation and Defection: Playing the Field in Virus Dynamics[END_REF]. DIPs are mutant viruses that lack most of the viral genome and cannot complete the infectious cycle by themselves [START_REF] Holland | Factors involved in the generation and replication of rhabdovirus defective T particles[END_REF][START_REF] Damayanti | Brome mosaic virus defective RNAs generated during infection of barley plants[END_REF] (see also [START_REF] Marriott | Defective interfering viruses and their potential as antiviral agents[END_REF] for a review). However, they can be replicated and encapsidated by the proteins translated from a helper virus coinfecting the same cell. The mechanisms ensuring the survival and persistence of DIPs are not entirely clear, as they behave as hyperparasites and usually get involve in an arms race with the full virus [START_REF] Horodyski | Properties of DI particles resistant mutants of vesicular stomatitis virus isolated from persistent infections and from undiluted passages[END_REF][START_REF] Depolo | Very rapid generation/amplification of defective interfering particles by vesicular stomatitis virus variants isolated from persistent infections[END_REF][START_REF] Depolo | Continuing coevolution of virus and defective interfering particles and of viral genome sequences during undiluted passages: virus mutants exhibiting nearly complete resistance to formerly dominant defective interfering particles[END_REF]. It is possible that the emergence of genomes shortened by deletions might confer an advantage in terms of replication speed compared with the full virus [START_REF] García-Arriaza | Evolutionary transition towards defective RNAs that are infectious by complementation[END_REF]. Furthermore, there is also evidence of a stronger form of interference whereby the DIPs genome competes more successfully for the viral replicative machinery [START_REF] Depolo | Very rapid generation/amplification of defective interfering particles by vesicular stomatitis virus variants isolated from persistent infections[END_REF][START_REF] Pattnaik | Cells that express all five proteins of vesicular stomatitis virus from cloned cDNAs support replication, assembly, and budding of defective interfering particles[END_REF].

As previously mentioned, several authors have theoretically investigated the dynamics of full viruses replicating together with DIPs. For example, Szathmáry [START_REF] Szathmáry | Natural Selection and Dynamical Coexistence of Defective and Complementing Virus Segments[END_REF] analyzed simple models based on mass action kinetics considering standard viruses and DIPs. In [START_REF] Szathmáry | Natural Selection and Dynamical Coexistence of Defective and Complementing Virus Segments[END_REF], structured deme models where developed to provide a description of the coexistence of virus segments considering standard virus and DIPs, sensitive and resistant viruses together with DIPs, covirus pairs (i.e., virus that exist as two or more separated particles all of which must be present for the complete replication cycle of a virus to occur), and virus-covirus systems. A deeper analysis of the model with standard virus and DIPs was later developed by Szathmáry [START_REF] Szathmáry | Co-operation and Defection: Playing the Field in Virus Dynamics[END_REF] considering cell populations infected by particles differing in number. Several different coexistence situations were shown to be possible by means of stable equilibria governed by fixed points or periodic orbits; also evidence was found of stable coexistence out of equilibrium governed by strange attractors.

Later, Kirkwood et al [START_REF] Kirkwood | Cycles, chaos, and evolution in virus cultures: A model of defective interfering particles[END_REF] developed a differential equations model to analyze a system formed by host cells, wild-type virions and DIPs. Such a model was able to explain several dynamic behaviors found in experiments, especially the fluctuations in virus titers on successive passages. Their model also gathered some interesting dynamical phenomena such as self-curing, which involves the extinction of the full virus together with the DIPs, and that had been previously observed to occur in vitro [START_REF] Jacobson | Determinants of spontaneous recovery and persistence in MDCK cells infected with lymphocytic choriomeningitis virus[END_REF]. They found that self-curing was associated to transient chaos, which might arise due to the presence of chaotic saddles in phase space.

More recently, Wilke et al [START_REF] Wilke | Phenotypic mixing and hiding may contribute to memory in viral quasispecies[END_REF], analyzed a simple model of complementation with differential fitness between master and mutant viruses due to an impaired ability to infect cells. The model studied by these authors predicted a strong influence of phenotypic mixing and hiding on population dynamics of viruses at high multiplicity of infection as well as important effects on the mutation-selection balance at low multiplicity of infection.

None of the aforementioned models, and to the ex- tend of our knowledge no one else, explores the effect of complementation on the error threshold predicted by the standard quasispecies theory. Therefore, neither the parametric regions in which the master sequence survives nor the bifurcations causing the error threshold have been characterized for systems of replicons in which complementation between mutant and master genomes may occur. We sought to cover this hole by studying two simple quasispecies models taking into account complementation during replication between master and mutant genomes.

The model studied in Section II considers that both master and mutant genomes encode, respectively, for wild-type and mutant replicases, with the mutant replicase having a lower catalytic activity. Hence, the fitness landscape is included at the replicase level and is given by the Swetina-Schuster fitness landscape [START_REF] Swetina | Self-replication with errors. A model for polynucleotide replication[END_REF]. However, since both types of replicases contribute to the cellular pool of proteins, both can replicate the master and mutant genomes. From a more biologically relevant perspective, the model can also give insights on the population dynamics of master and mutant viruses replicating within the same cell and may help to understand under which conditions master genomes can survive in the population.

The second model, developed in Section III, considers a problem of relevance in virology: the production of DIPs by some viruses and the effect they may exert on the full-length virus. In this model DIPs do not encode for replicase and have a faster replication rate and/or higher effective interaction with the wild-type replicase. Such a model will also provide some insights into the effect of competition between master genomes and DIPs, as well as into the possible persistence scenarios of the master sequence under the presence of DIPs. In both systems, replication is governed by a nonlinear, density-dependent growth. Actually, this system is equivalent to a two-member hypercycle with replicase-mediated replication (see Fig. 1). Some models of hypercycles coding for a replicase able to instruct the synthesis of other replicators can be found in [START_REF] Eigen | The Hypercycle. A Principle of natural Self-organization[END_REF].

II. QUASISPECIES MODEL WITH COMPLEMENTATION

We analyze the effect of complementation in the dynamics of RNA viruses with a simple quasispecies mean field model that describes the replication-mutation dynamics of RNA macromolecules in a flow reactor (see Fig. 1). Hence we study an unstructured model that only gathers the replication kinetics with a constant population, obviating the details of viral intracellular amplification. This is a standard method used in replicator theoretical models (see e.g., [START_REF] Tarazona | Error thresholds for molecular quasispecies as phase transitions: From simple landscapes to spin-glass models[END_REF][START_REF] Campos | Error propagation in the hypercycle[END_REF][START_REF] Silvestre | The information capacity of hypercycles[END_REF]). For mathematical convenience we studied the Swetina-Schuster single peak fitness landscape [START_REF] Swetina | Self-replication with errors. A model for polynucleotide replication[END_REF], where mutations are assumed to be largely deleterious and all mutant genotypes have identical low fitness. For simplicity we will not consider beneficial or neutral mutations. Our model considers the simplest scenario dividing the population of genomes in two types of sequences, the master and the mutant genomes, which are grouped into an "average" mutant sequence. The relative concentrations for such sequences are denoted x 0 ≥ 0 (master) and x 1 ≥ 0 (mutant), respectively. Such sequences are assumed to be positivesense strands that encode, respectively, for the wild-type (denoted as p 0 ≥ 0) and the mutant replicases (named p 1 ≥ 0). Hence, the master genomes are considered to be the sequences coding for the more efficient replicases. As a first approach we obviated the polarity of both sequences (i.e., genomic and antigenomic strands) not considering the effect of the intermediaries of replication.

A c c e p t e d m a n u s c r i p t

The model assumes a constant population (CP) constraint with 1 i=0 x i = 1 and

1 i=0 p i = 1.
Let us also assume that mutations occur from the master to the second compartment but not in the reverse sense. The enormous size of the sequence space makes this assumption a first good approximation. The model is then given by the following set of differential equations:

ẋ0 = (1 -μ)x 0 γ r 0 p 0 + r 1 p 1 -x 0 Φ(x), (1)
ẋ1 = μγx 0 + (1 -γ)x 1 r 0 p 0 + r 1 p 1 -x 1 Φ(x), (2) ṗi = k i x i -p i Ξ(p), (3)
with i = 0, 1. The outflux terms, which are computed from 1 i=0 ẋi = 0 and 1 i=0 ṗi = 0, are given by

Φ(x) = x 0 γ + (1 -γ)x 1 r 0 p 0 + r 1 p 1 ,
and Ξ(p) = 1 i=0 k i x i . Note that we use two different outfluxes because the RNA genomes and the replicase enzymes are competing for different resources (i.e., mononucleotides and amino acids, respectively). Here r 0 > 0 and r 1 > 0 are the replication rates due to the action of the wild-type and mutant replicases, respectively. We will consider that the mutant genome will always produce a mutant replicase whith a replication rate r 1 ≤ r 0 . Note that the term (r 0 p 0 + r 1 p 1 ) appears symmetrically in Eqs. ( 1)-( 2) because of complementation. Hence, both replicases, the master and the mutant, are able to bind and replicate both genome classes. The parameter μ, associated to the replication kinetics, is the average mutation rate. Finally, the parameter k i denote the translation rate of the replicases from both master and mutant genomes.

Finally, the parameter γ > 0 is the probability of interaction between the pool of replicases (i.e., wild-type and mutant replicases) and the master genomes, being (1-γ) the interaction probability between the mutant genomes and the pool of replicases. Our model allows to consider the effect of mutations at the level of the replicase and of the genomes, depending on the value of γ: with γ = 0.5 the replicase will interact with the same probability with the master and the mutant genomes. For this case, the deleterious effects of mutations are only incorporated at the level of the replicase. However, with γ > 0.5 the pool of replicases will undergo a better interaction with the master genomes, and the effect of mutations on the fitness will be also reflected at the level of the genotype. For the moment, we will not consider the case where the mutant genomes have more affinity for the pool of replicases (i.e., γ < 0.5).

The CP constraint allows reducing the previous model to a two-dimensional system because of the linear relations x 1 = 1 -x 0 and p 1 = 1 -p 0 . Hence, we can

:= (x 0 , p 0 ) ∈ R 2 : 0 ≤ x 0 , p 0 ≤ 1 . The reduced model is given by ẋ0 = (1 -μ)x 0 γ r 0 p 0 + r 1 (1 -p 0 ) -Φ(x)x 0 , ( 4
) ṗ0 = k 0 x 0 -p 0 k 0 x 0 + k 1 (1 -x 0 ) . ( 5
) Now with Φ(x) = x 0 γ + (1 -γ)(1 -x 0 ) r 0 p 0 + r 1 (1 -p 0 )
. This two-dimensional model has three equilibrium points. The first fixed point is the trivial one, P * 1 = (0, 0), which involves the extinction of the master genomes. Such a case would correspond to an absorbing state, in which the system remains trapped forever. The other two equilibria are

P * 2 = x * (2) 0 = 1 + γμ/(1 - 2γ), p * (2) 0 = k 0 (1+γ(μ-2))/(k 0 (1+γ(μ-2))-γk 1 μ) and
P * 3 = x * (3) 0 = k 1 r 1 /(k 1 r 1 -k 0 r 0 ), p * (3) 0 = r 1 /(r 1 -r 0 )
. We note that under the CP constraint, the fixed points for Eqs. (1-3) are fully determined and are given by: P * 1 = (0, 1, 0, 1),

P * 2 = x
, respectively. The Jacobian matrix for Eqs. (4-5) is given by

J = ⎛ ⎝ ∂ ẋ0 /∂x 0 ∂ ẋ0 /∂p 0 p 0 (k 1 -k 0 ) + k 0 x 0 (k 1 -k 0 ) -k 1 ⎞ ⎠ , with ∂ ẋ0 /∂x 0 = -(p 0 (r 0 -r 1 ) + r 1 )(1 -2x 0 + γ(-2 + μ + 4x 0 )) and ∂ ẋ0 /∂p 0 = -(r 0 -r 1 )x 0 (1 -x 0 + γ(-2 + μ + 2x 0 )).
In order to analytically find the critical mutation rate, μ c , responsible of the error-threshold, we study the stability of the trivial fixed point P * 1 , which, if stable, will involve the extinction of the master sequence. The Jacobian matrix evaluated at the trivial equilibrium takes the form

J(0) = ⎛ ⎝ r 1 (1 -μ)γ + (γ -1) 0 k 0 -k 1 ⎞ ⎠ .
The stability is studied from det |J(P * 1 ) -λI| = 0. We obtain two eigenvalues: λ (1) = (1 -μ)γr 1 + r 1 (γ -1) and λ (2) = -k 1 . Note that λ (2) is always negative, and thus the stability of this fixed point will entirely depend on λ (1) . The trivial fixed point, which involves the extinction of the master genomes, will be stable when μ > μ c or when γ < γ c , with μ c = 2 -γ -1 and γ c = (2 -μ) -1 , respectively. When γ = 1/2 the critical mutation rate is μ c = 0, and then, for any value of the mutation rate μ > 0, the master sequence is not able to survive. In Fig. 2 we display the dynamics obtained numerically from the full model with four variables. The critical boundary involving the error threshold and the extinction of the master genomes is shown in the phase diagram of Fig. 2(a). Here we computed the pairs of μ and γ where the master genome is able to persist for t ≤ 10 5 (white region). The solid line indicates the critical values of each of these parameters causing the error threshold phase transition. This curve perfectly fits with the analytical expressions of μ c and γ c previously derived from linear stability analysis. For this model, values of γ < 1/2 involve the extinction of the master genome independently of μ (results not shown). The dynamics with γ = 1/2 is illustrated in Fig. 2(b) by means of two-dimensional phase. For this case the master sequence extincts, being the entire phase plane (x 0 , p 0 ) a basin of attraction of the point (0, 0), which is globally stable. The exploration of the phase portraits using γ > 1/2 shows the non-trivial equilibrium of the master genome, and thus the coexistence of the master and the mutant genomes (see Figs.

2(c) and 2(d)).

If we now focus again in the reduced, two-dimensional model we can compute the values of the coordinates of the fixed points P * 2,3 at increasing μ and γ > 1/2. In Fig. 3(a) we show the values of such coordinates using μ as control parameter. We saw a linear decrease of the equilibrium concentration for the master genome, corresponding to the first coordinate of the fixed point P * 2 . This decrease was followed by the coordinate p * (2) 0 , because the production of wild-type replicase is proportional to the concentration of master genomes. The two coordinates of the equilibrium P * 3 remain negative for the parameter values used in Fig. 3. Similar results are shown in Fig. 3(b) using γ as control parameter. As γ decreases, the equilibria for the master genome and for the wild-type replicase decreased, and after the bifurcation value γ = γ c (with γ < γ c ) the master genomes are not able to persist.

The results shown in Figs. = (0, 0).

Hence, P * 2 = P * 1 for the critical values of μ and γ, indicating that such points collide and the bifurcation occurs. After the collision, the coordinate x * 0 of the fixed point P * 2 is negative, being out of Γ. Numerical computations show that after the collision, these two fixed points interchange the stability properties: with μ < μ c the fixed point P * 1 is unstable, while for μ > μ c becomes stable, as previously shown by the linear stability analysis. The stability of P * 2 works in an opposite way, that is, when μ < μ c , P * 2 is stable, and with μ > μ c is unstable. To better illustrate the effect of complementation in the behavior of the error threshold transition, we constructed the standard quasispecies model under the CP constraint (i.e., x 0 +x 1 = 1) and the Swetina-Schuster fitness landscape. The dynamical system, following similar notation than for the previous model, is given by

ẋ0 = f 0 (1 -μ)x 0 -x 0 η(x), ( 6
) ẋ1 = f 0 μx 0 + f 1 x 1 -x 1 η(x). ( 7
)
now with the outflow term given by η(x) = f i x i . Here f 0 and f 1 are the replication rates of the master (x 0 ) and the mutant (x 1 ) genomes, respectively. Equations ( 6)-( 7) can be reduced to a one-dimensional system, given by

ẋ0 = f 0 x 0 1 -μ - f 1 f 0 -x 0 1 - f 1 f 0 . ( 8
)
Actually, we can compute the time evolution of this single equation, which gives a logistic-like solution

x 0 (t) = ω 0 ω 1 1 + ω 0 /ω 1 -x 0 (0) x 0 (0) exp(-ω 0 f 0 t) -1
,

with ω 0 = 1 -μ -f 1 /f 0 and ω 1 = 1 -f 1 /f 0 .
The critical mutation rate for this simple quasispecies model without complemention, involving the extinction of the master genome is given by

μ c = 1 -f 1 /f 0 . ( 9
)
To compare the critical mutation rate with and without complementation we compute the equilibrium of the master genome by increasing mutation rate (indicated with a thick black line in Fig. 3a). Although both models with and without complementation are structurally different, we will use the same values of the replication rate studied in Fig. 3a. As expected, the results show that the critical mutation rate without complementation is higher. From expression (9), using f 0 = 0.9 and f 1 = 0.1, we obtain μ c = 0.888.... That is, considering complementation involves a lower mutation threshold because mutant genomes are taking advantage of the replicase translated from the master genomes and are, therefore, better competitors.

Finally, the equilibrium concentrations of the master and mutant genomes, computed numerically in the parameter space (μ, γ) using Eqs. (1-3), were in complete agreement with the previous results. We showed a decrease in the concentration of the master genomes at increasing mutation rate or decreasing γ, which is followed by an increase in the equilibria for the mutant genomes under such parameter changes (see Fig. 4).

In conclusion, we have shown that when complementation is incorporated into the picture, the critical mutation rate becomes independent of the replication rate of master and mutant genomes and only depends on the affinity of replicases for the master genome. From our results is

A c c e p t e d m a n u s c r i p t

also possible to conclude that complementation reduces the value of mutation rate at which transition to error catastrophe occurs.

III. QUASISPECIES MODEL WITH DIPS

Next, we focus our attention to a particularly interesting case in virology, namely the generation and amplification of defective interfering particles (DIPs) lacking the replicase gene. DIPs can only replicate in presence of a master genome that contributes with the wild-type replicase. This new model thus considers that the mutant genomes synthesized from the replication of the wild-type sequence produce a stronger interference because they are not producing mutant replicase which could contribute for replication of both genomes. As commented in the Introduction, experimental results suggest two different kinetic effects for DIPs: (i) the DIP, owed to its shorter genome length, replicates faster than the full-length master sequence and (ii) the DIP competes more successfully for interaction with the viral replicase machinery. Our model allows for the study of both effects: to simulate (i) we will set r 1 > r 0 , and by setting γ < 1/2, we will simulate (ii).

The DIPs model follows a similar notation than the one presented and analyzed in Section II. It is worth mentioning that, for this system with DIPs, the term master sequence should be redefined as the one coding for the replicase, even though this replicase synthesizes faster or has stronger affinity for the mutant genomes. The model is given by the following set of three nonlinear ODEs:

ẋ0 = (1 -μ)γr 0 x 0 p 0 -x 0 Φ(x), (10) ẋ1
= p 0 μγr 0 x 0 + (1 -γ)r 1 x 1 -x 1 Φ(x), (11) ṗ0
= k 0 x 0 -p 0 Ξ(p 0 ). (12)
Here we also use the CP constraint, with 1 i=0 x i = 1 and 0 ≤ p 0 ≤ 1. The outflux terms are now given by Φ(x) = p 0 (γr 0 x 0 + (1 -γ)r 1 x 1 ) , and Ξ(p 0 ) = k 0 x 0 . This model is close to the one analyzed by Szathmáry (see section 2.3 in [START_REF] Szathmáry | Natural Selection and Dynamical Coexistence of Defective and Complementing Virus Segments[END_REF]). This author analyzed the dynamics of standard and defective interfering viruses with a structured deme model. Although in such a work a model with mutation was also presented, the characterization of the critical parameters and the bifurcation involved in the error-threshold were not thoughtfully studied. As a difference, our model considers in an explicit way the product of the template (i.e., replicase) and differential probabilities of interaction between the DIP and the wild-type replicase. 10)-( 12) using μ and γ as control parameters. The equilibrium concentration for the master and mutant genomes are indicated with solid and dashed thick lines, respectively. The solid thin lines correspond to the population equilibria of the replicase. In (a) we use γ = 0.75 (red) and γ = 0.9 (black). In (b) we use μ = 0.05 (red) and μ = 0.2 (black). In both plots r 0 = 0.5, r1 = 0.9, k0 = 0.5 and initial conditions x0(0) = 1 and x 1(0) = p0(0) = 0. The arrows indicate the values of μc (a) and γc (b) for each analysis.

Let us analyze our model also reducing the system to a two-dimensional model using the linear relation obtained from the CP constraint, x 1 = 1 -x 0 . This system, also defined in the phase plane Γ (see Section II), is given by ẋ0 = x 0 p 0 Ψ(x 0 ),

= k 0 x 0 (1 -p 0 ). (13) ṗ0
with

Ψ(x 0 ) = (1 -μ)γr 0 -γr 0 x 0 + (1 -γ)(1 -x 0 )r 1 .
The points in the straight line {(0, p 0 )} are fixed points. For this case, the asymptotic state of the master sequence is the extinction. Hence, it can be interpreted as a ("onedimensional") absorbing state tied to the error catastrophe. This line of equilibria involves the extinction of the master genomes with the presence of wild-type replicase. This state is found because during the transient dynamics towards the extinction of the master genome, wild-type replicase molecules are still produced. As the wild-type replicases are not degraded, they remain in the system and, asymptotically, they reach a non-trivial value, although the master genomes become extincted. If (x 0 , p 0 ) is a fixed point with x 0 = 0, from Eq. ( 14) p 0 = 1, and then from Eq. ( 13) we obtain

x * 0 = (1 -μ)γr 0 -(1 -γ)r 1 γr 0 -(1 -γ)r 1 .
Then, assuming that 0 < x * 0 ≤ 1 we have another fixed point:

(x * 0 , 1) ≡ (1 -μ)γr 0 -(1 -γ)r 1 γr 0 -(1 -γ)r 1 , 1 .
The Jacobian matrix for this model is

J = ⎛ ⎝ p 0 Ψ(x 0 ) + x 0 Ψ (x 0 ) x 0 Ψ(x 0 ) k 0 (1 -p 0 ) -k 0 x 0 ⎞ ⎠ .
At the point (0, 0) both eigenvalues are zero (i.e., degenerate case). The linearized equation is unstable,

⎛ ⎝ u v⎞ ⎠ = ⎛ ⎝ 0 0 k 0 0 ⎞ ⎠ ⎛ ⎝ u v ⎞ ⎠ ,
u = 0 and v = k 0 u, and has solution u(t) = C 1 and v(t) = k 0 C 1 t + C 2 . These solutions with u(0) = C 1 = 0 go away from any neighborhood of (0, 0). At the points (0, p 0 ),

J = ⎛ ⎝ p 0 Ψ(0)) 0 k 0 (1 -p 0 ) 0 ⎞ ⎠ ,
and the two eigenvalues are λ 1 = 0 and

λ 2 = p 0 [(1 - μ)γr 0 -(1 -γ)r 1 ]
. Note that the points (0, p 0 ) will be the ones involved in the extinction of the master sequence for the model with DIPs. If

(1 -μ)γr 0 -(1 -γ)r 1 < 0, (15)
these points have a one-dimensional invariant stable manifold. From the previous inequation we can obtain the critical values responsible of the error threshold, which are given by

μ c = 1 -((1 -γ)r 1 /γr 0 )
and

γ c = r 1 /(r 0 (1 -μ) + r 1 ).
Then, if μ > μ c or γ < γ c the master genomes will disappear of the population, which will be completely dominated by the mutant DIPs. The previous critical values are in complete agreement with numerical simulations. In Fig. 5(a) we repeated the phase diagram shown in Fig. 2(a), where we represented (in gray) the pairs of the parameters (γ, μ) involving the extinction of the master genomes. Note that the critical boundary displayed with the solid black line perfectly fits with the values of γ c and μ c previously characterized. For the other fixed point (x * 0 , 1), Ψ(x * 0 ) = 0, and

J = ⎛ ⎝ x * 0 Ψ (x * 0 ) 0 0 -k 0 x * 0 ⎞ ⎠ .
The eigenvalues for this case are λ 1 = -kx * 0 < 0 and

λ 2 = x * 0 Ψ (x * 0 ), with x * 0 Ψ (x * 0 ) = -(1 -μ)γr 0 + (1 -γ)r 1 .
With this basic information obtained only from the linear study we can see that

(1 -μ)γr 0 -(1 -γ)r 1 , ( 16
)
plays an important role in the dynamics of the system. If expression ( 16) is negative, the points of the straight line {(0, p 0 )} are (weakly) attractive. However, with a positive value for expression [START_REF] Wilke | Phenotypic mixing and hiding may contribute to memory in viral quasispecies[END_REF] the point (x * 0 , 1) is attractive. The previous analytical results have been checked numerically. As previously discussed, the straight line {(0, p 0 )} contains fixed points, as it is shown in Fig. 5(b) in the scenario where {(0, p 0 )} has a stable invariant manifold (i.e., (1 -μ)γr 0 -(1 -γ)r 1 < 0). Note that different For these two cases μ < μ c and γ > γ c . The crossing of these critical values also involves transcritical-like bifurcations. This type of bifurcations, also found for the general model analyzed in Section II, is displayed in Fig. 6. The increase of mutation rate in the range μ < μ c involves a linear decrease in the concentration of the master genome and a linear increase of the mutant genomes concentration (Fig. 6(a)). Moreover, the decrease of γ involves a decrease in the equilibrium concentration of the master genomes, and when γ < γ c the population of master genomes is not able to survive at all (Fig. 6(b)). The previous results suggest that under the condition that the DIP has a higher affinity in the binding with the replicase, even for very small mutation rates (e.g., with μ = 0.05 as shown in Fig. 6(b) in red), these mutant DIPs will outcompete the master ones. Actually, the phase diagram of Fig. 5(a) shows that below γ < 0.5, the mutant DIPs outcompete the master one for any mutation rate (a similar result was also found in Section II).

We finally analyzed the effect of increasing the replication rate of the mutant genomes in the range r 0 < r 1 ≤ 1. The increasing of r 1 might correspond to the study of several DIPs with shorter genome lengths, assuming a linear relation between genome length and replicative speed for mutant genomes. These analyses were performed by computing numerically the equilibrium values of the master genomes in the parameter spaces (μ, r 1 ) (Fig. 7(a)) and (γ, r 1 ) (Fig. 7(b)). In both graphs it can be seen that the master genomes become extinct as the values of μ c and γ c are crossed. Moreover, as expected, we showed that an increase of the replication speed for the mutant DIPs involved an earlier extinction in both parameter spaces. That is, the values of μ c decreases as r 1 grows, and the values of γ c increase at increasing r 1 .

In conclusion, the previous results show that when DIPs are considered as a particular case of mutants, the critical mutation rate depends both on the replication rates and binding affinities. Furthermore, the more replicative advantage of the DIPs relative to the master genomes, the lower the critical mutation rate.

IV. DISCUSSION

The quasispecies theory [START_REF] Eigen | Self-organization of matter and the evolution of biological macromolecules[END_REF][START_REF] Eigen | Molecular quasi-species[END_REF][START_REF] Eigen | The Molecular Quasispecies[END_REF] has become the paradigmatic theoretical framework for the study of the dynamics and evolution of error-prone populations of replicators like RNA viruses. In the light of this theory, several interesting results have been reported, among them, the existence of the so-called error-threshold [START_REF] Eigen | Self-organization of matter and the evolution of biological macromolecules[END_REF][START_REF] Bull | Quasispecies Made Simple[END_REF][START_REF] Schuster | How do RNA molecules and viruses explore their worlds?[END_REF]. The error threshold has been defined as the value of the mutation rate below which populations equilibrate in a traditional mutation-selection balance and above which the population experiences an error catastrophe through excessive deleterious mutations [START_REF] Bull | Quasispecies Made Simple[END_REF]. The standard quasispecies model directly associates the fitness of a given replicator to its genotype. Translating this assumption to real viral populations may be cumbersome, since it has been amply demonstrated that the phenotype of a viral particle does not necessarily reflects the genotype it contains [START_REF] Novick | Virus strains of identi-cal phenotype but different genotype[END_REF][START_REF] Brenner | Genetic control and phenotypic mixing of the adsorption cofactor requirement in bacteriophages T2 and T4[END_REF][START_REF] Huang | Pseudotype formation between enveloped RNA and DNA viruses[END_REF][START_REF] Holland | Virus mutation frequencies can be greatly underestimated by monoclonal antibody neutralization of virions[END_REF].

Complementation between different viral genomes by means of trans interactions between different virus genotypes has been shown to be of importance in the dynamics of viral replication and cell infection, as shown in some theoretical studies [START_REF] Kirkwood | Cycles, chaos, and evolution in virus cultures: A model of defective interfering particles[END_REF][START_REF] Szathmáry | Natural Selection and Dynamical Coexistence of Defective and Complementing Virus Segments[END_REF][START_REF] Szathmáry | Co-operation and Defection: Playing the Field in Virus Dynamics[END_REF]. Previous modeling attempts of viral complementation, including interaction between DIPs and full viral genomes, focused in the ecological interactions [START_REF] Bangham | Defective interfering particles: effects in modulating virus growth and persistence[END_REF][START_REF] Frank | Within-host spatial dynamics of viruses and defective interfering partiles[END_REF][START_REF] Kirkwood | Cycles, chaos, and evolution in virus cultures: A model of defective interfering particles[END_REF] or population genetics [START_REF] Gao | Complementation and epistasis in viral coinfection dynamics[END_REF] of the different genotypes and the host cell. The main difference between these previous analyses and this work, is that here we studied two simple quasispecies models considering complementation during the replication of master (wild-type) and mutant genomes and the effect of such trans interactions on the error-threshold.

We first analyzed a model with master and mutant genomes coding for wild-type and mutant (functionally impaired) replicases. Due to complementation, both genomes can be replicated by both replicases, with differential interactions probabilities between master and mutant genomes with the whole pool of replicases (parametrized with γ) synthesized from both genome classes. Then, we analyzed a different model considering a particular case of complementation in which DIPs (see the review [START_REF] Marriott | Defective interfering viruses and their potential as antiviral agents[END_REF]) replicate faster due to their shortened genomes. Moreover, we can also incorporate into this model differential binding affinities for the unique, wild-type replicase.

For both models we characterized the coexistence scenario between master and mutant genomes, as well as the critical values of mutation and γ involving the entry into error catastrophe. That is, the values of such parameters with the full population dominated by the mutant genotypes. We found that the error-catastrophe scenario in the parameter space (μ, γ) was enlarged in the complementation model compared with the standard quasispecies model. Furthermore, we also found that the parameter space in which error-catastrophe occurs was even larger for the DIPs model than for the simple complementation one. We have also characterized the errorthreshold transition when crossing μ = μ c and γ = γ c as transcritical-like bifurcations in both systems under study, suggesting the presence of continuous phase transitions in the qualitative shift from viable master genomes to error catastrophe. Such a phase transition was characterized in previous works for the standard quasispecies model [START_REF] Tarazona | Error thresholds for molecular quasispecies as phase transitions: From simple landscapes to spin-glass models[END_REF]. However, other works with quasispecies models defined on sharp, single-peak landscapes have shown the presence of first-order phase transitions [START_REF] Pastor-Satorras | Field theory for a reaction-diffusion model of quasispecies dynamics[END_REF] as well as the possibility of both transitions occurring in the same system at increasing mutation rates [START_REF] Saakian | Phase diagram for the Eigen quasispecies theory with a truncated fitness landscape[END_REF]. As mentioned

A c c e p t e d m a n u s c r i p t 3 1 FIG. 2 :

 312 FIG.2: Dynamics generated by the model given by Eqs.[START_REF] Nowak | What is a quasispecies?[END_REF][START_REF] Domingo | RNA virus mutations and fitness for survival[END_REF][START_REF] Domingo | Quasispecies and RNA virus evolution: Principles and Consequences[END_REF]. (a) Phase diagram showing the region in the parameter space (γ, μ) in which the master sequence extincts (gray region, with parameters k 0 = k 1 = 0.4). The solid line indicates the errorthreshold boundary, which can be achieved by increasing μ or decreasing γ. This critical curve perfectly fits with the analytical values of μ c = 2 -γ -1 and γc = (2 -μ) -1 , obtained from linear stability analysis. We also display some phase portraits showing several trajectories (the arrows show the directions of the flows) on the phase plane (x0, p0) with: (b) k0 = k1 = 0.4 and γ = 1/2; (c) k 0 = k1 = 0.4 and γ = 0.75. In (d) we show the dynamics of the master and the mutant genomes in the simplex (x 0, x1), here with parameters k0 = 0.8, k1 = 0.4 and γ = 0.75. In all the plots we used r0 = 0.9, r1 = 0.1 and μ = 0.3. Thin trajectories correspond to the biologically meaningful initial conditions with 0 < x 0(0) ≤ 1 and x1(0) = p0(0) = p1(0) = 0. In (b) and (c) we used x1(0) = p1(0) = 0; in (d) we used p0(0) = p1(0) = 0. All the numerical results are obtained with the fourth-order Runge-Kutta method with a constant time step-size δt = 0.1.

FIG. 3 :

 3 FIG.3: Values of the coordinates of the fixed points, P * i (i = 1...3), for the two-dimensional model at increasing μ (a) and γ (b) with k 0 = 0.8, k1 = 0.4, r0 = 0.9 and r1 = 0.1 (in black we show the coordinate x * 0 and in red the coordinate p * 0 ). We indicate P * 2 and P * 3 with solid and dashed lines, respectively. Note that P * 3 remains negative in both panels. The intersections of the two lines (i.e., black and red) between fixed points are bifurcations. In (a) we use γ = 0.75 and in (b) we use μ = 0.3. The arrows indicate the stability of the extinction equilibrium P * 1 , which is indicated with a solid thin line. For comparison, we show the prediction of the error threshold for a simple quasispecies model without complementation, computing numerically the concentration of the master genomes (solid thick line) from Eqs. (6)-(7) using f 0 = 0.9, f1 = 0.1 and x0(0) = 1 and x1(0) = 0 as initial conditions. The critical values (e.g., μ c) involving the transition from viable to non-viable master genomes are indicated with a vertical dotted line.

FIG. 4 :

 4 FIG. 4: Equilibrium concentration of the master (a) and mutant (b) genomes numerically computed from the four-dimensional system (1-3) in the parameter space (μ, γ), using r 0 = 1, r1 = 0.1, k0 = 0.8 and k1 = 0.4. Initial conditions were: x0(0) = 1 and x1(0) = p0(0) = p1(0) = 0.

FIG. 5 :

 5 FIG.5: Same as in Fig.2but now considering the population dynamics with mutant genomes being DIPs and using Eqs. (10)-[START_REF] Eigen | Viral quasispecies[END_REF]. (a) Phase diagram with the critical boundary in the parameter space (γ, μ) involved in the error threshold (black line). This boundary perfectly matches with the critical parameter values analytically derived, and given by μ c = 1 -((1 -γ)r1/γr0) and γ c = r1/(r0(1 -μ) + r1). The phase portraits are displayed for: (b) μ = 0.4 and γ = 0.6; (c) and (d) μ = 0.2 and γ = 0.9. In all the plots we use r0 = 0.5, r1 = 0.9 and k0 = 0.5. In (b) and (c) x1(0) = 0 and in (d) p0(0) = 0.

FIG. 6 :

 6 FIG.6: Numerical values of the equilibria obtained from Eqs. (10)-(12) using μ and γ as control parameters. The equilibrium concentration for the master and mutant genomes are indicated with solid and dashed thick lines, respectively. The solid thin lines correspond to the population equilibria of the replicase. In (a) we use γ = 0.75 (red) and γ = 0.9 (black). In (b) we use μ = 0.05 (red) and μ = 0.2 (black). In both plots r 0 = 0.5, r1 = 0.9, k0 = 0.5 and initial conditions x0(0) = 1 and x 1(0) = p0(0) = 0. The arrows indicate the values of μc (a) and γc (b) for each analysis.

FIG. 7 :

 7 FIG. 7: Equilibrium concentration of the master genomes obtained numerically from Eqs. (10)-(12) in the parameter spaces (μ, r 1) (a), and (γ, r1) (b). In (a) we use γ = 0.75 and in (b) μ = 0.2. In both plots we set k0 = 0.5 and r0 = 0.5. Initial conditions were: x 0(0) = 1 and x1(0) = p0(0) = 0.

A c c e p t e d m a n u s c r i p t

  initial conditions achieve different asymptotic values in such a straight line, in which the master genomes become extinct. The other possible scenario, in which the master genomes survive and coexist with the DIPs, is shown in Figs. 5(c) and 5(d).

A c c e p t e d m a n u s c r i p t

with solid and dashed lines, respectively. Note that P * 3 remains negative in both panels. The intersections of the two lines (i.e., black and red) between fixed points are bifurcations. In (a) we use γ = 0.75 and in (b) we use μ = 0.3. The arrows indicate the stability of the extinction equilibrium P * 1 , which is indicated with a solid thin line. For comparison, we show the prediction of the error threshold for a simple quasispecies model without complementation, computing numerically the concentration of the master genomes (solid thick line) from Eqs. (6)-(7) using f 0 = 0.9, f1 = 0.1 and x0(0) = 1 and x1(0) = 0 as initial conditions. The critical values (e.g., μ c) involving the transition from viable to non-viable master genomes are indicated with a vertical dotted line.

Acknowledgments

The authors want to thank the useful comments and suggestions of Ernest Fontich. This work has been funded by the Human Frontier Science Program Organization grant RGP12/2008 and by the Spanish Ministerio de Ciencia e Innovación grant BFU2009-06993. We also acknowledge support from the Santa Fe Institute.

 A c c e p t e d ma n u s c r i p t 10 in the Introduction, our system might be equivalent to a hypercycle system with replicase-mediated replication.

The consideration of replicases, as we have shown, might involve the presence of a smoothed phase transition, which is in contrast with the governing bifurcations of ribozyme-like hypercycles, in which the extinction of the replicons is governed by saddle-node bifurcations [START_REF] Sardanyés | Error threshld ghosts in a simple hypercycle with error prone self-replication[END_REF][START_REF] Nuño | Influence of external fluctuations on a hypercycle formed by two kinetically indistinguishable species[END_REF] and by absorbing first-order phase transitions [START_REF] Sardanyés | Bifurcations and phase transitions in spatially-extended two-member hypercycles[END_REF]. Focusing in our first studied model, we have compared the critical mutation rate under complementation with the critical mutation obtained from a simple standard quasispecies model, also under the Swetina-Schuster fitness landscape. We have shown that the error threshold occurs at a lower mutation rate with complementation. We have also determined an error-threshold mediated by the binding probability of genomes with the replicases. Our results suggest that even for viruses that may be replicating far away from the mutation threshold (i.e., μ << μ c ), the population can be forced to error catastrophe by the existence of mutant RNA genomes for which the replicase may show an increased affinity.

For mathematical convenience, we have used the Swetina-Schuster single peak fitness landscape [START_REF] Swetina | Self-replication with errors. A model for polynucleotide replication[END_REF] in our models, despite its clear lack of biological realism. In a previous study, we have shown that different fitness landscapes, with variable degrees of epistasis, have substantially different effects on the error threshold [START_REF] Sardanyés | Replication Mode and Landscape Topology Differentially Affect RNA Virus Mutational Load and Robustness[END_REF]. Therefore, a logical follow up of the present study would be to explore the combined effect of epistasis and complementation among mutations on the properties of the error threshold. Gao and Feldman [START_REF] Gao | Complementation and epistasis in viral coinfection dynamics[END_REF] have recently developed a population genetic model of the effect of trans complementation with epistasis in viral dynamics. These authors found that the election of the epistasis function (i.e., synergistic or antagonistic) have profound effects on the genetic composition of complementing viral genomes. For instance, synergistic epistasis tends to strengthen the selection against the combinations of less-fit mutations, while frequent coinfection generally weakens the selection against all mutations, with the result that the master and single mutants constitute a major proportion of the virus population, and the stronger the complementation is, the lower the frequency of master genomes.

By contrast, antagonistic epistasis among deleterious mutations produce a much more complex picture and depends on whether antagonistic epistasis takes the form of magnitude or sign epistases [START_REF] Weinreich | Sign epistasis and genetic constraint on evolutionary trajectories[END_REF]. With sign antagonistic epistasis (e.g., if both the master and high-order mutants are fitter than the single mutants), master genomes tend to disappear and the population is dominated by genotypes with large numbers of mutations. However, with magnitude antagonistic epistasis (e.g., the strength of epistasis depends on the genetic background but the sign is unconditional), the master, single and double mutants dominate the population. These analytical results suggest that with complementation the error threshold may occur at lower critical mutation rates with sign antagonistic epistasis than with magnitude antagonistic epistasis and than with synergistic epistasis.


B M Jensen

P M H Damgaard

U Heegaard

C Halekoh

Lauridsen

Mrs Sine

DVM Nygaard Langerhuus

DVM Sine Nygaard Langerhuus
email: sine.langerhuus@agrsci.dk

MD, DMsc Else Kirstine Tønnesen

Professor; Karin Hjelholt Jensen
email: karinh.jensen@agrsci.dk

MSc In Biology

PhD ; Birthe

DVM M Damgaard
email: birthem.damgaard@agrsci.dk

PhD; Peter Heegaard

Professor; Ulrich Halekoh
email: ulrich.halekoh@agrsci.dk

MSc in Statistics, PhD; Charlotte Lauridsen
email: charlotte.lauridsen@agrsci.dk

MSc In Animal Nutrition

Marie Birthe

Damgaard

Peter Mikael

Helweg Heegaard

Biomarkers of aortic vascular prosthetic graft infection in a porcine model with Staphylococcus aureus

Keywords: biomarker, C-reactive protein (CRP), vascular prosthetic graft infection, porcine model, Staphylococcus aureus biomarker/C-reactive protein (CRP)/vascular prosthetic graft infection/porcine model/Staphylococcus aureus AVPGI: Aortic Vascular Prosthetic Graft Infection, CD+number: Cluster of differentiation + number, CFU: Colony Forming Units, CRP: C-Reactive Protein, ELISA: Enzyme Linked Immuno Sorbent Assay, LPS: Lipopolysaccharide, MFI: Mean Fluorescence Intensity, MHC II: Major Histocompatibility Complex II, NK cells: Natural Killer cells, SPF: Specific Pathogen Free S. aureus: Staphylococcus aureus, TNF-α: Tumor Necrosis Factoralpha, WBC: White Blood Cells

  L'archive ouverte pluridisciplinaire

Dear Alex van Belkum, Thank you very much for giving us the chance to revise the manuscript "Biomarkers of aortic vascular prosthetic graft infection in a porcine model with Staphylococcus aureus (Ref.: Ms. No. EJCMID-D-10-00212)". We hope the manuscript is reduced as much as required for a brief report, and that the reductions are performed so that it does not impair the readers' ability to reflect on the work. The list concerning the points being raised by the reviewer is given below. We hope the additional work concerning the rewriting of the manuscript will lead you to reconsider your previous decision, and publish the work in the European Journal of Clinical Microbiology and Infectious Diseases.

With kind regards, Sine Nygaard Langerhuus

Responses to reviewer's comments Point 1 and 4: We have fulfilled the reviewer's suggestion on changing the manuscript into a brief report and we hope that the compression of the manuscript is satisfactory and that no vital information to reflect on the study has been lost. We agree the specificity for CRP is low, and that this is caused by data from one animal. We believe that if more uninfected animals were included specificity would increase. Another way to increase specificity is to include the additionally evaluated biomarkers as a diagnostic set to rule out false positive animals, for the animal in question this would have been adequate.

Point 2:

We stopped using jugular catheters in this model after the reported study due to both difficulties in maintaining functionality as well as the potential problem of adding additional infection to the model. With the compression of the text to a brief report the catheter information is left out.

Point 3:

We also agree with the reviewer on the problem concerning the use of NSAID as analgesic in an inflammatory study. We are aware that NSAIDs may mask or suppress differences in inflammatory markers among groups in this study. The only reason for emphasising that all groups received the same analgesic treatment was that all groups then at least received the same inhibitory stimuli. Changes in inflammatory markers observed despite administered NSAIDs do have some importance, because NSAIDs are widely used analgesics both in human and veterinary medicine. We admit that in this study it would have been more elegant to leave them out, improving the chance of finding differences in evaluated biomarkers.

Discussion

Discussion has been incorporated in the result paragraph. All repetition of results in the discussion has been removed and much speculation has been removed or compressed.

Introduction

An 1-3% incidence of aortic vascular prosthetic graft infection (AVPGI) with a 30-day mortality of 6-13% has been reported [START_REF] Swain | Management of infected aortic prosthetic grafts[END_REF], and accountable for 43-64% of positive aortic graft cultures is S. aureus [START_REF] Mussa | Prevention and treatment of aortic graft infection[END_REF]. AVPGI can be diagnosed by imaging technology and surgical exploration. The limitations with imaging in the very early post-operative period are in distinguishing graft infection from post-operative changes in peri-graft tissues [START_REF] Swain | Management of infected aortic prosthetic grafts[END_REF].

Patients with vascular prosthetic graft infection show leukocytosis, elevated CRP, and positive bacteriologic blood cultures [START_REF] Mussa | Prevention and treatment of aortic graft infection[END_REF].

The purpose of the present study was to evaluate clinical signs and potential early biomarkers of S. aureus AVPGI in a porcine model.

It was hypothesized that plasma CRP, plasma fibrinogen, WBC, MHC II density and the lymphocyte CD4:CD8 ratio would increase, and TNF-α in vitro responsiveness after whole blood LipoPolySaccharide (LPS) stimulation would decrease in response to AVPGI.

Materials and Methods

The porcine model of AVPGI is described by Gao et al. [START_REF] Gao | Laparoscopic diagnosis and treatment of aortic vascular prosthetic graft infections in a porcine model[END_REF]. 16 Clinical signs of AVPGI, e.g. impaired general appearance, impaired movement, fever, and anorexia, were registered in the post-operative period. S. aureus growth was confirmed for the pigs in the LOW and HIGH groups at the end of the study period as described by Gao et al. [START_REF] Gao | Laparoscopic diagnosis and treatment of aortic vascular prosthetic graft infections in a porcine model[END_REF].

Blood was sampled prior to surgery and on day 2, 5, 7 and 14 post-operatively.

Plasma CRP was analyzed as described by Heegaard et al. [START_REF] Heegaard | A robust quantitative solid phase immunoassay for the acute phase protein C-reactive protein (CRP) based on cytidine 5'diphosphocholine coupled dendrimers[END_REF], and fibrinogen was determined according to standard procedures in an auto analyser, ADVIA 1650  Chemistry System (Bayer Corporation, Tarrytown, NY 10591, USA). Heparin stabilized whole blood (1mL) was stimulated with 50g LPS (Escherichia coli O111:B4, Sigma-Aldrich, Denmark) as described by Carstensen et al. [START_REF] Carstensen | Determination of tumor necrosis factor-alpha responsiveness in piglets around weaning using an ex vivo whole blood stimulation assay[END_REF] and TNF-α responsiveness was analyzed with a commercial porcine TNF/TNFSF2 Immunoassay (R&D Systems, Inc.).

The data concerning clinical signs of infection were analyzed as quasi-binomial distributed data and the biomarkers from blood samples were analyzed as normal distributed data. Generalized linear models were fitted to the data with the four level group effects as independent variable. The CRP and TNF-α data were log transformed prior to analysis, and repeated measurements were taken into account. Statistical significance was stated at the 5% level. Data were analyzed using the software R.

Results and discussion

No statistical difference in clinical signs of infection among the four interventions was observed. The raw data are given in Table 1.

As expected CRP was higher in the HIGH group compared to the SHAM and CLEAN group on day 2 and 5, whereas CRP in the LOW group did not differ significantly from the remaining groups. A CRP cut-off value above the estimated confidence interval for the CLEAN and SHAM group at 60µg/mL plasma was suggested for day 2 and 5 post-operatively. Plasma fibrinogen concentration however did not differ among groups, but a time effect was observed.

WBC was also higher in the HIGH group compared to the SHAM and LOW group on day 2 post-operatively. A cut-off value above the confidence intervals for WBC in the SHAM, CLEAN and LOW group at 25.5*10 9 /L blood on day 2 post-operatively was suggested.

MHC II density on peripheral leukocytes tended to be higher in the HIGH group compared to the CLEAN group on day 5 (P=0.06) and day 7 (P=0.06), respectively. This may indicate an up-regulation of MHC II due to S. aureus AVPGI, as was observed for human NK cells in response to S. aureus superantigens [START_REF] Orazio | Human natural killer (NK) cells present staphylococcal enterotoxin B (SEB) to T lymphocytes[END_REF]. Despite of lacking significance a cut-off value for MHC II density (MFI) above the confidence intervals for the SHAM and CLEAN group at 9900 MFI was suggested on day 5 and day 7 post-operatively.

The CD4:CD8 ratios did not however differ among groups, and no time effect was observed. CD8 is expressed on a broad range of porcine lymphocytes, including cytotoxic T-cells, NK-cells and CD4 positive lymphocytes, which may disturb the ratio interpretation as biomarker of the adaptive immune response. Furthermore, the TNF-α responsiveness was reduced in the HIGH group compared to the SHAM, CLEAN and LOW group on day 2 post-operatively. This was in agreement with a reported decreased TNF-α responsiveness in humans in relation to surgery and sepsis [START_REF] Ertel | Downregulation of proinflammatory cytokine release in whole blood from septic patients[END_REF][START_REF] Haupt | Monocyte function before and after surgical trauma[END_REF][START_REF] Kumpf | Pre-and post-operative cytokine release after in vitro whole blood lipopolysaccharide stimulation and frequent toll-like receptor 4 polymorphisms[END_REF]. However, the lack of effect in the LOW group contradicts this hypothesis. A cut-off value for TNF-α below the confidence interval of the estimated TNF-α responsiveness in the SHAM, CLEAN and LOW groups at 0.9µg/mL plasma was suggested on day 2 post-operatively. Data from the CLEAN, LOW and HIGH groups was used for sensitivity and specificity calculations using the suggested cut-off values as diagnostic criteria (Table 2). Using these criteria pigs in the CLEAN group were considered AVPGI negative, and pigs in the LOW and HIGH groups with a confirmed positive bacteriological evaluation was considered AVPGI positive.

Of all biomarkers evaluated CRP was superior for diagnosing S. aureus AVPGI in pig with a sensitivity of 0.86 and specificity of 0.75. Using three positive biomarkers out of the four (CRP, WBC, TNF-α in vitro responsiveness, and MHC II density) as diagnostic criterion increased the specificity but at the expense of the sensitivity.

In conclusion, AVPGI in this porcine model was not accompanied by clinical signs, but resulted in immunological changes comparable with the signs of AVPGI in humans. CRP was superior as an individual biomarker of AVPGI in this porcine model.

  pigs of (Landrace*Yorkshire)*Duroc genetic background weighing 62.4±6.2kg (s.d.) were assigned to one of four groups: SHAM (n=4): Pigs had anesthesia and laparotomi, and the infra renal part of aorta was exposed for the full duration-time of a graft insertion; CLEAN (n=4): Pigs had a graft inserted; LOW (n=4): Pigs had a graft inserted, and 5x10 4 CFU of S. aureus were inoculated on the graft; or HIGH (n=4): Pigs had a graft inserted, and 1x10 6 CFU of S. aureus were inoculated on the graft. The first four days post-operatively flunixin (Finadyne®Vet, 150mg IM 24-hour intervals) and buprenorfine (Temgesic®, 0,6mg IM with 8-hour intervals) were administered for analgesia. No antibacterial treatment was given post-operatively. The study was approved by the Danish Ministry of Justice, Animal Research Inspection.

  Differential counts were conducted by a Cell-Dyn R 3500R (Abbott, USA). Flow cytometry was conducted on a BD FACSCanto™. MHC II was stained with mouse anti SLA CLAS II DQ (MCA 1335, Abd serotec) and secondarily IgG:FITC (MCA Star 9B, Abd serotec). CD4 and CD8 were stained with mouse anti pig CD4a:FITC (MCA 1749F, Abd serotec) and mouse anti pig wCD8a:PE (MCA 1223PE, Abd serotec), respectively.

  Reference substitute brief AVPGI model description, Table 1 is removed Clinical signs of infection: Clinical signs evaluated are shortly mentioned and postoperative pain treatment and lack of antibacterial treatment are mentioned remaining text is removed Bacteriological evaluation: Reference substitute description Blood sampling: Reduced to contain sampling times CRP: Reference substitute description Fibrinogen: Compressed WBC: Compressed TNF-α in vivo responsiveness: Reference substitute stimulation description and the remaining text is compressed MHC II density and CD4:CD8 ratio: Compressed Statistical analysis: Compressed Results Subdivision is removed. Figure 1, 2, and 3 are removed leaving only Table 2 and 3 from the original manuscript. The results have been compressed excluding all model estimates and most P-values.

	The abstract has been compressed
	Introduction
	Paragraph 3 and 4 has been removed and the remaining text has been compressed
	Material and Methods
	Subdivision is removed
	Animals: Compressed
	Groups:

Acknowledgements

The surgical part of the study was conducted in collaboration with the Vascular Research Unit, Department of Vascular Surgery and Department of Clinical Microbiology, Viborg Hospital, 8800 Viborg, Denmark.

Table 1. Clinical signs of infection and bacteriological evaluations. The clinical results are given as the total number of days the individual pigs in the four groups (SHAM, CLEAN, LOW and HIGH) showed the following clinical signs of infection, decreased well-being, decreased movement, fever, and anorexia in the 14-day post-operative period. The results obtained from inoculating graft surface and peri-graft tissue swabs onto Colombia agar are given as confirmed bacterial growth (+) or disconfirmed bacterial growth (-)


M Kettlitz
email: kettlitz@inp-greifswald.de

J Rarbach
email: j.rarbach@vorschaltgeraete.de

J Zalach

Investigation of acoustic resonances in high power lamps

Keywords: numbers: 52.35.Dm, 52.77.Fv, 52.80.-s, 51.50.+v, 52.80.Mg light source, acoustic resonances, photo diode, high power lamp

High power, medium pressure, mercury containing lamps are used as UV sources for many industrial applications. Lamps investigated in this paper are driven with an electronic ballast with a non-sinusoidal current waveform at a fixed frequency of 20 kHz and a maximum power output of 35 kW. Instabilities can occur if the input power is reduced below 50%. The reason is identified as acoustic resonances in the lamp. Comparison of calculated and measured resonance frequencies shows a good agreement and explains the observed lamp behavior. This led to a development of a new ballast prototype which is able to avoid instabilities by changing the driving frequency in dependence on the applied power.

Introduction

Plasma light sources are widely used in industry as ultraviolet radiators for UV curing of ink, varnish and plastics, for water purification and for the interlacing of adhesives. The UV radiation initiates the chemical reactions for these processes, and a high radiation output is needed for an effective production process. In many cases, high-and mediumpressure mercury lamps are used for these purposes because of their good performance, high luminous flux and efficiency.

Generally, inductive ballasts are used to drive these lamps, but in the last years there was a change towards electronic ballasts. In order to keep them compact and small in size medium working frequencies are preferred. For both types of power supply, a good match of lamp and ballast is important to guarantee a stable and secure operation. The performance of plasma lamps can be reduced by several effects. One of them is an instability of the lamp plasma caused by acoustic resonances, which limits the parameter space for secure lamp operation. Instabilities are a potential safety risk and can destroy the lamp. Some less annoying but also unwanted effects which may occur are light flickering, color change and a reduced lifetime.

High power lamps are often designed to operate at one fixed power and frequency. Of course, it is desirable to extend the working parameters and, for instance, to reduce the lamp power without the risk of lamp destruction. Unfortunately, changing the lamp power will change temperature profiles in the plasma and therefore its properties. The transition between two applied electric power levels may pass the lamp's specific region of acoustic resonances which can lead to arc instability. There are several ways to omit instabilities which are already described in literature [START_REF] Bhosle | Proc. of the LS-10[END_REF][START_REF] Wada | [END_REF][3][START_REF] Bhosle | Fourtieth Industry Applications Conference Annual Meeting[END_REF][START_REF] Jongerius | [END_REF][6][7][8]. First, do not change the power, just operate the lamp at the fixed, stable power level. Second, change the lamp geometry to suppress the development of acoustic resonances, and third, adjust the electronic ballast to avoid the region of acoustic resonances during the change of operating parameters.

Actually the electric ballast has been identified as a very important factor to avoid acoustic resonances [START_REF] Bhosle | Proc. of the LS-10[END_REF]. Therefore electronic ballast engineering has become increasingly important for the development of new lamps and lamp systems. Other possibilities to avoid acoustic resonances are the variation of the lamp tube geometry as described in [START_REF] Wada | [END_REF].

In this work, acoustic resonances in high power lamps are investigated. With 35 kW lamp types this was done for the first time. In contrast to most other experiments the current waveform was non-sinusoidal with values of roughly 30 A at 35 kW. Another challenging aspect of this work was the combination of plasma diagnostics, namely optical emission spectroscopy and acoustic resonance diagnostics represented by the Fourier transform of the light fluctuation spectrum, which together delivered data for the computation of resonance frequencies.

The purpose of this work is not to develop a new strategy to avoid acoustic resonances. But to analyze already described methods of changing the driving frequency for its applicability on real ballast driving a particular type of high power lamp.

Acoustic Resonances

Acoustic resonances, like all instabilities, are problematic because they can lead to several unwanted effects. When established, the obvious symmetry of the straight, symmetric arc is broken and flickering, color change, acoustic noise, traveling arc structures and arc bending can be observed. In the last case the deflected arc can locally heat the wall, increasing its temperature and weaken its mechanical stability. This can result in lamp destruction due to deformation or a structural break of the glass body. Which of these effects occur and at which intensities depends on the lamp type used. In contrast to this, acoustic resonances can also exhibit advantageous aspects. By intentionally inducing well defined modes it is possible to reduce demixing effects and straighten the arc as demonstrated by Stockwald [START_REF] Stockwald | IEEE ICOPS[END_REF]. Furthermore, acoustic resonances can be used as a diagnostic tool to estimate the total gas pressure inside the discharge vessel [START_REF] Käning | Proc. of the LS-10[END_REF] The mechanism of acoustic waves occurring in lamps has already been described for example in [START_REF] Wada | [END_REF][3][START_REF] Bhosle | Fourtieth Industry Applications Conference Annual Meeting[END_REF][7][8][START_REF] Stockwald | IEEE ICOPS[END_REF][START_REF] De Groot | The high-pressure sodium lamp[END_REF] and is shortly summarized here.

Plasma light sources are often designed as glass or ceramic tubes. Gas filled cylinders define geometrical boundary conditions for pressure waves. An excitation mechanism is always included in form of alternating current at the electrodes. This leads to a periodical power input to the plasma.

Witting [3] computed the resonance frequencies analogously to standing pressure waves in cylinders (organ pipes) using the lamp length l, inner radius r and the speed of sound v s as parameters:

ω i = πi v s l (longitudinal mode) (1)
ω mn = a mn v s r (transversal modes) ( 2
)
ω 2 mni = ω 2 mn + ω 2 i (mixed resonance) (3) v s = γ RT m ( 4
)
The variables have the usual meaning -gas constant R, mass m, temperature T and coefficients a mn which are derived from the Bessel function by applying boundary conditions. For the numerical values of a mn Witting refers to Rayleigh [START_REF] Strutt | The Theory of Sound MacMillan[END_REF], they are also given in more recent publications, for example by Jahnke et al [START_REF] Jahnke | Tafeln höherer Funktionen Teubner Studienbücher[END_REF]. The term transversal modes means all wave modes propagating perpendicular to the arc. Spatial fluctuation pattern for several modes are graphically shown i.e. in [START_REF] Collatz | Differentialgleichungen Teubner Studienbücher[END_REF]. The speed of sound v s is usually calculated for a mean gas temperature with γ = 5 3 for monoatomic gases used as filling in lamps analyzed in this work. This model of standing pressure waves ignores all complex interactions of gas properties. Particularly, it takes not into account the radial temperature profile T (r) as it occurs in lamps. In this case an average temperature is considered. The temperature profile affects of course the speed of sound as well as other relevant plasma and gas properties, for example viscosity, convection, heat transfer, collision frequencies and electrical conductivity. This means that acoustic resonances are a complex phenomenon. Also, it can be seen in the fact that the observed hysteresis as a characteristic behavior can not be explained within the picture of standing pressure waves [START_REF] Bhosle | Fourtieth Industry Applications Conference Annual Meeting[END_REF].

This complexity has obvious consequences for lamp operation. An intentional change of the power input affects the temperature profile and therefore many plasma and gas properties along the discharge tube and along the inner radius which in turn change the temperature dependent resonance frequencies -the lamp can become unstable. These properties effectively limit the parameter region of safe lamp operation.

A more detailed look at the pulsating power input mentioned above leads to the electrodes. During a power cycle the electrodes become alternatively an anode or cathode -in dependence on the phase. The change between anode and cathode results in a periodical occurrence of a hot spot, which locally heats the gas. This corresponds to a recurring small 'gas explosion', which travels along the discharge tube as a pressure fluctuation. Usually, the electrode region is the origin of the instability as shown in figure 1. Geometrical boundary conditions given by the shape of the discharge volume define the tube specific resonance frequencies and give rise to the possibility of standing pressure waves. The pressure wave amplitude can grow if the frequency of the fluctuations is in phase with the frequency of the power input and lies close to a resonance frequency of the tube. The resonance frequencies are defined not only by the lamp geometry, but also by the temperature profile T (r) and gas type. Increasing amplitude of pressure waves can disturb the discharge path and its shape, the whole set of discharge parameters may change and strong instability with the above mentioned effects can occur.

The instability grows on time scales of 0.01-1 s, which is much longer than the period of power input. This means that the acoustic resonance generating conditions have to last for a relatively long time compared to the duration of one HF cycle to make growing of the pressure wave amplitude possible. Otherwise the arc can still be stable despite the presence of all ingredients for acoustic resonance [START_REF] Wada | [END_REF]6]. This property was used by Wada [START_REF] Wada | [END_REF] to suppress acoustic resonances by inducing a high frequency component to the driving voltage which distorted the growing process. A good example of the different time scales was given by Olsen [START_REF] Olsen | [END_REF], where a driving frequency of 60 Hz produced an acoustic resonance at 24 kHz. Furthermore, not all acoustic modes are in the same way sensitive to the excitation level, some need much stronger stimulation to develop than others [START_REF] Bhosle | Fourtieth Industry Applications Conference Annual Meeting[END_REF]. An example of a non-destructive acoustic resonance occurrence is shown in figure 2. The initial stage of a not fully developed instability can be seen in figure 1. An ideal electronic ballast should suppress the growth of acoustic resonances i.e. by changing the driving frequency in dependence on the applied power, in order to avoid the region of acoustic resonance frequencies.

Setup

The investigated medium pressure high power lamps were horizontally operated by a digital high-speed power supply with a control cycle below 10 ms (Teck-Tiger, Diez). It used a non-sinusoidal fixed frequency of 20 kHz with a maximum power output of 35 kW. Current and voltage values reached levels of roughly 30 A and 1 kV respectively. Some current waveforms for different power levels are shown in figure 3.

The used waveform is responsible for very low electromagnetic interference (EMI) and fulfills the CE (Conformité Européenne) specifications. Together with a total mass of about 80 kg, which lies greatly below typical weights of inductive ballasts, it was easy to integrate the compact device into the experimental setup, which was placed on an optical workbench shown in figure 4.

The lamps were cylindrically shaped, quartz made, without outer bulb and with tungsten coil electrodes. The inner diameter ranged between 18 and 25 mm and the length between 1.1 and 1.7 m. The filling consisted of some mg of Hg and about 20 mbar Ar as ignition gas. Figure 5 shows the end part of such a lamp.

Forced air cooling was necessary during lamp operation. In contrast to their industrial application a special lamp housing was used allowing good access for diagnostics and good thermal and UV radiation shielding. For diagnostic purposes special windows for the IR (ZnSe) and UV (sapphire) spectral range were inserted into the housing.

Spectroscopic plasma diagnostic measurements [16,17] were carried out in the middle of the lamp with a 0.5 m-spectrograph (SpectraPro 500i, Acton) combined with a 2048x2048 CCD-camera (Pixis, Roper). Two different gratings with 600 l/mm and 2400 l/mm were used in order to obtain overview spectra and high resolution lineshapes necessary for optical temperature measurements. The plasma temperature was determined with the Bartels method [18,19,20] from the radiance of the self-reversed maxima of the optically thick Hg spectral line at 546 nm. The electrodes and the lamp ending were observed in the visible and infrared spectral range. Especially the use of a fast 12 bit camera (SensiCam, PCO) with a minimal exposure time of 1 µs allowed phase resolved observations of the arc attachment mode on the electrode (figure 7) as well as the movement of the arc near the electrode region. The arc in the middle of the lamp was observed with a 8 bit camera (FlashCam, PCO) (figure 2).

Additionally, the electrode region was imaged using a simple optical system. At the image plane a photo diode (IPL-10531, IPL) with external low noise power supply was installed and the optical signal from different positions around the electrode was recorded and analyzed with a fast spectrum analyzer (FSEA20, Rohde&Schwarz) to detect plasma fluctuations. With this method introduced by Olsen [START_REF] Olsen | [END_REF] the resonance frequencies of the lamp for a given power could be recorded at a glance.

A photo diode measurement produces a signal caused by the integrated effect of the light flux originating from its line of sight. Spatial information about the fluctuation intensities which are coupled to the local temperature of the radial profile T (r) is lost.

So it was not possible to distinguish from which radial position the signal came from. Measurements using acoustic detection methods (eg. microphones) as suggested by Stromberg and Schaefer [6] were not possible due to strong vibrations and high noise level caused by the air cooling of the lamp.

Results and Discussion

Measurements on high power lamps for varying input power were performed on standard lamps used in industrial applications as well as on a limited number of specially prepared lamps with changed geometry and filling. All lamps operate stably at the nominal power of 35 kW, but a reduction of the input power by 50% (or less) leads in most cases to arc constriction and instability. In these cases arc distortion occurred combined with flickering and arc rotation. This was easily observable as increasing light flickering and arc deflection when the operating parameter approaches the unstable region.

The observed instabilities were identified to be caused by acoustic resonances because of its typical hysteresis which is not present i.e. in helical instabilities [21,22] exhibiting a similar appearance. Also the local magnetic field induced at the electrodes by currents of up to 30 A for maximum lamp power seem not to play a dominant role as the arc is stable at such high currents. The observed instabilities occurred at much smaller currents than its maximal value.

Plasma temperature profiles were measured in the middle of the lamp with temperature values ranging from 5500 K to 7000 K in the plasma core.

The temperature profile can be retraced in the high temperature central region of the arc where the self absorption of the optically thick line results in a self reversion. By assuming a parabolic like shape a corresponding profile can be fitted using:

T parabolic (r) = T (0) -(T (0) -T wall (R)) r R m ( 5
)
Here T (0) is the arc core temperature, R the inner radius of the discharge volume and T wall (R) the glass wall temperature. It should be mentioned that a true parabolic profile is obtained only for an exponent value m = 2 which is not necessarily the best fit parameter for a given profile. Usually values around 2 are found to give best profile fits. With the temperature profile a mean gas temperature can be computed using

< T mean >= T (r)rdr rdr . ( 6
)
Resonance frequencies were computed for every lamp as shown exemplarily for the reference lamp (1400 mm length) in figure 9. This figure also shows measured values for comparison which were taken from photo diode measurements presented in figure 6. The measurements were carried out with the photo diode near the electrode and in the lamp center. A Fourier analysis of its signal delivered a light 'fluctuation spectrum' which reveals the resonance frequencies at a glance as described by Olsen [START_REF] Olsen | [END_REF]. Furthermore, in figure 6 it is evident that the light fluctuation amplitudes are not homogeneously distributed along the lamp axis.

The measured spectrum reflects also the spectrum of the driving signal. The prominent peak at 40 kHz is situated at the double driving frequency which corresponds to the frequency of power input. The peak at 20 kHz corresponds to the frequency of mode switching at the electrode (compare with figure 7). In our case the change between cathode and anode at low powers also showed a change of the arc attachment mode which alternated between a spot-and diffuse mode. This can be seen in figure 7 for one electrode observed during different phases of a current period. Due to horizontal lamp operation heat convection along the lamp axis was obviated leading to expected identical thermal conditions at the electrodes. For this reason the second electrode was not observed and is expected to behave the same, similar to observations in high pressure mercury lamps [23]. The lamp observed in figure 6 was driven with 25 kW power. This assured always a stable operation. Any increase of input power increased the temperature, any decrease could lead to instability. In this sense the lamp was operated near the borderline between the regions of stability and instability. Once an instability was established, it was not possible to get rid of it by further decreasing the power due to a strong hysteresis as shown in figure 8. Only a jump to much higher input power eliminated the instability. A computation of resonance frequencies in dependence on temperature should agree with the measured peaks found in figure 6 and hence help to understand the lamp behavior with respect to stability. A comparison of measurement and computation is presented in figure 9. In this diagram the resonance frequencies are calculated using equations (1-4) for the reference lamp as a function of the temperature. Applying Gaussian error propagation to equations 1 and 2 the sensitivity between frequency and temperature in terms of relative errors can be derived giving ∆ rel f = C • 0.5 • ∆ rel T 0.5 , with C being a constant for longitudinal or transversal modes respectively. However, the more relevant source of uncertainty is the steep slope of the curves in the diagram.

The measured resonance frequencies in figure 9 are indicated by vertical lines. The slope of the calculated curves for different modes depends on the frequency -it is steep at low frequencies and becomes flat as the frequency increases. Due to the strong slope of the first resonance curve around 8 kHz the corresponding temperature interval is initially broad but it becomes relatively narrow if combined with the temperature intervals defined by the slope of the following resonance curves. Wall temperature measurements delivered values within an interval from 700 K to 1000 K leading to a mean temperature ranging from 3300 K to 3500 K.

Longitudinal resonances are located at much lower frequencies than transversal resonances of the same mode number. The first computed longitudinal mode starts at roughly 180 Hz at the corresponding mean temperature, whereas the transversal modes begin above several kHz. For this reason the longitudinal modes are not included in the diagram. This huge frequency spread between longitudinal and transversal modes is caused by the aspect ratio of this lamp type as the length and inner radius determine the frequencies in equation 1 and 2. Longitudinal modes seem to be dominant in the low frequency region which is consistent with the observation of lamps operated at 50 and 60 Hz, where the arc was stable in the first case and unstable in the second case as the slightly higher driving frequency or its harmonics were able to excite a resonance. On the other hand transversal resonances become important when the driving waveform contains a significant high frequency part [START_REF] Jongerius | [END_REF].

For every intersection point between a measured resonance and the computed resonance curve a hysteresis of ± 500 Hz was assumed. This value is justified by experience operating HID lamps where a typical hysteresis had a width of roughly 1 kHz. Usually the hysteresis was much higher near strong resonances. With the assumed hysteresis a temperature interval linked to the respective resonance was defined. For subsequent resonance modes along the frequency axis the width of this temperature interval was successively reduced due to the decreasing slope of the resonance curves. Consequently, a narrow temperature interval with a mean temperature of about 3450 K was found which corresponds to an assumed wall temperature of 920 K.

On the other hand, if one would like to determine the mean temperature and therefore the pressure in the lamp it is useful to take more than one frequency for a precise calculation.

Many of the measured resonances were found to be located in the marked mean temperature interval calculated from the measured temperature profile. Ovals indicate measured frequency peaks for which such a computed resonance can not be found in the corresponding temperature interval. Beside the possibility of electrical distortions coupled to the signal a more satisfactory explanation can be given if wave coupling is taken into account, as a coupling of two waves can transfer energy to other frequencies [24]. In this case the missing correlation between measurement and computation disappears if the driving frequency and its harmonics couple to the first transversal resonance mode. Within the relevant temperature interval its frequency lies roughly at 12 kHz. A coupling of the 20 kHz driving frequency leads to beat tones at the difference and sum frequencies with values of 8 kHz and 32 kHz respectively. The further coupling of the harmonics at 40 kHz and 60 kHz with the same mode leads to beat tones at 28, 52, 48 and 72 kHz covering all missing consistency.

The important finding for lamp operation is a computed resonance for T = 3450 K situated very near to the 20 kHz driving frequency. This is the main factor for lamp distortion and explains the observed behavior of flickering. The hysteresis in the lamp dysfunction as shown in figure 8 might be explained in such a way, that an increase of input power pushes the mean temperature upward increasing the distance to the resonance curve assuring stable operation. Reducing the input power lowers the mean temperature and the resonance approaches the driving frequency. For this particular lamp geometry this resonance curve is always crossed at the same input power (mean temperature) marking the transition between the regions of stable and unstable operation.

For the sake of completeness the formal possibility to correlate all measured peaks to resonance curves calculated for a much lower temperature T low of roughly 1270 K should be mentioned. For this temperature all observed peaks in the measured spectrum could be assigned, neglecting wave coupling. However, the required colder mantel around the hot central plasma, to allow the modes to exist, cannot be assumed as it would lead to an unrealistically low pressure for the particular lamp volume, mean temperature and fill.

Conclusion

The influence of length, inner radius and filling was investigated for high power lamps concerning their stable operation under power variation. All lamps behaved in a way corresponding to the picture of standing pressure waves. Best results were obtained with lamps of reduced inner radius. In this case the number of possible radial modes was effectively reduced and therefore the region of stable lamp operation considerably extended.

Peaks extracted from the measured light fluctuation spectra correspond well with transversal standing waves computed for the mean gas temperature. Considering wave coupling between the first transversal mode with the driving frequency and its harmonics beat tones occur which explain all additional measured peaks well.

Arc instability starts usually at the electrodes and is amplified when the resonance condition is reached, the excitation extends a certain level and lasts long enough to allow growth of the instability amplitude. In frequency space the resonance exhibit a Lorentzian profile [25,26] with side bands which in principle allow a coupling of the resonance to distant excitation frequencies. This means that a lamp should be operated in a region far away from its acoustic resonances.

Our calculations revealed a resonance free window for frequencies below 20 kHz. This encouraged us to develop a new electronic ballast which should be able to vary the driving frequency in order to omit the lamp specific acoustic resonance frequencies. Consequently such a ballast (DIMAC, Diez) was built. With this prototype the driving frequency can be varied between 10 and 20 kHz with a maximum output power of 13 kW. Much higher driving frequencies in the range of some 100 kHz at which acoustic resonances should not occur [START_REF] Wada | [END_REF][START_REF] Siessegger | Proc. of the LS-11[END_REF] were not considered as well as frequencies below the frequency of the lowest transversal mode due to special electrical properties of the ballast. Its high efficiency is optimized for the used frequency band.

First experiments with the new ballast are promising. Even though the development is not finished yet, it shows a good matching with the lamp and is able to avoid instabilities by changing the driving frequency.

Despite these good results it should be reminded that both the ballast and the lamp shape have to be adjusted in order to effectively avoid acoustic resonances as already stated in [START_REF] Wada | [END_REF][START_REF] Kettlitz | IEEE Industry Applications Society Annual Meeting[END_REF].

Finally, our approach of looking for 'resonance free windows' seems to be a handy tool to estimate actual lamp behavior with respect to stability. The correlation between computed resonance frequencies, the mean temperature and geometry can be checked easily and adjusted in such a way, that acoustic resonances are not excited. This was successfully tested in additional experiments on different commercial 150 W type HID

Figure 1 .

 1 Figure 1. Stable and unstable arc near the electrode of a high power lamp. The fluctuations are concentrated at the electrode region and in most cases instabilities began there.

Figure 2 .

 2 Figure 2. Distorted arc of a high power lamp during unstable operation caused by acoustic resonances.

Figure 3 .

 3 Figure 3. Current waveform of the power supply (Teck-Tiger) for different power levels at a fixed frequency of 20 kHz.

Figure 4 .

 4 Figure 4. Experimental setup for high power lamp diagnostics. The electronic ballast was integrated into the setup.

Figure 5 .

 5 Figure 5. One of the analysed lamps -quartz glass tube ending with electrode, electrical connection and infrared reflection coating.

Figure 6 .

 6 Figure 6. Measured light fluctuation spectrum at two different positions in the reference lamp driven with 25 kW. The strongest peaks correspond to the frequency of power input and its harmonics. The peak at 20 kHz lies close to a computed resonance frequency in the corresponding temperature region.

Figure 7 .

 7 Figure 7. Phase resolved measurements of the arc attachment at the electrode of a high power lamp. The anode phase is shown on the left, the cathode phase on the right-hand side. The working frequency of 20 kHz corresponds to a period of 50 µs. The images were taken with an exposure time of 10 µs at 15 kW input power.

Figure 8 .

 8 Figure 8. Hysteresis of the 35 kW reference lamp.

Figure 9 .

 9 Figure 9. Calculated resonance frequencies compared with measured ones for the reference lamp at 25 kW. The resonance frequencies taken from figure 6 are indicated by the dotted vertical lines. Two mean temperatures result from two wall temperatures of 700 K and 1000 K respectively. A hysteresis of ± 500 Hz is assumed for the intersection points of calculated and measured frequencies at the mean temperature.

Acknowledgments

We would like to thank the neoplas GmbH for all support and the good collaboration. This work was promoted within the PRO INNO II project of the Federal Ministry of Economics and Technology.

lamps, where acoustic resonances were intentionally established or avoided.