Benedikt Obermayer
email: obermayer@physik.lmu.de

Erwin Frey
email: frey@physik.lmu.de

Error thresholds for self-and cross-specific enzymatic replication

Keywords: origin of life, quasispecies theory, enzymatic self-replication, higher-order catalysis

The information content of a non-enzymatic self-replicator is limited by Eigen's error threshold. Presumably, enzymatic replication can maintain higher complexity, but in a competitive environment such a replicator is faced with two problems related to its twofold role as enzyme and substrate: as enzyme, it should replicate itself rather than wastefully copy non-functional substrates, and as substrate it should preferably be replicated by superior enzymes instead of less-efficient mutants. Because specific recognition can enforce these propensities, we thoroughly analyze an idealized quasispecies model for enzymatic replication, with replication rates that are either a decreasing (self-specific) or increasing (cross-specific) function of the Hamming distance between the recognition or "tag" sequences of enzyme and substrate. We find that very weak self-specificity suffices to localize a population about a master sequence and thus to preserve its information, while simultaneous localization about complementary sequences in the cross-specific case is more challenging. A surprising result is that stronger specificity constraints allow longer recognition sequences, because the populations are better localized. Extrapolating from experimental data, we obtain rough quantitative estimates for the maximal length of the recognition or tag sequence that can be used to reliably discriminate appropriate and infeasible enzymes and substrates, respectively.

Introduction

The acclaimed experimental finding [START_REF] Cech | Nobel lecture. Self-splicing and enzymatic activity of an intervening sequence RNA from tetrahymena[END_REF] that RNA not only stores genetic information but also provides catalytic function has inspired the RNA world scenario [START_REF] Gilbert | Origin of life -the RNA world[END_REF], a hypothesis for the starting point of Darwinian evolution at the origin of life through self-replication of RNA-like polynucleotides. Substantial progress over the last decades demonstrates the capability of RNA enzymes (ribozymes) to catalyze diverse chemical reactions [START_REF] Doudna | The chemical repertoire of natural ribozymes[END_REF][START_REF] Lilley | Structure, folding and mechanisms of ribozymes[END_REF][START_REF] Joyce | Forty years of in vitro evolution[END_REF], among them the polymerization of as many as 20 nucleotides to a template molecule [START_REF] Johnston | RNA-catalyzed RNA polymerization: accurate and general RNAtemplated primer extension[END_REF][START_REF] Zaher | Selection of an improved RNA polymerase ribozyme with superior extension and fidelity[END_REF], and even replication through template-directed ligation reactions involving short RNA oligomers as building blocks [START_REF] Lincoln | Self-sustained replication of an RNA enzyme[END_REF]. However, the currently known systems are not yet capable of Darwinian evolution, lacking either the ability to replicate molecules as long and complex as themselves or to introduce heritable variation.

For theorists, the focus has mainly been on whether the information content of a self-replicating molecule can be maintained in the presence of replication errors, usually employing Eigen's well-known quasispecies theory [START_REF] Eigen | The molecular quasi-species[END_REF] for the self-replication of L-nucleotide sequences. In this model, replication errors occur with an error probability μ per single nucleotide, and the replication rates are taken as functions only of the template sequence. This mapping from genotype to replication rate (or fitness) presents a considerable challenge for evolutionary models, and even generic structural features of such fitness landscapes, such as the extent of ruggedness and epistasis, are under debate. Valuable insight comes from experiments on mutagenized ribozymes [START_REF] Kun | Real ribozymes suggest a relaxed error threshold[END_REF] and from computer simulations using RNA secondary structure as phenotype [START_REF] Huynen | Smoothness within ruggedness: The role of neutrality in adaptation[END_REF][START_REF] Takeuchi | Phenotypic error threshold; additivity and epistasis in RNA evolution[END_REF], which indicate a significant degree of neutrality around an optimal sequence, reducing the deleterious effects of mutations. Other theoretical studies used a large number of different idealized fitness landscapes, e.g., with a single peak at one fittest "master" sequence [START_REF] Swetina | Self-replication with errors -a model for polynucleotide replication[END_REF][START_REF] Leuthäusser | An exact correspondence between Eigen's evolution model and a two-dimensional Ising system[END_REF][START_REF] Schuster | Stationary mutant distributions and evolutionary optimization[END_REF][START_REF] Woodcock | Population evolution on a multiplicative single-peak fitness landscape[END_REF][START_REF] Galluccio | Exact solution of the quasispecies model in a sharply peaked fitness landscape[END_REF][START_REF] Hermisson | Mutation-selection balance: Ancestry, load, and maximum principle[END_REF][START_REF] Peliti | Quasispecies evolution in general mean-field landscapes[END_REF][START_REF] Saakian | Exact solution of the Eigen model with general fitness functions and degradation rates[END_REF][START_REF] Saakian | Phase diagram for the Eigen quasispecies theory with a truncated fitness landscape[END_REF], with some rather generic results [START_REF] Wiehe | Model dependency of error thresholds: the role of fitness functions and contrasts between the finite and infinite sites models[END_REF][START_REF] Jain | Adaptation in simple and complex fitness landscapes[END_REF]: the population in sequence space is characterized by a broad mutant distribution (a quasispecies) localized about the master sequence for mutation probabilities smaller than a critical value μ c (the error threshold), while it consists of random sequences (it is delocalized) for larger values. Because the error threshold μ c ∼ 1/L is usually inversely proportional to sequence length, the problem arises whether the maximally sustainable complexity of a self-replicator suffices to perform the complex task of self-replication [START_REF] Eigen | The hypercycle -a principle of natural selforganization. Part B: The abstract hypercycle[END_REF].

In a prebiotic context, it is important to emphasize that using a fitness landscape where rates depend only on the template sequence pertains to nonenzymatic rather than enzymatic replication, because in the latter case the replication rates also depend on the concentrations and the characteristics of involved enzymes. For RNA, the potential for non-enzymatic replication is questionable, given that template-directed polymerization or ligation seems limited to short molecules with rather specific sequences [START_REF] Von Kiedrowski | A self-replicating hexadeoxynucleotide[END_REF][START_REF] Acevedo | Non-enzymatic transcription of an oligodeoxynucleotide 14 residues long[END_REF][START_REF] Wu | Nonenzymatic template-directed synthesis on oligodeoxycytidylate sequences in hairpin oligonucleotides[END_REF][START_REF] Orgel | Prebiotic chemistry and the origin of the RNA world[END_REF]. Moreover, experimental and simulation studies demonstrate a strongly disadvantageous tendency for elongating side-reactions at the cost of replication [START_REF] Fernando | A stochastic model of nonenzymatic nucleic acid replication: "elongators" sequester replicators[END_REF]. Hence, although there have been speculations [START_REF] Pace | RNA catalysis and the origin of life[END_REF], it remains unclear how a single more complex RNA should literally copy itself [START_REF] Joyce | Forty years of in vitro evolution[END_REF][START_REF] Szostak | Synthesizing life[END_REF].

Enzymatic replication is more plausible [START_REF] Orgel | Molecular replication[END_REF], but raises the question whether high replication efficiency (high fitness) is a property of the substrate or the enzyme. In the latter case, a superior replicase does not enjoy a selective advantage, because it replicates non-functional mutant templates just as well as itself, while it is not guaranteed in the former case that a superior template is functional at all. Likewise, mutations generate substrates that are replicated less efficiently, but they also produce less-efficient enzymes [START_REF] Smith | Hypercycles and the origin of life[END_REF], thus affecting the replication rates of all potential substrates. On theoretical grounds, one should expect that a superior replicator is a good enzyme and a good substrate at the same time. As enzyme, it should therefore replicate only functional substrates, and as substrate, it should be replicated preferably by efficient enzymes. It has long been realized that in a competitive environment these propensities are crucial for the emergence, improvement and perpetuation of replicase activity, which is essentially an altruistic trait that is not by itself selected for [START_REF] Michod | Population biology of the first replicators -on the origin of the genotype, phenotype and organism[END_REF]. The commonly proposed solution to render enzymatic replication evolutionarily stable is to impose a form of group selection, e.g., via compartmentalization in vesicles, in order to keep similar molecules closely together [START_REF] Szathmáry | Group selection of early replicators and the origin of life[END_REF][START_REF] Alves | Group selection models in prebiotic evolution[END_REF][START_REF] Fontanari | Coexistence and error propagation in pre-biotic vesicle models: A group selection approach[END_REF], even though this requires a simultaneous and coordinated emergence of replicators and protocells [START_REF] Szostak | Synthesizing life[END_REF]. However, known ribozymes act with moderately or even strongly substrate-specific efficiency [START_REF] Joyce | Forty years of in vitro evolution[END_REF], such that specific recognition may also have a significant influence. Since unspecific reactions require sophisticated substrate-binding properties that could well have been a rather late invention in prebiotic evolution [START_REF] Johnston | RNA-catalyzed RNA polymerization: accurate and general RNAtemplated primer extension[END_REF], it seems natural to assume that replication efficiency should depend strongly on the interaction between enzyme and substrate.

In this paper, we analyze a general model of enzymatic replication accounting both for varying degrees of specificity and the characteristically broad quasispecies distributions in order to address the consequences for the error threshold. Similar to models for the evolution of regulatory DNA motifs [START_REF] Gerland | On the selection and evolution of regulatory DNA motifs[END_REF], we assume that specificity depends on the quality of binding to some recognition or tag sites [START_REF] Weiner | Transfer RNA-like structures tag the 3' ends of genomic RNA molecules for replication -implications for the origin of protein-synthesis[END_REF]. Idealizing this condition, we use replication rates that depend on the Hamming distance between these sequence regions of enzyme and substrate via a decreasing (self-specific) or increasing (cross-specific) function. After formulating the model and discussing our analytical approach in a methods section, we show for these two scenarios results from stochastic simulations and numerical solutions of deterministic rate equations. Using our analytical toolbox, we discuss the resulting localization conditions, error thresholds and the phase diagram. In our conclusions, we use experimental values for polymerization rates to obtain simple estimates for the maximum number of nucleotides that can be used for recognition.

Model

In the framework of quasispecies theory, each molecule is characterized by its sequence S i = (σ

(i) 1 . . . σ (i) L ) of L binary nucleotides σ (i) ∈ {0, 1}
. In an infinitely large population, its concentration X i evolves according to the deterministic rate equations [START_REF] Eigen | The molecular quasi-species[END_REF])

Ẋi = j M ij R j X j -X i j R j X j . (1)
Here,

M ij = μ dij (1 -μ) L-dij
is the mutation probability between sequences S i and S j with Hamming distance

d ij = |σ (i) -σ (j)
|, where μ is the error probability per single nucleotide (usually called "mutation rate"), and the replication rate R i of sequence S i is given by:

R i = A i + j B ij X j . (2)
Whereas the non-enzymatic rate A i depends only on the genotype S i , the second term implies frequency-dependent selection and makes our model intrinsically nonlinear. It encodes the catalytic interactions of two molecules: B ij measures how well S j catalyzes the replication of S i . The second term in Eq. (1) ensures the normalization j X j = 1, and a degradation term -D i X i therefore drops out of Eq. ( 1) since we assume that the decay rate D i ≡ D is sequenceindependent for simplicity. Note that actual RNA sequences replicate via a complementary intermediate, while our idealized model assumes direct replication in a single step. It has been shown that these two approaches are essentially equivalent for the symmetric situations considered here [START_REF] Stadler | Complementary replication[END_REF].

To capture the pertinent features of a situation where replicase enzymes prefer to replicate themselves instead of their competitors, we assume that the quality of specific recognition influences catalytic rates more strongly than the actual genotypes of enzyme and substrate. Hence, we effectively only model the recognition regions of ribozymes, which are often clearly separated from the catalytic domains [START_REF] Lilley | Structure, folding and mechanisms of ribozymes[END_REF]. Of course, the proper function of the latter region is indispensable, but in our idealized model we neglect the influence of mutations: as we have shown previously in a simple model, their effect on the error threshold is largely independent from the more interesting consequences of mutations in the recognition region [START_REF] Obermayer | Escalation of error catastrophe for enzymatic self-replicators[END_REF]. The sequence length L is thus restricted to the number of nucleotides that take part in recognition. Probably mediated via specific base-pairing interactions [START_REF] Doudna | The chemical repertoire of natural ribozymes[END_REF], the quality of recognition can be taken as function of the number of mismatches between the binding sites of enzyme and substrate, and we let the catalytic matrix B ij therefore depend via a specificity function f (d) only on the Hamming distance d ij between enzyme and substrate. Further, because rate enhancements through ribozyme catalysis can be substantial [START_REF] Doudna | The chemical repertoire of natural ribozymes[END_REF], such that non-enzymatic replication rates are comparably small (if nonzero at all), we neglect their genotype dependence altogether and choose a flat fitness landscape for A i :

A i ≡ α, B ij = βf (d ij ). (3)
Because we are interested in the stationary state, the parameter α (if nonzero) merely sets the time scale while β measures the selection strength.

We will first analyze a scenario for self-specific replication, where replication rates increase with similarity of enzyme and substrate. To have the degree of specificity explicitly tunable via a parameter p, we use the specificity function

f s (d) = (1 -d/L) p . ( 4
)
As a contrasting example, we will then analyze a similarly defined specificity function where replication rates increase with complementarity between enzyme and substrate:

f c (d) = (d/L) p . (5)
The numerical and analytical approach presented in the next section can readily be applied to other functional forms of the specificity function.

Methods

For reasonably large sequence length, the full 2 L -dimensional system Eq. ( 1) can only be analyzed using stochastic simulations in a finite population of N sequences. Here, we employ the straightforward stochastic simulation algorithm used by [START_REF] Wilke | Dynamic fitness landscapes in molecular evolution[END_REF]. At time t each sequence S k , present in N k copies, has a probability p 0,k = N k / i N i (1 + R i ) to be copied without mutations into the population at time t + 1, and a probability

p mut,jk = M jk R k N k / i N i (1 +
R i ) to be selected and mutated into sequence S j . Following initialization, our observables of interest are measured by averaging over time after reaching a stationary state.

In order to derive analytical results, we exploit that the stationary states of Eq. ( 1) are localized about a particular "master" sequence S * , which is necessary to preserve its information content. In contrast, delocalization indicates that such a sequence S * cannot be maintained due to replication errors. Since the possibility of non-trivial dynamics such as periodic orbits cannot be excluded for general replicator-mutator equations like Eq. ( 1) [START_REF] Stadler | Dynamics of small autocatalytic reaction networks. 2. Replication, mutation and catalysis[END_REF], there might also be other reasons for the absence of localization. However, we did not find any signs of complex dynamical behavior in our simulations, and one can easily convince oneself that the delocalized state, where all sequences have the same concentration and therefore identical replication rates, can lead to localization: since replication rates are essentially proportional to concentration, stochastic concentration fluctuations imply higher rates. Unless hindered by excessive mutations, this can induce a transition to a state localized about some randomly chosen master sequence, which then also has the highest replication rate (see [START_REF] Obermayer | Escalation of error catastrophe for enzymatic self-replicators[END_REF] for a visualization). Such "fixation" events are very similar to the phenomenon of consensus formation, e.g., in language dynamics [START_REF] Blythe | Generic modes of consensus formation in stochastic language dynamics[END_REF]. Note that our idealized replication rates do not predetermine any specific master sequence for localization. This symmetry would be broken in a full model where replication rates depend on the full genotypes of enzyme and substrate (and not just the Hamming distance between their recognition regions).

Given localization about S * , we can significantly reduce the dimensionality of Eq. ( 1) by lumping all sequences S i with a Hamming distance k to S * together into "error class" k. Without loss of generality, we assume that S * = (00 . . . 0).

This well-known procedure [START_REF] Schuster | Stationary mutant distributions and evolutionary optimization[END_REF][START_REF] Woodcock | Population evolution on a multiplicative single-peak fitness landscape[END_REF] permits to formulate reduced rate equations formally equivalent to Eq. ( 1) in terms of new variables x k denoting the concentration of error class k in the population:

ẋk = ji m kj [a j δ ji + b ji x i ]x j -x k ji [a j δ ji + b ji x i ]x j . (6)
Here and in the following, we use lowercase letters for all variables in the reduced system. For our model Eq. ( 3), the non-enzymatic rates are given by a i ≡ α. The accordingly reduced mutation matrix and the catalytic matrix depend only on the Hamming distance between pairs of sequences in different error classes (measured with respect to the master sequence), which allows us to combinatorially assess all possibilities for their relative distance. The total probability of distributing 0 ≤ kj + 2 ≤ L mutations to move a sequence from error class j into error class k has been derived previously as [START_REF] Woodcock | Population evolution on a multiplicative single-peak fitness landscape[END_REF])

m kj = L -j k -j + j (1 -μ) L-(k-j+2 ) μ k-j+2 . ( 7
)
The replication rate Eq. ( 2) also depends on the frequency of each sequence in each error class. With the homogeneity assumption that all L j sequences in class j are equally populated, this complication can be resolved, and the reduced matrix b ij reads analogously

b ij = β n L -i j -i + n i n f (j -i + 2n) L j . ( 8
)
Numerical solutions to the (L + 1)-dimensional rate equations given in Eq. ( 6) can now easily be found by means of standard algorithms.

A more detailed understanding of these solutions can be obtained from the population distribution's normalized first moment a = k /L = σ , which as the mean Hamming distance k = k kx k to the master S * = (00 . . . 0) characterizes the width of the distribution and measures the mean value σ of each sequence's binary nucleotides. Writing down an equation for the first moment of Eq. ( 6) requires a hierarchy of expressions for higher moments, which can be truncated by means of a moment closure technique. For the self-specific case, we assume that the stationary Hamming distance distribution is approx-

imately binomial, x k ≈ L k a k (1 -a) L-k
, because this reproduces the expected distribution in the limits a → 0 (complete localization about one sequence) and a → 1/2 (the delocalized state, where the binary nucleotides are random numbers). Moreover, it solves the rate equations Eq. (1) exactly for linear fitness landscapes without epistasis [START_REF] Woodcock | Population evolution on a multiplicative single-peak fitness landscape[END_REF] and for an extension of the quasispecies model to a game theory setting [START_REF] Lässig | Evolutionary games and quasispecies[END_REF]. With this binomial ansatz, a is the population distribution's only parameter, and for our model of the replication rates Eq. (3), it obeys the equation

(1 -2a) μL [α + βS(1 -2a(1 -a))] -a(1 -a)(1 -2μ)βS (1 -2a(1 -a)) = 0. (9)
This equation for a, which is one of our main analytical results (see Appendix A.1 for a derivation), holds for any specificity function f (d), which enters via the function

S(x) = k L k f (k) x L-k (1 -x) k . ( 10
)
An intuitive interpretation of this auxiliary function derives from recognizing that the quantity α + βS(1 -2a(1a)) measures the mean replication rate (or mean fitness) j R j X j of the population if a solves Eq. ( 9). Such solutions a(μ) can be obtained explicitly only in special cases where S(x) attains a simple form (see Appendix A.2), but we can easily solve for μ(a) and invert graphically to obtain a bifurcation diagram:

μ(a) = 2 + α + βS(1 -2a(1 -a)) βa(1 -a)S (1 -2a(1 -a)) -1 . ( 11
)
For cross-specific replication with the specificity function Eq. ( 5), we expect two equivalent subpopulations localized about complementary sequences, which corresponds to a superposition of binomial distributions:

x k ≈ 1 2 L k a k (1 - a) L-k + (1 -a) k a L-k .
To obtain an equation similar to Eq. ( 9) for their mean widths a, we cannot use the first moment of the reduced rate equations (it vanishes by construction, because the distribution is symmetric about

a = 1/2), but use the second moment Δk 2 = k (k -k ) 2 x k ,

Results and Discussion

Self-specific replication

Our first scenario is concerned with self-specific replication, where replication rates increase with similarity of enzyme and substrate through the specificity function Eq. ( 4). As argued in the preceding section, stochastic fluctuations increasing the frequency of one particular sequence also increase its replication rate, such that the population can localize about this master sequence S * . Fig-

ure 1 shows the stationary Hamming distance distribution x k , which measures the concentration of sequences with k mutations relative to S * , for different degrees of specificity from the linear case p = 1 to complete self-specificity p → ∞, where enzyme and substrate have to be identical. We compare results from stochastic simulations to numerical solutions of the reduced rate equations, Eq. ( 6). The excellent agreement between simulation and deterministic theory justifies the homogeneity assumption made in symmetrizing the specificity matrix (see Eq. ( 8)). In the simulations shown in Fig. 1, we initialized all sequences at a predetermined master sequence in order to avoid noise from the intrinsically stochastic "fixation" events, but we also tested other initial conditions with modest inhomogeneities, which were quickly "washed out" and did not give rise to measurable differences in the stationary state.

The limit p → ∞ with f s (d) → δ d,0 , depicted in Fig. 1(c), leads to a generalized Schlögl model of auto-catalytic replication, which has been partly analyzed by [START_REF] Stadler | Dynamics of small autocatalytic reaction networks. 2. Replication, mutation and catalysis[END_REF]. In this limit we can employ the well-known "error-tail" approximation [START_REF] Schuster | Stationary mutant distributions and evolutionary optimization[END_REF]: we define x 0 as the concentration of the master sequence, α + βx 0 its replication rate and (1μ) L the probability not to have a mutation. All other sequences are lumped together in the error tail with concentration 1x 0 and replication rate α (the concentration of suitable replicase enzymes is so small that the frequency-dependent term in the replication rate does not contribute). Neglecting back mutations from the error tail into x 0 (corresponding to the large-genome limit), we obtain the simple

equation ẋ0 = (α + βx 0 )(1 -μ) L x 0 -x 0 r, ( 12
)
with r = (α + βx 0 )x 0 + α(1x 0 ) the mean replication rate. In the stationary state, we easily find that the delocalized state x 0 = 0 is stable for all μ, while a branch of solutions with nonzero x 0 emerges for (1μ) L > 2( α(α + β)α)/β through a discontinuous transition (see also [START_REF] Campos | Error propagation in the hypercycle[END_REF]; Obermayer and Frey (2009); [START_REF] Wagner | Second order catalytic quasispecies yields discontinuous mean fitness at error threshold[END_REF] for similar results in related models).

Whereas for the somewhat related sharply-peaked fitness landscape, where only the master sequence has a higher replication rate, the error threshold arises through a continuous bifurcation (Baake and Wiehe, 1997), the discontinuity observed here expresses the qualitatively different behavior we previously termed "escalation of error catastrophe" [START_REF] Obermayer | Escalation of error catastrophe for enzymatic self-replicators[END_REF]: as the mutation rate grows, the proportion of fittest sequences, i.e., of the necessary replicase enzymes, is diminished and therefore their replication rate. This in turn reduces their concentration, until at the error threshold the concentration of enzymes

x 0 is not large enough to have them replicate with an efficiency sufficient for localization.

Although Eq. ( 12) approximates the exact result for x 0 very well (see the dashed line in Fig. 1(c)), it is valid only for L → ∞ (because back mutations are neglected) and p → ∞ (because the replication rate of the error-tail is taken as concentration-independent). In order to gain a more general perspective, we use the analytical solutions of Eq. ( 9) for the mean Hamming distance to the master (the population distribution's first moment a = k /L) obtained from the bifurcation diagram Eq. ( 11). A comparison between the exact solution a(μ)

resulting from the reduced rate equations Eq. ( 6) via numerical continuation 1 and Eq. ( 11) is shown in Fig. 2 for L = 8, α = β = 1 and different values of p in the specificity function Eq. ( 4). Recalling the symmetry of the original model, Eq. ( 3), namely that the population can localize about any sequence, the remaining reflection symmetry a → 1a about the delocalized solution a = 1/2 indicates that after symmetrization localization is only possible about the master sequence or its complement. The stable branches associated with localized solutions start for μ = 0 at a = 0 (or a = 1), i.e., full localization about one sequence (or its complement), and higher mutation rates give rise to broader distributions with larger mean, until these localized regimes disappear at critical mutation rates μ c denoted by circles. The delocalized solutions, on the other hand, gain stability at finite mutation rates μc (denoted by crosses), and we find bistability for μc < μ < μ c . A similar situation is encountered in formally related models of grammar evolution [START_REF] Nowak | Evolution of universal grammar[END_REF][START_REF] Komarova | Replicator-mutator equation, universality property and population dynamics of learning[END_REF].

Explicit expressions for the error threshold μ c = max μ(a) are available from Eq. ( 11) when the auxiliary function S(x) (defined in Eq. ( 10)) has a simple form (see Appendix A.2). For instance, μ c = β/(4αL + 2β(L + 1)) for p = 1, where specificity decreases linearly with the distance between enzyme and substrate.

1 AUTO software package available via http://indy.cs.concordia.ca/auto/.

Although an error threshold is absent for linear fitness landscapes without epistasis [START_REF] Woodcock | Population evolution on a multiplicative single-peak fitness landscape[END_REF], here the intrinsically nonlinear model gives a sharp transition even in this apparently similar case. For very strong specificity p → ∞ we find μ c = W[β/(eα)]/(2L) using Lambert's W-function, and we recover the result μ c = ln(β/α)/(2L) that can also be obtained from the error-tail approximation in the limit β α. Further, the delocalized state a = 1/2 is the only solution of Eq. ( 9) in the complete absence of specificity (f s (d) ≡ 1), supporting the intuition that unspecific replication does not suffice to preferentially maintain the information of one particular sequence. Interestingly, taking p → 0 in Eq. ( 4) gives the finite even though exponentially small value μ c = β/(2 L+1 (α + β)). This result implies that limited localization is possible even for very weak specificity (if p = 0 in the specificity function Eq. ( 4), enzymes replicate everything except their exact complement, because always

f s (L) = 0).
From the bifurcation diagram obtained via Eq. ( 11), we easily read off exact results for the value μc = μ(1/2) where the delocalized state gains stability. E.g.,

for the generalized Schlögl model p → ∞, we get μc = β/(α2 L+1 + 4β) [START_REF] Stadler | Dynamics of small autocatalytic reaction networks. 2. Replication, mutation and catalysis[END_REF]. This exponentially small yet finite value is consistent with our previous conclusion that the delocalized regime is stable for all values of μ within the error-tail approximation, Eq. ( 12), which holds for L → ∞. Further, we find that the two critical values μc and μ c are identical for p = 0, 1, 2.

Recognizing that Fig. 2 describes a pitchfork bifurcation at a = 1/2 and μ = μc , we infer that the two critical mutation rates are equal (μ c = μc ) whenever the pitchfork is supercritical, whereas bistability between localized and delocalized states for intermediate mutation rates μc < μ < μ c is possible in the subcritical case, leading to the discontinuous transition observed in Fig. 1. The bistability regime vanishes as the curvature μ (1/2) in the bifurcation diagram changes sign, which gives from Eq. ( 11) an approximate expression for the corresponding critical value of β:

β * = α S 2 (1/2) S (1/2) -2S (1/2) -S(1/2) -1 . ( 13
)
This general result can readily be evaluated for any specificity function entering the auxiliary function S(x) defined in Eq. ( 10). It predicts bistability for all Nevertheless, we emphasize that the result μ c = ln(β/α)/(2L) in the limiting case p → ∞ agrees with the value obtained from the error-tail approximation for large β. The remarkable performance of the binomial closure approximation can be appreciated in more detail from the projected phase diagrams shown in Fig. 4.

values β < β * = α 2 L-1 (L -2) if p → ∞ in
A noticeable feature of these phase diagrams is that the error threshold μ c increases for stronger specificity p (see Fig. 4(b)), which implies that higher mutation rates can be tolerated, i.e., that longer sequences can be maintained.

This seems at first counter-intuitive, because weaker specificity constraints on the recognition sequence should allow more mutational "freedom". However, as shown in Fig. 2 and Fig. 1, smaller values for p lead to much broader distributions: mutants are still reasonably well replicated by master enzymes, but the master is only moderately well (but not quite as efficiently) replicated by the mutants. The resulting broadening of the distribution effectively reduces the replication rate of the master and escalates the error catastrophe [START_REF] Obermayer | Escalation of error catastrophe for enzymatic self-replicators[END_REF]. Thus, the necessity for an enzymatic replicator to discriminate not only between functional and non-functional substrates, but also between efficient and unproductive enzymes, is again emphasized.

We finally want to remark on the case α = 0, i.e., no background level for the non-enzymatic replication rate. Most of the above results obtained from the binomial closure approximation can be simply evaluated for α = 0 (note that then β sets the timescale and drops out), but the strong specificity limit p → ∞ deserves extra attention. The error-tail approximation indicates that the error threshold vanishes (μ c → 1), but from Eq. ( 9) we find the exact result

μ c = μc = 1/6 for L 1, i.e., a macroscopic yet finite value (see Appendix A.2).
This remarkable result can be explained by recalling that the traditional result [START_REF] Eigen | The molecular quasi-species[END_REF][START_REF] Wiehe | Model dependency of error thresholds: the role of fitness functions and contrasts between the finite and infinite sites models[END_REF] for the error threshold depends on the replication advantage r of the master relative to a possibly small but finite value for the mutants. In our case, rates are directly proportional to concentration, and because the master sequence has a concentration of order 1, while in an infinitely large population distant mutants have concentrations of order 2 -L , this relative advantage itself is of order 2 L , and cancels the length dependence of the error threshold. In the corresponding non-enzymatic case, results for so-called "truncation" fitness landscapes have lead to some debate about the applicability of the error threshold concept in the presence of lethal mutations [START_REF] Wilke | Quasispecies theory in the context of population genetics[END_REF][START_REF] Summers | Examining the theory of error catastrophe[END_REF][START_REF] Takeuchi | Error-threshold exists in fitness landscapes with lethal mutants[END_REF][START_REF] Saakian | Phase diagram for the Eigen quasispecies theory with a truncated fitness landscape[END_REF]. Accordingly, we should cautiously note that our results for α = 0 will probably be affected when accounting for the effects of finite populations and the full dependence of the replication rates on the genotypes of enzyme and substrate.

μ c ≈ ln r/L (

Cross-specific replication

To increase the information content of replicating systems beyond the limited complexity of a single replicator, auto-catalytic reaction networks such as hypercycles [START_REF] Eigen | The hypercycle -a principle of natural selforganization. Part B: The abstract hypercycle[END_REF][START_REF] Stadler | Dynamics of small autocatalytic reaction networks. 2. Replication, mutation and catalysis[END_REF] have been proposed, where different molecular species catalyze each other's replication in a possibly complex interaction graph. Only very little is known for these systems regarding the issue of reaction specificity and the cross-interactions of each species' mutant clouds. The simplest conceivable networks are 2-member cross-catalytic hypercycles. While the two members (which actually replicate via complementary intermediates) need not be strictly complementary, we assume that they are sufficiently distinct that an idealized model requiring complementary recognition regions captures the essentials. This suggests to analyze the specificity function Eq. ( 5) where replication rates increase with Hamming distance between enzyme and substrate. In this case, we expect the formation of two sub-populations localized about complementary sequences, each catalyzing the replication of the other. The main question to be answered is how specificity affects coexistence.

In the standard Eigen model, quasispecies coexistence is prevented by competitive exclusion except in degenerate cases, because the "fittest" individuals take over the population [START_REF] Swetina | Self-replication with errors -a model for polynucleotide replication[END_REF]. Interactions between sub-populations, e.g., based on complementarity of binary traits (de Oliveira and [START_REF] De Oliveira | Complementarity and diversity in a soluble model ecosystem[END_REF], are known to enable coexistence. In our case, we expect novel conditions for coexistence: even though each subpopulation depends on the presence of the other for efficient replication, it is unclear how the possibly broad mutant distributions influence each other.

From numerical solutions to the reduced rate equations and simulation results, where we initialized the population split between the master sequence and its complement (see Fig. 5), we find that localization is only possible for p > 1.

This important result is confirmed through the binomial moment closure approx- In particular, the negligible interaction between the respective mutant clouds allows one to employ the error-tail approximation assuming independent species and error tails as in [START_REF] Campos | Error propagation in the hypercycle[END_REF][START_REF] Obermayer | Escalation of error catastrophe for enzymatic self-replicators[END_REF]. Further, we find that the pitchfork bifurcation is always subcritical, i.e., that μ c > μc .

We summarize our main results by plotting projected phase diagrams of μ c and μc as functions of p-1 and β in Fig. 7. This confirms that both critical mutation rates μ c and μc vanish linearly with p -1, because μ(a) ∝ p -1 as p → 1, which gives the sharp bound p > 1 for coexistence of two populations. Finally, the case α = 0 of zero non-enzymatic replication rate is similar to the self-specific case: the slight chance of distant mutants to find an appropriate enzyme gives an enormous replication advantage to the mainly populated master sequences and therefore macroscopic values for the error threshold (see Appendix B.3 for details).

Conclusion

Because enzymatic replication provokes the necessity for enzymes to favor functional substrates and for substrates to prefer efficient enzymes, we analyzed the effects of specific replication on the error threshold. In our model, which accounts for the full quasispecies mutant distributions, replication rates depend on the Hamming distance between enzyme and substrate via an adjustable degree of specificity. Combining stochastic simulations, numerical solutions of reduced rate equations and analytical solutions to a binomial closure approximation, we could analyze the entire phase diagram and assess how mutation rate μ, selection strength β and specificity degree p influence the localization about a master sequence in order to preserve its information content. We found that for self-specific replication very weak specificity suffices for localization, whereas stronger specificity gives more tolerance against mutations but leads to bistability with the delocalized regime of random sequences. In particular, the binomial closure approximation permits to obtain analytical expressions for the bifurcation diagram and an upper bound β * for the bistability regime, which can be evaluated for any specificity function f s (d). Apart from our special choice, Eq. ( 4), a mesa-shaped function would also be conceivable, in correspondence to fitness landscapes for transcription factor binding allowing for some "fuzziness"

or neutrality in the binding sequence [START_REF] Gerland | On the selection and evolution of regulatory DNA motifs[END_REF]. Preliminary results indicate that in this case the binomial closure approximation gives at least qualitative agreement as well. While our approximation is not restricted to large L, this limit can probably be more systematically be described using the maximum principle employed previously in quasispecies theory [START_REF] Hermisson | Mutation-selection balance: Ancestry, load, and maximum principle[END_REF][START_REF] Saakian | Exact solution of the Eigen model with general fitness functions and degradation rates[END_REF]. In the case of cross-specificity, we found that coexistence of subpopulations localized about complementary sequences is possible only if replication rates increase faster than linearly with Hamming distance.

Although our model is based on idealizing assumptions, we can extrapolate from currently available experimental data to obtain rough quantitative esti-mates for the maximal length L c of the recognition or tag sequences that can be used by replicase enzymes to specifically and reliably discriminate appropriate and useless templates (and vice versa). Considering that non-enzymatic template-directed polymerization rates are on the order of several hours to days per base [START_REF] Acevedo | Non-enzymatic transcription of an oligodeoxynucleotide 14 residues long[END_REF][START_REF] Wu | Nonenzymatic template-directed synthesis on oligodeoxycytidylate sequences in hairpin oligonucleotides[END_REF], while ribozymecatalyzed polymerization gives rates in the hour range [START_REF] Johnston | RNA-catalyzed RNA polymerization: accurate and general RNAtemplated primer extension[END_REF][START_REF] Zaher | Selection of an improved RNA polymerase ribozyme with superior extension and fidelity[END_REF], we can estimate the ratio β/α somewhere near 5-20 if polymerization is the rate-limiting step. Assuming a self-specific enzymatic replicator with a mutation rate on the order of 3% as in [START_REF] Johnston | RNA-catalyzed RNA polymerization: accurate and general RNAtemplated primer extension[END_REF], we obtain a critical length L c = 11-15 for weak specificity p = 1, and a larger value L c = 18-33 for p → ∞, because stronger specificity constraints allow longer sequences due to better localization. These values are significantly smaller than the lengths of, e.g., the 154-nucleotide specificity domain of Bacillus subtilis

RNase P [START_REF] Lilley | Structure, folding and mechanisms of ribozymes[END_REF] or the tRNA-like structures supposed to act as "genomic tags" for the replication of RNA viruses [START_REF] Weiner | Transfer RNA-like structures tag the 3' ends of genomic RNA molecules for replication -implications for the origin of protein-synthesis[END_REF]. Many of the nucleotides in these instances have a structural role, which makes them effectively redundant or neutral [START_REF] Kun | Real ribozymes suggest a relaxed error threshold[END_REF], and only a minority is actually involved in recognition. Also, recent research indicates that "stalling" of polymerization after mismatch incorporation might significantly reduce the error threshold [START_REF] Rajamani | Effect of stalling after mismatches on the error catastrophe in nonenzymatic nucleic acid replication[END_REF]. Nevertheless, our result suggests that the error threshold puts hard constraints on the information content of enzymatic replicators as well. T

1 = jk k m kj a j x b j T 2 = ijk k m kj b ji x b i x b j T 3 = - k k x b k j a j x b j T 4 = - k k x b k ij b ji x b i x b j .
(A.1)

The reduced mutation matrix m kj and the catalytic matrix b ij are given in Eqs. ( 7) and ( 8). Note that although not immediately obvious, the catalytic matrix b ij = b ji is symmetric, because the binomial in the denominator of Eq. ( 8) normalizes it to the single sequence level. We further keep in mind that n k = 0 if k < 0 or k > n if n and k are integer, so we do not need to keep track of the summation limits in the following calculations.

While it is straightforward to find T 1 = αL(a + μ(1 -2a)) and T 3 = -αaL, we concentrate first on T 4 , which reads after performing the summation over k

T 4 = -βaL ijn L -j i -j + n L j j n f s (i -j + 2n)a i+j (1 -a) 2L-(i+j) . (A.2) Replacing i = i -j + 2n and rearranging L-j i -n L j j n = L-i j-n L i
i n , we can sum over j and n, and are left with

T 4 = -βaL i L i f s (i )(1 -2a(1 -a)) L-i (2a(1 -a)) i = -βaLS(1 -2a(1 -a)). (A.3)
The term T 2 can after similar rearrangements be written as the product of two generalized Vandermonde matrices: (A.9)

T 2 = β i L i f s (i ) jk k ⎡ ⎣ μ 1 -μ 1 -μ μ ⎤ ⎦ L-k,j ⎡ ⎣ (1 -a) 2 a(1 -a)
Adding up T 1 + T 2 + T 3 + T 4 = 0 gives Eq. ( 9).

Appendix A.2. Solutions of Eq. (9)

This equation can be solved whenever S(x), defined in Eq. ( 10), assumes a simple form. For f s (d) = (1d/L) p as in Eq. (4), S(x) is a polynomial of order L, except for integer 0 < p < L, where it is of order p. A few instances for L > 2 are given by

S(x) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 -(1 -x) L , p = 0, x, p = 1, 1 L x + L-1 L x 2 , p = 2,
x L , p= ∞.

(A.10)

To find solutions a(μ) = 1/2, we write x = a(1a) and solve Eq. ( 9) for x: (A.11)

x k = x c k ≡ 1 2 L k [a k (1 -a) L-k + a L-k (1 -a) k
] is a sum of two complementary binomials, because in this case the first moment vanishes by construction. The four terms T 1 -T 4 are defined and evaluated analogously to Eq. (A.1), and after some algebra we find: where we take the negative sign and the positive root to obtain values μ near zero (values near unity imply complementary replication and give equivalent results for cross-specific replication).

T 1 = 1 2 αL L -2(L -1)(a(1 -a)(1 -2μ) 2 + μ(1 -μ)) (B.
We readily find that μ (1/2) > 0 if β > 0, which implies that the pitchfork bifurcation described through Eq. (B.7) is always subcritical.

  leading to a lengthy result explicitly given in Appendix B.1 and a corresponding expression for the bifurcation diagram (see Eqs. (B.5) and (B.7)).

  our choice Eq. (4) of the specificity function, and no bistability for weak specificity because the critical value of the coupling constant is negative (β * ≤ 0) for small p < p min = 2 + O(L -1 ). Since the mean Hamming distance k = aL is discontinuous at the error threshold only in the bistability regime β < β * (p), we argue that the observation of a discontinuous mean replication rate (mean fitness)α + βS(1 -2a(1a))depends on the specificity of second order catalysis, which generalizes the result of[START_REF] Wagner | Second order catalytic quasispecies yields discontinuous mean fitness at error threshold[END_REF].The results of the binomial closure approximation Eq. (9) are summarized in the phase diagram Fig.3, where the two critical mutation rates μ c and μc are shown as functions of the selection strength β and the specificity degree p for α = 1. The thick line denoted β * indicates the boundary of the bistability regime μ c > μc . Fig.2demonstrates that the binomial approximation is quantitatively excellent in the supercritical situation β > β * (in particular, it gives exact results for μ c and μc ), and qualitatively correct otherwise, where the values for μ c are somewhat underestimated: near the error threshold, the variance of the population distribution is considerably larger than that of a binomial.

  Fig.1(c) once we replace β → β/2. The "coupling constant" is only half as large because only one half of the population is available as enzymes for the other.

Figure 1 :

 1 Figure 1: Solutions for self-specific replication with fs(d) = (1d/L) p . Stationary Hamming distance distribution x k as function of mutation rate μ from numerical solutions to the reduced rate equations Eq. (6) (straight lines) and a stochastic simulation of the full system Eq. (1) in a population of N = 10 4 sequences of length L = 32 (dots) for α = 1, β = 5 and (a) p = 1, (b) p = 5, and (c) p = ∞ (here, the dashed line shows the error-tail approximation for x 0 ). The insets depict the average Hamming distance a = k /L to the master sequence.

Figure 2 :

 2 Figure 2: Bifurcation diagram of Eq. (6) for the normalized mean Hamming distance a = k /L with L = 8, α = 1, β = 1 and different values of the specificity degree p. Thick lines indicate stable branches, dashed lines unstable branches and thin dotted lines the results of the binomial closure approximation Eq. (9), which is barely visible for p = 0.1 and p = 1. Circles indicate critical mutation rates μc where the localized regime vanishes, crosses show values μc where the delocalized state changes stability (horizontal branches corresponding to delocalized states are drawn slightly shifted for visualization).

Figure 3 :

 3 Figure 3: Phase diagram of localization regimes in the parameter space of mutation rate μ, selection strength β and specificity degree p (log-log-log scale) obtained from Eq. (9) for L = 8 and α = 1: below the upper plane μc, a localized solution exists, while above the lower plane μc the delocalized state becomes stable. Bistability (μc > μc) is possible only for β < β * (p).

Figure 4 :

 4 Figure 4: Projected phase diagrams of localization regimes for self-specific replication with L = 8 and α = 1 obtained from Eq. (6): (a) as function of β with p = 10 fixed; (b) as function of p with β = 10 fixed. Below μc (thick line), a localized solution exists, while above μc (thick dashed line) the delocalized state becomes stable. The dotted lines denote the result for μc obtained via the binomial closure approximation. Bistability (μc > μc) is possible only for β < β * (p) (or p > p * (β)).

Figure 5 :

 5 Figure 5: Hamming distance distribution x k as function of mutation rate μ as in Fig. 1, but for cross-specific replication with α = 1, β = 10, L = 32, and (a) p = 5 or (b) p = ∞. The insets show the variance Δk 2 = k (kk ) 2 x k (in a symmetric population, k = L/2).

Figure 6 :

 6 Figure 6: Bifurcation diagram of Eq. (6) as in Fig. 2, but for cross-specific replication in terms of the population parameter a, which denotes the width of the two subpopulations and is calculated from the variance Δk 2 in the Hamming distance distribution. Parameters are L = 8, α = 1, and β = 1.

Figure 7 :

 7 Figure 7: Projected phase diagrams of localization regimes as in Fig. 4, but for cross-specific replication with L = 8 and α = 1: (a) as function of β with p = 10 fixed; (b) as function of p with β = 10 fixed.

  the L th Vandermonde matrix with parameters a, b, c, and d is defined as⎡ -i-j b j-c i-d . (A.5)This allows us to use a nice multiplication identity[START_REF] Rawlings | On the metamorphosis of vandermonde's identity[END_REF]μ(1 -2a) a(1a) (1a) 2μ(1 -2a) a(1a) a + μ(1 -2a))S(1 -2a(1a)) (A.8) -(1 -2a)a(1a)(1 -2μ)S (1 -2a(1a)) .

  2αμL β e 2μL , p = ∞.

  L -2(L -1)(a + μ(1 -2a))(1aμ(1 -2a)))C(1 -2a(1a)) + 2a(1a)(1 -2a) 2 (1 -2μ) 2 (L -1)C (1 -2a(1a)) + 2(a(1a)(1 -2a)(1 -2μ)) 2 C (1 -2a(1a)) (B.2) T 3 = -1 2 αL[L -2(L -1)a(1a)] (B.3) T 4 = -1 2 βL[L -2(L -1)a(1a)]C(1 -2a(1a)). (k)+ f c (Lk) x L-k (1x) k . Adding up T 1 + T 2 + T 3 + T 4 = 0 gives the condition (1 -2a) 2 μ(1μ)L(L -1) α + βC(1 -2a(1a)) βa(1a)(1 -2μ) 2 (L -1)C (1 -2a(1a)) + a(1a)C (1 -2a(1a)) = 0. (B.5) Most importantly, Eq. (B.5) reads for p = 1 in the specificity function f c (d) = (d/L) p :-(1 -2a) 2 L(L -1)μ(1μ)(α + β/2)

Figure

Acknowledgements

We gratefully acknowledge helpful discussions with Joachim Krug and Irene Chen and financial support by the German Excellence Initiative via the program "Nanosystems Initiative Munich" and by the Deutsche Forschungsgemeinschaft through SFB TR12.

Appendix A.1. Derivation of Eq. (9)

To obtain Eq. ( 9), we compute the first moment k k ẋk of the reduced rate equations Eq. ( 6) under the assumption that x k = x b k ≡ L k a k (1a) L-k is binomially distributed. This gives four terms For the last case, we approximated (1 -2x) L ≈ e -2xL and μ/(1 -2μ) ≈ μ for the important asymptote L 1 with μL fixed, and used Lambert's Wfunction. A more complicated expression is obtained for p = 2. The solution a = 1 2 1 ± √ 1 -4x is then easily computed, and the error thresholds μ c follow from evaluating the condition x = 1/4 (for p = 0, 1, 2), or from requiring the W-function to give a real result (for p = ∞):

Finally, it is easy to evaluate the bifurcation diagram Eq. ( 11) at a = 1/2 to get the critical mutation rate μc = μ(1/2):

For our specificity function

Note that μ c = μc for p ≤ 2.

While most of these results can be evaluated also for α = 0, the case p → ∞ is special. Here, Eq. ( 9) gives x = μ/(1-2μ) if we again approximate (1-2x) L ≈ e -2xL , hence the error threshold is μ c = 1/6, independent of β and L.

Appendix B. Cross-specific replication Appendix B.1. Derivation of an equation for a

To obtain an equation for the parameter a, we compute the second moment k k 2 ẋk of the reduced rate equations, Eq. ( 6), under the assumption that Appendix B.3. Solutions of Eq. (B.5)

The auxiliary function C(x) can be evaluated for small integer p as in Eq. (A.10). There is no solution to Eq. (B.5) except a = 1/2 for p = 1, and the expression for p = 2 is quite lengthy. For p = ∞ and α > 0, we get

which is exactly the result for the self-specific case if we replace β → β/2.

Accordingly, we get μ c = W[β/(2eα)]/(2L) for the error threshold.

Observing that C (1/2) = 0, the critical mutation rate μc = μ(1/2) is given by μc L → ∞, we conclude that μ c → 1/4 in this limit.


Alexandre Cornu

Steven Derrien

Dominique Lavenier

HLS Tools for FPGA: Faster Development with Better Performance

Keywords: ImpulseC, High Level Synthesis Tool, Hardware Accelerator, FPGA

Designing FPGA-based accelerators is a difficult and timeconsuming task which can be softened by the emergence of new generations of High Level Synthesis Tools. This paper describes how the ImpulseC C-to-hardware compiler tool has been used to develop efficient hardware for a known genomic sequence alignment algorithms and reports HLL designs performance outperforming traditional hand written optimized HDL implementations.

Introduction

FPGA density is increasing exponentially in such a way that the number of gates is approximately doubling every two years. Consequently, very complex designs can consequently be integrated into a single FPGA component, which can now be considered as high computing power accelerators. However, pushing processing into such devices leads to important development time and design reliability issues. Recently, many efforts have been done to help FPGA-targeted application designers to deal with such huge amount of resources. In particular, Electronic System Level tools provide higher levels of abstraction than traditional HDL design flow. Several High Level Languages (HLL) for modeling complex systems, and corresponding High Level Synthesis (HLS) Tools to translate HLL-designs into HDL synthesizable projects are now available.

Most of them are based on a subset of C/C++ [START_REF] Holland | Survey of C-based application mapping tools for reconfigurable computing[END_REF] generally extended with specific types or I/O capabilities. This paper focuses on ImpulseC [START_REF] Pellerin | Practical FPGA Programming in C[END_REF] and its associated design flow proposed by Impulse Accelerated Technologies. It also gives feedback in the context of high performance hardware accelerator design.

Experimentations have been conducted on a specific family of algorithms coming from bioinformatics and which are known to have highly efficient parallelization on FPGA. More specifically, in this work, genomic sequence comparison algorithms are considered. As a matter of fact, many efforts have been done to parallelize these algorithms, providing many optimized implementations which can be used as references [START_REF] Oliver | Multiple Sequence Alignment on an FPGA[END_REF][START_REF] Yilmaz | An Optimized System for Multiple Sequence Alignment[END_REF][START_REF] Allred | Smith-Waterman implementation on a FSB-FPGA module using the Intel Accelerator Abstraction Layer[END_REF][START_REF] Aung | C-Based Design Methodology for FPGA Implementation of ClustalW MSA[END_REF].

In this paper, we detail how ImpulseC has been used to quickly implement parallel systolic architectures. Even if here, only a single application is considered (genomic sequence comparison), the methodology we followed can obviously be extended to many other algorithms with efficient implementation on systolic array, or more generally, implementation on massively parallel architectures. To achieve reasonable performance, standard optimizations such as loop pipelining or process splitting need however to be done. The use of HLL provides an easy way to perform these optimizations, together with a fast development cycle time. However, to obtain high performance, designers have to bypass ImpulseC restrictions and perform complex program transformations.

The experiments presented in this paper show that an HLL-designed accelerator can outperform optimized HDL designs. This can be first achieved by rapidly exploring several architectural variations without great efforts compared to HDL specifications. Second, this is also achieved by the use of high level code transformations allowing the designer to generate code which better fit to HLS tool input.

The paper is organized as follows: the first section briefly describes the HLS tool we have used : ImpulseC and gives some background on the parallelization scheme used for our target genomic sequence comparison algorithm. Section 3 presents our design strategy. Performances are finally detailed in section 4, in terms of code transformation efficiency, hardware accelerator performance, and design process.

HLS Tool, Algorithm and Parallelization

ImpulseC

ImpulseC is a high level language based on ANSI C. It has a few restrictions, mainly on structure and pointer usage. On the other hand, it includes libraries to define constructor functions, bit-accurate data types and communication functions.

Two levels of parallelism are available: (1) coarse grain parallelism, by implementing several ImpulseC processes that can communicate through streams, signals or registers; (2) fine grain operator parallelism, within one process or one process loop, through the use of instruction pipelining and data flow parallelism.

Each process can be set as hardware process, meaning it will be hard-wired, or as software process, meaning its sequential code will be executed on a processor. Implementation of streams between hardware and software processes are managed by specific PSP (Platform Support Package). Here, two different PSP have been used : the Nios2 softcore and the XD2000i development platform.

The development environment (IDE) is called CoDevelopper. Designer can perform software simulation, generate HDL for a specific platform through the use of ImpulseC compiler, analyze ImpulseC compiler report through the Stage Master Explorer tool or generate HDL simulation testbench with the CoValidator tool.

Algorithm

The goal of genomic sequence comparison algorithms is to identify similarities between genomic sequences for discovering functional, structural or evolutionary relationships. Similarities are detected by computing sequence alignments which generally represent the main time consuming part. A score is attached to each alignment and represents the number of matches between 2 sequences. The exact method for computing alignments between 2 sequences is based on dynamic programming (DP) algorithms. These algorithms have been abundantly described in the literature and we just give the main idea. Basically, DP algorithms consist in computing a matrix of size N × M where N and M are respectively the sizes of the genomic sequences. The following recurrent equation explains how to fill the matrix:

S(i, j) = max ⎧ ⎪ ⎨ ⎪ ⎩ S(i -1, j -1) + sub(x i , y j ) S(i -1, j) -gap penalty() S(i, j -1) -gap penalty() 1 ≤ i ≤ N, 1 ≤ j ≤ M
x i and y j are respectively elements of sequences X{x 1 . . . x N } and Y {y 1 . . . y M }; sub() refers to a substitution cost. The final score is given by S(N, M ). From this basic equation many variations can be elaborated to compute different types of alignments. But, the main structure remains and the same parallelizing scheme can be applied to all.

Parallelization

Data dependencies are illustrated in figure 1. Note that computation of cells that belong to the same diagonal are independent and can be computed in the same time. This very regular and fine grained parallelization leads to a systolic array architecture: each Processing Element (PE) is affected to a sequence character X; then sequence Y crosses the array, character by character; after M steps, all N × M values have been computed. Again, a large literature is available on this topic. Interested readers can refer to the work of Lavenier [START_REF] Lavenier | Bioinformatics Applications, in Reconfigurable Computing: Accelerating Computation with Field-Programmable Gate Arrays[END_REF] for a complete overview of systolic architectures for genomic sequence comparison.

Using systolic architectures requires two separate phases: a initialization step and a computational step. The first step pushes a query sequence into the systolic array in such a way that one PE store one character. The second step pushes several sequences from a bank, character by character, to the systolic array.

Design Methodology and Architecture Description

The overall architecture is based on very simple two-processes structure (Figure 2(1)): A first ImpulseC process (Server) generates sequences; the second process (Master) computes alignments. Both processes exchange data through The Server process has been first hardwired (release 1), then implemented on the Nios2 softcore processor (release 2, 3) and, latter on, to the XtremeData Xeon processor (cf. subsection 4.1).

The initial Master process has been split into two entities: a systolic array and a controller entity (still called Master), which has been iteratively split into many simple but optimized processes. Finally, as implementing too many small independent processes was leading to a waste of resources, PEs have been clustered. This architecture evolution is detailed in the following subsections.

Systolic Array Design

Processing Element (PE) have been designed to constitute a 1-D systolic array. Each PE is described as an independent ImpulseC process and performs computation of one column of the matrix (one cell per PE per loop iteration). The Master process feed the first PE of the array and get results from the last one (Figure 2(2)). A PE has 3 pairs of stream interfaces, allowing it to read values from previous PE and forwarding values to the next one. Streams score i , score o , are carrying scores from previous PE and to next PE. Both sequence stream pairs x i , x o and y i , y o are carrying sequences, character by character. For sake of clarity, two different streams are shown. Actually, a single time-sharing stream is implemented as there are no concurrent read/write operations.

To compute one cell (current score : score C ), one PE need to have: the upper cell's score (score U ), the left cell's score (score L ), the upper-left cell's score (score UL ), and the two characters X C and Y C .

Algorithm 1 describes PE behavior. This description is the input of the Im-pulseC compiler to provide an FPGA implementation. Timing isn't specified and is automatically deduced by the ImpulseC compiler from data dependencies. Actually, real code is a little bit more complex, as it contains initialization (not represented here) and it considers more complex alignment computation.

Design Optimizations

This section explains the various optimizing steps which have been done to reach the final architecture. Speeding-Up Processes. Once the overall multi-process system has been functionally validated, each process has been separately optimized. This has been done through simple source code transformations or by adding compiler directives (pragmas) such as the CO PIPELINE pragma, which will lead to a pipelined execution on the most inner loop. For example, in Algorithm 1, both inner loops are important loops to pipeline. In addition, code transformations may be necessary to help the ImpulseC compiler to provide better efficient scheduling. In our case, both loops have a rate of 1, meaning that one loop iteration can be started every clock cycle.

Splitting Slow Processes. In some cases, processes exhibits codes which are too complex to achieve efficient schedule or pipeline. For example, Master in figure 2(2) must perform all the following tasks: reading 32-bit data on incoming streams from host; serializing these data to push them character by character into the systolic array; storing sequences that need to be compared to multiple sequences; controlling systolic array for different sequences; etc. To execute these tasks, the Master code uses several nested multi-conditional controlled loops.

One of the ImpulseC restrictions is that nested loops cannot be pipelined. While Master is waiting for the last result coming from the systolic array, it could have started to send the next sequence, as first PEs are ready to work. One way to overcome this issue is to merge these loops. However, this transformation leads to complex code and poor scheduling performance.

Then, splitting processes with loops that can be executed in parallel can help to improve the overall system efficiency. In the accelerator system, many processes have been created to simplify the Master process code (Figure 2(3)). A garbage process reads useless data from the systolic array. A switch process, configured by the Master, carries results either to the server software process or to the Master when needed for the next computation. A cache system reads alignment structure coming from the Server, stores sequences which need to be sent again, and does the sequence serialization. This splitting mechanism provides an efficient electronic system, minimizing process idle periods and generating simple code which is better suited for ImpulseC scheduler capabilities.

Minimizing Resources.

As the systolic array represent the major resource utilization, a lot of efforts have been done to minimize individual PE and stream resources. Streams are implemented as FIFOs and have been reduced to their minimal size, allowing memory block FIFOs to be advantageously replaced by register FIFOs. PE arithmetic has been limited to 14 bits which is the theoretical upper bound of score values. Some operators have also been transformed to simple ones when possible. All these optimizations can easily be specified at ImpulseC level and do not generate additional complexity in the design process.

Clustering PEs

Specifying one PE by one process provides efficient hardware (cf. section 4.2): the architecture is fully distributed. Signals (control and operand) remain local, leading to small fanout and easy synthesis. Nevertheless, it uses a lot of resources: each PE duplicates identical control logic and several intermediate streams.

Reducing the duplication of control logic (and streams) resources would enable a larger systolic array (section 4.2). A way to achieve this objective is to cluster PEs. At ImpulseC level, this is done by specifying the behavior of several PEs in a single process (Figure 2(4)). As an example, algorithm 2 shows the code of one cluster of 2 PEs. Clustering PE doesn't decrease loop rate, nor increases critical path length (max unit delay), at least for datapath. However it may increase control signals fanout as one state machine has to control more operators.

Algorithm 2. Dual Processing

Reducing the number of FIFO stages decreases the systolic array latency, but increases cluster latency. A small systolic array latency is adapted to process sequences larger than the systolic array size: computation are done in several passes and results of the last PE are needed as soon as it is possible for the next computation. On the other hand, a small PE cluster latency is preferable when the query size fit the systolic array size: the first character of a bank sequence can be pushed as soon as the cluster has processed the last character of the previous bank sequence. Depending of applications, a tradeoff between the size of the query sequence, the total number of PEs and the size of a PE cluster need to be found.

Performances

Platform Description

The platform used for our experiments is the XD2000i Development platform from XtremeData, Inc. [4]. This hardware holds a dual processor socket motherboard where one of the initial processor chip has been replaced by a XD2000i module. This In-Socket Accelerator embeds 2 Altera EP3ES260 FPGA components.

Design Variations Comparison

Clustering PEs had two goals: (1) share control logic and (2) decrease the number of inter-PEs FIFOs. Table 1 shows FPGA resource utilization (Adaptive Logic Module and M9K memory blocks) equivalent to 8 PE stages for different sizes of clusters. As the number of PE per ImpulseC process increases, resource utilization of equivalent partition of the design is reduced by 2.5, allowing the systolic array to be larger. For each size of clusters, we tried to fit as many PEs as possible on the FPGA, provided the design could sustain a frequency of 200MHz. Measured BCUPS (Billion Cells Updates Per Seconds) are also reported in table 1. This unit represents the number of matrix cells (cf. section 2.2) which are computed in one second, and is the standard unit to measure genomic sequence comparison algorithm performance.

It could be interesting to investigate if better performance could be reached by increasing the cluster size over 8, but, at this point, manually writing the code becomes too complex. This situation reveals the current tool limits.

Comparison with Hand-Written HDL Code

In [START_REF] Allred | Smith-Waterman implementation on a FSB-FPGA module using the Intel Accelerator Abstraction Layer[END_REF] the authors present a Smith & Waterman HDL-implementation on the same platform. The Smith & Waterman algorithm is one of the most popular genomic sequence comparison and is usually used as reference. In this hand-written code implementation, performance reaches 9 BCUPS. Table 2 summarizes performance of this implementation (HDL), of the first ImpulseC optimized release (HLL(1)) and of the highly-optimized ImpulseC release (HLL(4)).

Our final design, HLL(4), is 4.2x more powerful than HDL-implemented design.

State of the Art

While lots of paper show great speed-ups of genomic sequence comparison over sequential processor, it still remains very difficult to provide fair comparisons. Hardware platforms and algorithms are constantly evolving leading to huge design differences. Thus, we restrict our comparison to the well known Needleman & Wunsch and Smith & Waterman algorithms managing protein sequences and affine gap penalty. As in our implementation, these algorithms compute every cells of the matrix and don't make any optimization such as the diagonal-bandlimited optimization which limits computation near the main diagonal of the matrix.

In [START_REF] Oliver | Multiple Sequence Alignment on an FPGA[END_REF] and [START_REF] Yilmaz | An Optimized System for Multiple Sequence Alignment[END_REF], both design targets Virtex II XC2V6000 FPGA and respectively get 3.1 GCUPS and 10.54 GCUPS. XC2V6000 are now outdated FPGAs that contain only 33K slices. It would be interesting to see how the proposed architectures scale with current devices.

Two more recent designs, [START_REF] Allred | Smith-Waterman implementation on a FSB-FPGA module using the Intel Accelerator Abstraction Layer[END_REF], on XD2000i platform, and [START_REF] Zhang | Implementation of the Smith-Waterman algorithm on a reconfigurable supercomputing platform[END_REF], on XD1000 platform, exhibit respectively 15 GCUPS theoretical performance (9 GCUPS actually measured) and 25.6 GCUPS pick performance. These implementations seems currently to be the best references.

On conventional processors, Farrar [START_REF] Farrar | Striped Smith-Waterman speeds database searches six times over other SIMD implementations[END_REF] holds actually the best performance (up to 3 GCUPS) with the use of SIMD processor instruction set. A multithreaded implementation on a 4 core Intel Xeon processor (as the processor available on the XD2000i platform) would achieved 12 GCUPS. Our hardware accelerator could also benefit from multi-threaded opportunity. Currently, only half of the FPGA resources of the XD2000i platform is used. As our architecture is not bandwidth limited, using both FPGAs shouldn't alter the design efficiency, providing up to 76.8 GCUPS, or a x6 speed-up over optimized multi-threaded SIMD software.

Fast and Reliable Design Process

Designing hardware accelerator with High Level Language significantly speeds up design process. In this experiment, a first systolic array has been created from scratch within a few hours, while it could have taken days in HDL. First of all, using HLL allows designers to have a behavioral description at a higher level compared to HDL, avoiding lots of signal synchronization problems that occur when dealing with low level signals (data ready, data eos and data read for example). Instead, the use of ImpulseC stream read and write functions allows the designer to specify identical behavior in a very safe way. Inside a process, operation scheduling are deduced from ImpulseC code and from data dependencies. No other timing constraint need to be specified. With ImpulseC, changing operation order, for example, is thus quickly done, while it can be error-prone, or at least time consuming, in HDL.

Working with HLL makes implementation step faster, as shown in figure 3. Software simulation is even faster and is sometimes a sufficient test to pass to the next design process iteration, skipping HDL simulation which can be a more difficult way of debugging. In the same way, an intermediate model is generated during HDL generation, that can be explored with the Stage Master Explorer tool and from which designer can have a view of process performance (scheduling, pipeline rate,. . . ).

Finally, this fast design methodology allows designers to try several design options, or add extra features that can be simply crippling with HDL languages. For example, in the design considered in this paper, adding features such as variable sequence length management or dynamic load of substitution matrix is much easier than it could be with HDL.

Reusable and Portable Design

HLL makes design specification more understandable, easier to maintain and also easier to reuse. Using corresponding Platform Support Package to port an ImpulseC design from a platform to another one means that application designers would not have to manage architecture dependent issues. This role is devoted to API level and to the platform designers.

Unfortunately, the PSP we experimented is still in its early development stage and was not able to provide all standard features. We had to rewrite part of our software and slightly adapt hardware processes, when moving from Nios2 to the host Xeon processor. We can expect that the future PSP releases will make this transition transparent.

However, moving an accelerated application from an FPGA architecture to a very different one (with different RAM block size or a higher memory bandwidth for example), could require design modifications in order to get maximum performance. Unchanged design should at least provide functional designs and technology modification would be easier to do within HLS tool environment.

Obtaining Real Performance

One typical claim against HLS tools is that the reduction of design time comes at the price of a loss of performance, as these tools prevent designers from optimizing the details of the architecture. We believe such claims do not hold when it comes to the parallelization of high performance computing applications on FPGAs. As a matter of fact, we even believe that using HLL can actually help designer improving performance as it allows designers to abstract away all non critical implementation details (ex: optimizing the RTL datapath of a Processing Elemenent [START_REF] Allred | Smith-Waterman implementation on a FSB-FPGA module using the Intel Accelerator Abstraction Layer[END_REF]) and instead focus on performance critical system level issues (ex: storage/communication balance).

Indeed, thanks to the expressivity of HLL, designers effectively explore a much larger design space, leading them to optimal architectural solutions they would have never considered in the first place. Even though such system level optimizations could also be implemented at the RTL level using HDLs, it turns out they rarely are so. The explanation is quite straightforward : the impact of these optimizations on performance is often discovered at posteriori (or too late in the design flow) and reenginnering the whole design at the HDL level would require too much effort.

Last, thanks to the drastic reduction in design time to get an efficient implementation, designers can afford to implement several non critical design optimizations. Even if the impact of one of such optimization alone is generally limited (often 10%-30% performance gain), their combined effect on global performance is generally multiplicative/cumulative and ultimately boosts global performance, for very little design effort.

Toward User Assisted Source to Source Transformation Tools

It is clear that HLS do not provide a direct algorithmic to gate design path. In a similar way high performance architecture programmer do, designers need to modify and rewrite their source code so as to get the tool to derive the hardware architecture they actually want for their application. Even though such a rewrite process is much easier to carry with a HLL, it still leaves a lot of burden to the designer.

As a consequence, we feel that there is a need for HLS oriented source to source refactoring tools, which would help designers re-structuring their programs by providing them with a code transformation toolbox tailored to their needs. Such a toolbox would of course encompass most of the transformations that can be found in semi automatic parallelizing source to source compilers (loop and array layout transformations for example), but should also provide domain-specific code transformations specifically targeted at HLS tools users (process clustering could be one of them).

Fig. 1 .

 1 Fig. 1. Data Dependency for Similarity Matrix Computation

Fig. 2 .

 2 Fig. 2. Design Releases

Fig. 3 .

 3 Fig. 3. ImpulseC design flow

Algorithm 1 .

 1 Simple Processing Element Pseudo-Code

	ProcessingElement(xI , xO, yI , yO, scoreI , scoreO)
	loop
	xC ← stream read(xI )
	while xTMP = EOS do
	xTMP ← stream read(xI )
	stream write(xTMP ) → xO
	end while
	for n = 1 to N do
	yC ← stream read(yI )
	stream write(yC ) → yO
	scoreL ← stream read(scoreI )
	scoreC ← compute score( scoreL, scoreUL, scoreU , , xC, yC )
	stream write(scoreC ) → scoreO
	scoreU ← scoreC
	scoreUL ← scoreL
	end for
	end loop

  Element Cluster Pseudo-Code

	ProcessingElement(xI , xO, yI , yO, scoreI , scoreO)
	loop
	xC1 ← stream read(xI )
	xC2 ← stream read(xI )
	while xTMP = EOS do
	xTMP ← stream read(xI )
	stream write(xTMP ) → xO
	end while
	for n = 1 to N do
	#pragma CO PIPELINE
	yC ← stream read(yI )
	scoreL ← stream read(scoreI )
	scoreC1 ← compute score( scoreL1, scoreUL1, scoreU1, xC1, yC1)
	scoreU1 ← scoreC1
	scoreUL1 ← scoreL1
	scoreL2 ← scoreC1
	yC2 ← yC1
	scoreC2 ← compute score( scoreL2, scoreUL2, scoreU2, xC2, yC2)
	scoreU2 ← scoreC2
	scoreUL2 ← scoreL2
	stream write(yC2) → yO
	stream write(scoreC2) → scoreO
	end for
	end loop

Table 1 .

 1 Impact of clustering PEs on Resources and Performances

	Design	HLL(1) HLL(2) HLL(3) HLL(4)
	Cluster size (x)	1	2	4	8
	8 PEs ALM util.	4336	3680	2246	1578
	8 PEs M9K util.	8	4	4	4
	Design Logic Util.	88 % 88 % 74 % 81 %
	Max number of PE (N) 224	256	288	384
	BCUPS	20.5	23.7	30.5	38.4

Table 2 .

 2 HLL versus HDL design performance

	Design	HDL HLL(1) HLL(4)
	Systolic Array Size 300	224	384
	BCUPS	9	20.5	38.4

Acknowledgment

We would like to thank Impulse Accelerated Technologies for their valuable help. This work has been supported by the French ANR BioWIC (ANR-08-SEGI-005).