Kinetic Selection of Template Polymer with Complex Sequences

Emergence and maintenance of polymers with complex sequences is a major question in the study of origins of life. To answer this, we studied a model polymerization reaction, where polymers are synthesized by stepwise ligation from two types of monomers, catalyzed by a long polymer as a template. Direct stochastic simulation and dynamical systems analysis revealed that the most dominant polymer sequence in a population successively changes against the flow rate of monomer to the system. The slower the flow, the more is the complex sequence selected. We discuss the relevance of this kinetic selection of sequence by the non-equilibrium flow rate to the origin of complex polymers.

The emergence and maintenance of polymers with complex sequences pose a major question in the study of the origin of life. To answer this, we study a model polymerization reaction, where polymers are synthesized by stepwise ligation from two types of monomers, catalyzed by a long polymer as a template. Direct stochastic simulation and dynamical systems analysis reveal that the most dominant polymer sequence in a population successively changes depending on the flow rate of monomers to the system, with more complex sequences selected at a lower flow rate. We discuss the relevance of this kinetic sequence selection through nonequilibrium flow to the origin of complex polymers.
Schrödinger, in his celebrated monograph "What is Life" [1], recognized that the essence of a heredity-carrier lies in aperiodic crystals, as was soon confirmed by the elucidation of DNA structure. In general, biopolymers consist of different kinds of monomers, constituting an aperiodic complex sequence, as is important not only for genetic information, but also for the catalytic function.
With regard to the origin of life or self-replication, catalytic polymers are required, as has also been investigated in recent experiments [2][3][4][5][6]. In their study, polymers with complex sequences had to be synthesized and preserved, to encode a large amount of information and for catalytic functions [7,8]. Hence, it is important to uncover the condition that allows for the generation of a polymer with a complex monomer sequence.
If there exist multiple catalytic polymers with sufficient lengths, then the autocatalytic process for their synthesis can select one of such polymers. Here we study the dependence of the selected sequence on the flow rate of monomers. The replication of polymers needs a supply of monomers to compensate polymer degradation. Such a nonequilibrium open condition could be provided by a hydrothermal vent at the origin of life and by a chemostat condition in a laboratory experiment [27]. Despite the recognition that such nonequilibrium flow is essential to the emergence and maintenance of life, its influence on the selection of sequence has not been systematically investigated.
The selection of a specific polymer by catalytic reactions studied so far [8,[16][17][18][19][20][21][22][23] does not show any dependence on external conditions. Ref [24] reported the dependence of the selected polymer sequence on the bias in the monomer components. In contrast, here, we study a template-catalyst polymerization model by a complementary sequence, and demonstrate that the selected sequence of the template polymer depends critically on the flow rate, and more complex sequences are selected as the flow rate is decreased. This selection mechanism is the kinetic, rather than energetic; selection of complex sequences (that include a variety of subsequences, which are defined later) without any specific design of energy dependence, as a result of changing the rate of external supply of monomers. As will be analyzed by dynamical systems and combinatorial analysis, complex sequences that are synthesized via diverse kinetic pathways from monomers are selected for a low flow rate. The generality of this selection mechanism will be discussed with possible relevance to prebiotic polymer synthesis.
We consider a synthetic reaction of polymers, which consist of two kinds of monomers, denoted as 0 and 1, where each polymerization is catalyzed by a long template polymer. Thus, all polymer species are described by the binary integer s of length l (for example, 01, 0101, 00000 etc.). The polymerization progresses in a container with volume V under a well-mixed (homogeneous) condition. Only the two monomer species flow into the container from outside at the same constant rate 1 2 f V , while all the molecules in the container diffuse out at rate d.
For simplicity, only the ligation reaction between the polymers and monomers is considered here. The polymer has directionality, and the monomer can ligate from both the left and right sides of a polymer (e.g., 1 + 00 → 100 and 00 + 1 → 001 ). For simplicity and to focus on the kinetic aspect, all the polymer species of the same length are assumed to have the same chemical potential. Thus, in equilibrium, the concentrations of all the polymer species with the same length are equal.
A template polymer serving as a catalyst accelerates both the forward and backward reactions, without changing the equilibrium condition. The catalytic reaction works if the template includes a subsequence of the bit inversion of the product. For example, the reaction 111 + 1 → 1111 is catalyzed by the template containing 0000 as a subsequence, i.e., 00000, 10000, or 00001. For simplicity, only the longest polymers with the given length L, independent of their sequence, can serve as the template in ligation reactions. Since all the ligation reactions do not change the total number of monomers, the ratio of monomer inflow to the outflow by the diffusion of the polymers determines the total number of monomers in the system, i.e., the total monomer concentration is given by f /d.
In summary, the chemical reactions are given as follows by denoting the polymer of the length l with the sequence s = {0, 1} l as A l,s : where m is a monomer 0 or 1 and sm and ms are the sequences with m added to the left side or right side of s, respectively. The polymerization can occur spontaneously at a rate ǫ, which is set to be quite small.
Under the large volume limit, the concentration of polymer A l,s , x l,s follows the deterministic rate equation.
Here, x l−1,s and x (2) l−1,s are the concentrations of the (l − 1)-mer, which is a subsequence of A l,s on both sides [30], κ l,s = 1 2 (ǫ + (l ′ ,s ′ )∈T l+1,s n l ′ ,s ′ /V ), and the T l,s is the set of all template molecules that can catalyze the reactions to produce the polymer species A l,s . This equation has multiple attractors: As the pair of the longest polymer of a given sequence and its complementary pair catalyze the reaction to synthesize them, the concentrated population for each of the complementary pairs can be an attractor. It should be noted that since the complementary sequence works as a template for the polymerization, the complementary pair (e.g., 0010 and 1101) always coexists in an equal fraction.
If the system size is finite, fluctuations around the rate Eq. 1 induce switching among the attractors, and specific attractor(s) are selected dominantly. We investigated this selection by numerically solving the original stochastic reaction model [31]. In Fig. 1(a), the time series of the concentration of template species is plotted. A specific complementary pair of sequences is dominant for some time span, followed by switching to a state with a different dominant pair. As the concentrations of complementary sequences are equal, only one of them is plotted throughout the paper. Next, the temporal average of the concentration of each template sequence was computed to examine its dependence on the flow rate, f ( Fig. 1(b)), by fixing f /d = 1 (hence, the total monomer concentration f /d is kept constant). As shown, the concentration of each sequence, and in particular, that of the dominant species depend on f . Note that a sequence and its reverse form (e.g. 0010(1101) and 0100(1011) ) have the same fraction over a long time average due to symmetry, although at each instance, one of them is selected in the system. For L = 4, 0101 (and 1010) is selected for a large f , while 0011 is selected for a small f . For L = 6, as is shown in Supplemental Material[32] Fig. S1(a), the dominant sequence changes as 010101 → 011001 → 001001 → 000100 with a decrease in f .
These changes in the dominant sequence against the decrease in the flow rate are interpreted as the increase in the sequence complexity. Here, a "complex" sequence implies that it contains more subsequences including their complementary ones. For example, 0101 includes only 01 and 10 dimers, while 0010 includes 00, 01, 10, and 11, and is thus more complex. In general, let C l (2k) be the set of template sequences that have 2k species of l-mers as subsequences (note that the complementary sequence is included in the counting). For example, 0101 belongs to C 2 (2) because it includes two dimer subsequences 01 and 10, and 0110 belongs to C 2 (4) as it includes 01, 11, 10, and 00. To study the selection of a complex sequence with f , we computed the sum of the average residence time at each attractor corresponding to the sequences belonging to C l (2k). It is defined by x l (2k) = s∈C l (2k) x L,s . In Fig. S1(b) in the Supplemental Material[32], x l (2k) is plotted as a function of f . For L = 4, the dominant concentrations change from x 2 (4) to x 2 (2), and for L = 6, they change in the order x 3 (8), x 3 (6), x 3 (4), x 3 (2) with an increase in f , implying that a more complex sequence is selected for a smaller f .
To understand this change in the dominant sequence, we note that the residence time for each attractor, under noise due to the finiteness in the molecule number, is larger if the attraction toward it is stronger. To understand this change in the dominant sequence, we first study the change from the sequence with C 2 (4) to that with C 2 (2). Since the polymerization progresses from a dimer to a trimer, and then to a tetramer, the transition first occurs between tetramers that share the trimers, i.e., those with only a one-bit difference. The transition diagram is shown in Fig. 2(a), where the change from C 2 (2) with 0101 to C 2 (4) with a decrease in f is mediated by the transition to sequence 1101 (0010). Hence, we first discuss the competition between the attractor with the dominance of the 0101 and 1010 pair (denoted as B) and that with the 0010 and 1101 pair (as A) by focusing only on these two pairs of templates, while neglecting other template polymers [36]. By assuming that the number of shorter molecules changes faster than do the templates, the reduced rate equation is written only in terms of the concentrations of these templates, as follows: Here, r A (r B ) is the synthesis rate of the corresponding template for each type A (B) for a given x A , x B , respectively, whose form is obtained from Eq. S3, as shown in Supplemental Material [32]. This rate equation has two fixed-point attractors, . The vector field for dx A /dt, dx B /dt is thus obtained as in Fig. 3. Here, the attractor close to the unstable fixed point is each kicked out by noise.
It is shown straightforwardly that the attraction speed to attractor A is given by , from Eq. S1, S2 and the eigenvalues of the Jacobian matrix of the dynamics, and that to ). As shown therein v B − v A < 0 for f ∼ 0 and > 0 for f ≫ 1, explaining the transition from B to A by noise with the decrease in f . This flow rate dependence is also interpreted by the diversity in reaction pathways to replicate templates A and B. While only one reaction pathway exists for synthesizing B from the monomers with double the speed because of symmetry, various pathways exist for synthesis of A from the monomers because they have more types of subsequences. For large f (= d), double autocatalytic paths to B work efficiently, while existence of diverse subsequences for the pathways to A causes a larger loss of shorter polymers, leading to r B > r A . For small f , the concentrated use of the same dimer for B results in its deficiency, and r A is larger since A can use The above argument holds for all the template sequences belonging to C 2 (2) and C 2 (4). Therefore, the total concentration of the sequence x 2 (2)(x 2 (4)) is large when f is large(small). Furthermore, it is valid for the L-mer template sequence A belonging to C l (2k + 2) and B belonging to C l (2k) sharing subsequences in part. Here again, the replication speeds of A and B, i.e., r A and r B , determine the attraction speed to the fixed points A and B, and the attractor A with complex sequence is dominant for a low f . In Fig. 4, the ligation reactions from the monomer to 8-mer templates are drawn. The most complex sequence with C 4 (10) has a variety of pathways, as compared with the simple sequence with C 4 (2). For a low f , the sequence with C 4 (10) is expected to be selected, while for a high f , that with C 4 (2) is selected. The Reaction pathway to synthesize 8-mer templates from monomers for (a) a complex sequence (01101110 ∈ C4(10)) and (b) a simple sequence (01010101 ∈ C4 (2)). Black and white circles represent monomers 1 and 0, respectively. Each arrow indicates ligation reaction between a shorter polymer and a monomer, which are catalyzed by product templates.  (2)), 00110011 (C4(4)), 01011010 (C4(6)), to 01001110 or 01101110 (C4 (10)) with a decrease in flow rate.
fraction of each sequence obtained from direct simulation is plotted in Fig. 5 as a function of f . Since direct simulation of cases with larger L is timeconsuming, we studied a reduced model, where only the transition between the attractors with a certain dominant template was considered, with the transition probability determined by the difference between the replication rates of the two templates (see Supplemental Material[32] for details). A complex sequence with a large number of subsequences was selected as the flow rate decreased.
We have confirmed the generality of the present result, i.e., successive increase in the complexity of dominant sequence with the decrease in the flow rate, for the following cases (see Supplemental Material[32]): (i) inclusion of variation in the catalytic activity by each template, if it is not too large (ii) change in the rules for catalytic strength (iii) inclusion of the ligation between polymers longer than monomers (iv) difference in the concentrations of monomers (0,1): In the study here, complementary polymers take the same concentrations, irrespective of the monomer concentrations, and thus, the bias in monomer cannot select one of the complementary molecules. (v) increase in spontaneous ligation rate ǫ. Finally, in the study presented here, the number of monomer types is only 2, in contrast to 4 in RNAs and 20 amino acids in proteins. Still, the increase in polymer complexity by the diversity of reaction pathways holds for a larger number of monomer types (see also [38]).
To sum up, we have shown that in a polymerization process with template-based catalytic reactions, the dominant sequence selected depends on the flow rate of the monomers, and its complexity increases with the decrease in the rate. The selection is based on two basic mechanisms; (I) the attractor selection of a higher growth rate due to stochasticity in reactions, and (II) preference of a complex sequence to have a variety of pathways to avoid the deficiency of sub-part shorter polymers.
(I) The finiteness of the system size provides stochasticity (noise) in reaction, which gives rise to switches among attractors with different dominant sequences. Attractors with a higher growth rate have larger stability against noise (see for example Fig. 3), and are selected dominantly under the presence of noise. If the noise strength is too small, no transitions occur among the attractors, whereas if it is too large, no specific sequence dominates[39] (see also [40,41] for attractor selection of a higher cellular growth, and [14,42,43] for relevance of noise to chemical evolution).
(II) When the flow rate is limited, the polymerization of simple (say 01-periodic) sequence suffers from the deficiency of subsequence polymers of shorter length. The polymerization of a complex sequence of diverse pathways can keep a variety of subsequence polymers (see also [44] for an increase in component diversity in the protocell under low nutrient flow), while under a larger d value, the decomposition of many components is disadvantageous for the growth.
This variety of subsequence leads to the definition of the complexity in sequence we adopted here [45]. Of course, this is one possible definition, whereas a related definition is adopted to compare the subsequences in the genome in bioinformatics [46][47][48][49][50][51]. The relationship between sequence complexity and structural or functional property of polymers has been discussed in several studies [52]. In our study, such a complex sequence is kinetically selected in the low flow rate region.
The results presented here will provide a novel perspective on kinetic conditions for the origin of life. Although a nonequilibrium condition is needed for life, our result implies that an excessively strong nonequilibrium condition will damage the sequence complexity. For the emergence and maintenance of prebiotic autocatalytic systems, op-timal nonequilibrium conditions are required. Although it is difficult to estimate the actual flow rate for the appearance of complex sequences, the flux rate dependence of the complexity of synthesized polymers will be important to uncover the possible condition for the origin of life and for the laboratory construction of prebiotic selfreplication systems.

Stationary distribution of the case with hexamer
We calculated the average concentration of the template species following the same procedure as that for the tetramer case shown in Fig. 1(b), in the main text. The result showed a stepwise switch of the dominant sequence species. Note that as the flow rate reaches zero, the system approaches thermal equilibrium, where the frequency of each template sequence tends to be uniform (by considering also the backward reaction). The number of sequences that belong to C 3 (6) is the largest among all the hexamers, and therefore, in Fig. S1(b), the frequency x 3 (6) increases as the flow rate approaches zero. We set L = 6, V = 500, ǫ = 5.0 × 10 −6 . The most dominant sequence is switched stepwise with the decrease in flow rate, as 010101 → 011001 → 001001 → 001000, which correspond to C3(2), C3(4), C3(6), and C3 (8), respectively. For the flow rate 10 −7 , the distribution approaches that of thermal equilibrium.

Derivation of dynamics of two templates
Considering the switch between the dominant sequences of 0010(1101) and 0101(1010), as mentioned in the main text, the rate equation is reduced to the concentrations of these templates and the corresponding trimers, dimers, and monomers.ẋ Here, x m represents the concentration of 0 or 1, x a is 00 or 11, x b is 01 or 10, x a ′ is 001 or 110, x b ′ is 010 or 101 , x A is 0010 or 1101, and x B is 0101 or 1010. Note that the respective complementary pair has the same concentration. The reaction coefficient in the above equations is proportional to the number of corresponding reactions. For example, 00 (a) is produced in a single reaction 0 + 0 + 0010, and is consumed in the reaction 00 + 1 + 0010 by the template 0010 (A). On the other hand, 01 (b) is produced in four reactions 0 + 1 + 0010, 0 + 1 + 1101, 0 + 1 + 0101, and 0 + 1 + 1010 (see main text Fig. 2(b)), where B and A are used in two reactions, respectively, and consumed in reactions 0 + 01 + 0010, 01 + 0 + 0010, 1 + 01 + 1101, 01 + 0 + 0101, 01 + 0 + 1010, 1 + 01 + 0101, and 1 + 01 + 1010, with three reactions using A and four reactions using B.
The concentrations of a, b, a ′ , b ′ and m are adiabatically eliminated, and two-variable dynamics just for x A and x B are obtained as Eq. S1 and S2 with Using them, the explicit form of r A and r B can be written as (Note that x m is also the function of x A and x B but cannot be solved explicitly.)

Derivation of strength of attraction near the fixed point
The Jacobian of the rate equation at a fixed point corresponding to the A-abundant attractor is given by Here, r A and r B are the replication rates of A and B, respectively, and The eigenvalues of this Jacobian are ∂rA ∂xA (x * A , 0)x * A and r B (x * A , 0) − d, and the maximum eigenvalue of the Jacobian is given by Therefore, the replication rates of A and B, i.e., r A and r B , determine the maximal eigenvalue of the Jacobian.

Flow rate dependence of concentration of shorter polymers and relation to reaction pathways
For a large f (= d), x a = x A x 2 m /(2d) and x b = (2x A +2x B )x 2 m /(2d), so that x a < x b . The subsequence b is produced to a greater extent because b is shared by both A and B, while the consumption of b is irrelevant in the fast dilution case (in Eq. S3, the dilution term d in the denominator is much larger than the other term). If the supply of dimers is fast, the replication speed of the simple sequence B is higher, as the reactions are concentrated on those using b. Using this, the replication rates r A and r B become The attraction speed near the fixed point A is given by 2 ) < 0 and that near the fixed point B is given by On the other hand, when f (= m 0 d) is small, x a = x m and x b = 2xA+2xB 3xA+4xB x m , so that x a > x b . The replication rates are (S7) When f is small, the attraction speed toward the fixed point A is given by r A (x * A , 0) − r B (x * A , 0)(= 8 3 x 2 m ) > 0 and that toward the fixed point B is given by r B (0, x * B ) − r A (0, x * B )(= − 1 2 x 2 m ) < 0. In summary, the attraction to the fixed point A is faster (slower) than that to B when f is small (large).

Reduced description of the model
Since the original model requires excessive computational time, we examined the selection in a system size with larger L using the reduced description. Here, the reduction scheme is as follows.
In the reduced model, we consider only the transition between the attractors, in which a certain pair of templates is dominant. During the transition, we assume that the number of templates other than the pair is zero. This is guaranteed when the spontaneous ligation rate, ǫ, is small because a new template will appear only by a spontaneous reaction. For template length L, there are 2 L−1 attractors. We estimate the transition probabilities, T A→B , between A-abundant and B-abundant attractors by noise.
First, we consider the dynamics of concentration of a pair of templates A and B. By assuming that the dynamics of concentration of the shorter polymer is sufficiently faster than that for the template, the concentrations of polymers shorter than L − 1 are adiabatically eliminated and determined as functions of the concentration of A and B, i.e., x A and x B , so thatẋ This equation has two fixed points Then, the corresponding Langevin equations are given bẏ where We assume a one-dimensional slow manifold, in which the dynamics connecting the two attractors are restricted on the straight line between the fixed points A and B. We also assume that the replication rate, r A , and degradation rate, d, are constant along the line, so that G A = 2x A d and G B = 2x B d.
Since we assume that the dynamics are restricted to the one-dimensional line, the probability that the system escapes from attractor A and enters B is calculated by the Kramers formulae, as follows (for example, [1,2]): where the integral x is carried over from fixed point A to point x along the straight line connecting the fixed points A and B, and x is the value at which the integral takes the maximum value. Using the approximation x * A ∼ x * B , we can now obtain the transition probability explicitly.
∆ A→B is considered to be the highest barrier in the quasi-potential of the dynamics of x A and x B along the line (To calculate r A and r B explicitly, we also used x * A ∼ f /(2dL) by assuming that almost all the monomers (f /d) are included in the templates, and the free monomer concentration x m ∼ √ f , as derived from Eq. 1 presented in the main text, by assuming that the concentrations of polymers other than the templates are sufficiently low.) We calculated all the possible transitions between the attractors and constructed the transition matrix T . We assumed that transition occurs only between the pairs of L-mer templates that share the (L − 1)-mer subsequences. By calculating the eigenvector with the maximum eigenvalue of the transition matrix T , we could obtain the stationary distribution of the model.
We calculated the residence probability of the attractors in the stationary distribution for L = 4, 8, 12, as shown in Fig. S2, derived from the reduced model.
Although the obtained concentration of the template species differs quantitatively from that obtained by the original chemical reaction model for a small L, the stepwise increase in complexity in the dominant sequence with a decrease in flow rate was reproduced. In other words, at the highest flow rate, the periodic sequence (010101..., C 2 (2)) is selected. Here, we also examined the case with a longer L, e.g. L = 12. As shown in Fig. S2(c), an increase in complexity was seen with a decrease in f .

(i) Variation in catalytic activity
We tested the case in which the catalytic activity of each template L for the sequence s is not identical, given by k L,s , so that κ l,s in the main text is replaced by 1 2 (ǫ + (L,s ′ )∈T l+1,s k L,s ′ n l ′ ,s ′ /V ). As an example, we studied the tetramer as the template, by changing the catalytic activity of the template 0101 (and 1010), which is the dominant sequence for fast flow, so that the activity for the sequence is given by k while that for other sequences is set at unity.
In Fig. S3, we plotted the fraction of each sequence, as shown in Fig. 1b in the main text. The dominant templates change depending on the flow rate, whereas the transition flow rate and the frequency of dominant species changes so that the fraction of sequence 0101 is increased (decreased) with the increase in k. If k is too large (say ∼ 2) ( too small (say ∼ 0.5) ), 0101 remains (is not) the dominant species, respectively, against the change in the flow rate. The competition between the difference in the catalytic activities and sequence diversity can be compared to the competition between the energy and entropy in equilibrium thermodynamics.

(ii)Change in the ligation speed
To examine the generality of the results in the main text, i.e., the increase in sequence complexity with the decrease in flow rate, we studied an alternative model in which the rules for catalytic strength were changed from that in the original model.
In the alternative chemical reaction system, the ligation reaction rate of l-mer polymer with sequence s (A l,s ) κ l,s is changed as follows. where m l,s,L,s ′ is the number of l-mers with sequence s, which is contained in the template with sequence s ′ . Under this rule, for example, the tetramer sequence 0101 has two dimers 01 as subsequences and catalyzes the reaction 0 + 1 → 01 three times faster and 1 + 0 → 10 two times faster as compared to the catalytic strength in the original rule.
We calculated the flow rate dependence of the average concentrations of templates in the stationary state (Fig. S3). Although the dominant sequence at each flow rate differs from that of the original rule (e.g. in the high flow rate case, not only 0101 but also 0000 is selected), a stepwise increase in the complexity of the dominant sequence with a decrease in flow rate is observed, as in the original model.

(iii)Ligation between polymers beyond the monomer
We also simulated another model in which there is not only ligation between a monomer and a polymer but also ligation between two polymers longer than the monomers (e.g. 01 + 01 → 0101, 01 + 0101 → 010101), so that the ligation reaction in the original model is replaced by In this case also, the flow rate dependence of the complexity of the dominant sequence is observed, as in Fig.S5.
(iv) Difference between concentration of monomers 0 and 1 In the main text, we assumed that the ratio between monomer species 0 and 1 was equal. Here, we also calculated the case in which the supply rates of the two monomers are different, by varying the flow rate of monomers 0 given by f 0 , and that of monomers 1 given by f 1 , while keeping f 0 + f 1 = 1 fixed. The flow rate dependence of the averaged concentration of template species is shown in Fig. S6. If the bias between the two flow rates is not too large, qualitatively, the same result is obtained as in the main text Fig. 1b. This is because in our model, the templates replicate as a complementary pair between 0 and 1, in contrast with the model in ref. [3]. The same number of monomers (0 and 1) is consumed to make each complementary pair of polymers, so that complementary polymers (e.g., 0000 and 1111) take the same concentrations, irrespective, of the monomer concentrations. Thus, the bias in monomer cannot select one of the complementary molecules. As the bias increases, however, the synthesis of templates is slowed down because the minority monomer (i.e. monomer 1) is deficient. Then, the complementary pair hardly exists (i.e. an attractor at which a specific pair of templates exists is unstable) and the catalytic reaction stops. In this case, only the polymer consisting of the majority monomer (i.e. 0000) is synthesized by spontaneous (template independent) ligation reactions.
(v) Change in the spontaneous ligation rate ǫ We also studied dependence on the spontaneous ligation rate ǫ. The flow rate dependence of the average concentration of template sequences is plotted in Fig. S7, as in the main text Fig. 1b. As shown, the result in Fig. 1b ( i.e. the dominance of sequence 0101 in the fast flow region, and 0011 in the slow flow region.) is robust against the increase in ǫ up to ∼ 2.0 × 10 −3 . If ǫ is beyond this threshold, the attractor with a dominance of a certain template disappears and the template species just coexist (similar to the error catastrophe [4]), where the flow rate dependence is not observed.