Symmetry-breaking in chiral polymerisation

We propose a model for chiral polymerisation and investigate its symmetric and asymmetric solutions. The model has a source species which decays into left- and right-handed types of monomer, each of which can polymerise to form homochiral chains; these chains are susceptible to `poisoning' by the opposite handed monomer. Homochiral polymers are assumed to influence the proportion of each type of monomer formed from the precursor. We show that for certain parameter values a positive feedback mechanism makes the symmetric steady-state solution unstable. The kinetics of polymer formation are then analysed in the case where the system starts from zero concentrations of monomer and chains. We show that following a long induction time, extremely large concentrations of polymers are formed for a short time, during this time an asymmetry introduced into the system by a random external perturbation may be massively amplified. The system then approaches one of the steady-state solutions described above.


Introduction
Studies of the origins of life raise many associated fundamental questions. Among these, one is concerned with the origin and propagation of molecular handedness. It is well known that chirality is a signature of life as we know it. Nucleic acids contain only D-sugars while proteins are made only from L-amino acids (although D-amino acides do occur in Nature and even occasionally show up in some proteins (Jung, 1992). What leads to the synthesis of homochiral polymers, in which all the constituent monomers have the same handedness ? And what is responsible for the evolution of chiral purity, the more or less exclusive dominance of one macromolecular handedness over its mirror image? These are questions of great interest and importance and remain the subject of much discussion.
It is known that, in general, the addition of the correct enantiomer to a growing polymer chain is more favourable that the wrong one (Joshi et al. 2000). Indeed Joyce et al.(1984) showed that addition of the wrong-handed monomer to a growing oligonucleotide chain acts as a chain terminator, stopping all further reaction. For the case of proteins, there is also the driving force of the beneficial secondary structures, such as α-helices and β-sheets, that may arise from homochiral polymers. Given these assumptions, that wrong-handed monomers inhibit chain growth, our paper is concerned with whether, starting from a racemic mixture of monomers, it is possible to produce a system of homochiral polymers of a greater or lesser degree of chiral purity. Starting from an achiral substrate, we shall be concerned with whether it is possible to produce a system of homochiral polymers of high chiral purity by analysing some plausible kinetic models.
There is some discussion of related matters and experimental observations in the recent literature. Zubay (2000) provides a readable discussion of possible pre-biotic chemistry, while Colonna et al. (1994) describes a number of alternative self-reproducing systems. Sandars (2002) reviews the range and importance of chirality in biological systems, as well as the chemical processes which lead to achiral states. After summarising a possible historical order of events in the origin of life on Earth Sandars discusses the stages where a bifurcation to a chiral state may occur. He then applies existing knowledge of chirality on Earth to speculate on the question of extra-terrestrial life and its chirality.
Luisi's group at ETH in Zurich has studied various polymeric systems in which left-and right-handed monomers aggregate together to form larger than expected concentrations of homochiral polymers Hitz et al.2001). Hitz et al. present an analysis of the data and postulate that the exess of homochiral polymers is due to a high-order Markov process rather than the feedback mechanism which we analyse in this paper. It is readily understood that the rate coefficient governing polymer growth may depend on the handedness of the monomer to be attached to the chain and handedness of the monomer which currently terminates the chain; however, with a high-order Markov processes this rate coefficient may also depend on the handedness of the penultimate (and possibly the antepenultimate) monomers in the polymeric chain.
While oligopeptides spontaneously form homochiral sequences, Hitz & Luisi (2002) have shown that the presence of quartz promotes the production of a high yield of homochiral sequences. More recently (Hitz & Luisi, 2003) this has been quantified and the level of enantiomeric exess in a system has been amplified from 20% to over 70% and in some cases to 100%, by the presence of quartz.
The model of symmetry-breaking in chiral polymerisation which we explore in this paper is based on a model suggested by Sandars (2003). Sandars gives an account of the history of chemical discoveries leading up to the mechanisms of enantiomeric cross-inhibition and autocatalysis upon which his and our models rely. Sandars integrates the resulting system of ordinary differential equations numerically to explore the parameter regimes where symmetrybreaking solutions exist. He observes a 'phase-transition' type of phenomenon where a small change in the fidelity of the feedback mechanism leads to a large-scale change in the steady-state which the system as a whole converges to. Below a critical fidelity in the nonlinear feedback process the system approaches a symmetric state where equal amounts of left-and right-handed polymers coexist, whilst above the critical fidelity a homochiral state is approached in which one chirality of polymer dominates to the almost complete exclusion of the other.
The model studied here includes the inhibition of homochiral sequences of long chains. This bears some similarity with our modelling of cement hydration (Wattis & Coveney, 1997) in which larger clusters were susceptible to poisoning by another component. Here the poison is simply the monomer of opposite chirality, so each monomer can either act as an agent of growth (of polymers of the same handedness) or of inhibition (of polymers of the opposite handedness). This dual role leads to some subtle and interesting effects since an abundance of, say, right-handed polymers makes it unlikely for any left-handed polymeric sequences to form, and the majority of left-handed monomers produced will be consumed by inhibiting right-handed homochiral sequences -a 'double-whammy' effect. This form of competition is distinct to the models of nucleation involving competition we have analysed previously; for example in Wattis (1999) and in Bolton & Wattis (2004) there is only one monomer which assembles to form two morphologies of cluster. Competition is thus between the growth of one type of cluster and that of another. However, similar mechanisms are operative in both those examples and in the present paper, since both types of homochiral polymer sequence (left and right) are ultimately composed from, and hence competing for, the same source material.
We have investigated the growth of RNA chains in an earlier paper, (Wattis & Coveney, 1999) wherein we used a much more detailed and hence complicated model to assess the feasibility of long self-replicating RNA sequences forming within a realistic timescale. While that model of RNA polymerisation had no precursor species, it contained four types of nucleotide monomer (A, C, G and T), both autocatalytic and cross-catalytic polymerisation mechanisms and an important hydrolysis step that recycled growing RNA sequences. On the basis of plausible assumptions about the prebiotic soup, and by invoking a number of approximations and coarse-graining procedures, we were able to show that self-replicating RNA sequences are amplified in such mixtures, while less capable replicators are driven to extinction.
By comparison, in the present paper we are able to construct a model which is more directly tractable using standard methods of mathematical analysis. The paper is structured as follows. In section 2, we specify the basic kinetic model, while in section 3 we study its steadystate solutions. Section 4 considers the time-dependent achiral solutions and investigates their kinetic stability. Section 5 then considers the case of perfect chiral symmetry-breaking. It is followed by a discussion (section 6) and conclusions from our work are presented in section 7.

The kinetic model 2.1 Microscopic modelling
We aim to investigate spontaneous symmetry breaking in a system which allows both rightand left-handed chiral polymers to form. We assume there is some achiral source S which spontaneously transforms into right-and left-handed monomers at a slow rate ε; we also assume that the presence of longer chiral polymers (denoted L or R) accelerates the formation of monomers of the same chirality. Thus we shall be concerned with the following set of coupled chemical reactions: where L 1 , R 1 are the left and right monomer species respectively while Q (P ) represents some measure of the total concentrations of left-handed (right-handed) homochiral polymers in the system. The precise forms of these rate processes will be specified later (see equation (4)).
The rate coefficients of the reactions are ε, k. We shall assume that ε ≪ k, and that k depends on the length of the polymer L, R. The parameter f in these rate coefficients is the fidelity of the feedback mechanism; typically this will not be perfect, that is, f < 1 is likely in general. This model has some similarities with models we have studied in other areas of investigation; for example, the kinetics of micelle-formation in ethyl caprylate (Coveney & Wattis, 1996). In this system the breakdown of caprylate ester into monomer occurs spontaneously at some slow rate, but is massively accelerated by the presence of micelles, which have a catalytic role in the breakdown of the source species. In the present paper this mechanism is more complex since there are two monomers, and there is an additional fidelity parameter since long left-or right-handed sequences can promote the formation of the oppositely-handed (right/left) monomer as well as its own. In the caprylate system, unusual kinetic behaviour is observed when the system is initiated without any product present. The system is then effectively in a metastable state, and very little appears to happen for a long time. Following this induction time, the kinetics then proceed fairly rapidly. We expect to see similar behaviour in the current model in the case ε ≪ 1. Sandars' model (2003) differs from ours in that his imposes a maximum polymer length, typically set at five, and only polymers of this maximum length act catalytically in the breakdown of the source into monomers. We allow polymers to grow to arbitrary lengths, and chiral polymers of all lengths have some degree of efficacy in the autocatalytic feedback mechanism by which the source species decays to form chiral monomers.
The monomers will be allowed to combine to form chirally pure polymers, denoted by L n and R n according to We assume that the monomer of opposite handedness may attach to a growing polymer and so inhibit its growth; such inhibited sequences will be denoted by RL n for a polymer L n which has been terminated by an R 1 monomer, and LR n for the corresponding R n L 1 polymer. We denote the rate of such reactions by aχ. Thus we have the two rate processes However, the sequences RL n , LR n are treated as inert products which have no influence on the other rate processes; therefore their concentrations can be ignored in the mathematical modelling of the chemical reactions. Following Joyce et al. (1984) we assume no further growth of these products can occur. The system studied by Joshi et al. (2000) corresponds to 0 < χ < 1; however, we shall consider the full range of possible χ > 0.

Macroscopic model
To close the system of rate equations ensuing from the scheme (1)-(3) we shall assume that the source chemical, S, is added to the system at some constant rate S 0 . We define macroscopic quantities for the total concentration of all homochiral sequences of each type of polymer by L and R, and the mass of monomers in each set of homochiral sequences by Q, P in schema (1) by Applying the law of mass action to (1)-(3) and using the definitions of (4), we obtain the kinetic equations This constitutes a quantitative description of the polymerising system. Equations (7) and (9) hold for all integers n ≥ 2, that is for all polymer lengths, from dimers (n = 2) to infinitely long polymers. A diagrammatic summary of the rate processes involved is shown in Figure 1.  A considerable advantage of the present model is that it is possible to reduce the complexity of the corresponding infinite set of rate equations (5)-(9), resulting in a closed system of only seven equations which contains the dynamics of the full system. This is based on the quantities S, L 1 , R 1 , L, R, P, Q which evolve according to The rate processes described in terms of equations (10)-(16) are depicted in Figure 2. This reduction from an infinite system of coupled ordinary differential equations to a system of just seven ordinary differential equations is remarkable in that it is an exact simplification -no approximations have been made. Such an exact reduction would not be possible if the growth rate coefficients were dependent on polymer length. If appropriate at all it would then have to rely on approximations, and an accurate approximation may require not just the zeroth and first moments of the distributions L n , R n , but on higher moments as well. Figure 2: Illustration of our homochiral polymerisation scheme including the nonlinear feedback mechanisms. The quantities LR n , RL n are not included in the mathematical model since they play no role in the reaction scheme. The quantities L, R refer to the total numer of polymers of each handedness, and P, Q to the total mass of material in polymeric form. It is these latter quantities (P , Q) which determine the effectiveness of the catalytic breakdown of the source species into monomers.

Transforming the system of kinetic equations
Before we analyse the steady-state kinetics, it is useful to recast the system of seven ordinary differential equations in an alternative form. Instead of describing the concentration of each species separately, we assign variables for the total concentrations of source material (S), monomers (µ), sequences (N) and sequence mass (M), and a set of variables δ, θ, η for the proportions of right-handed molecules in each of µ, N, M respectively These quantities relate to the chiral purity of the system: δ is the chiral purity of the monomers (µ), θ and η describe the chiral purity of the homochiral polymer chains, θ being a numberweighted measure (corresponding to N) and η a mass-weighted measure (corresponding to M). We transform the kinetic equations (10)-(16) into the variables relating to the concentrations of polymers given by (17)-(18) and hence obtain the seven coupled equations The advantage of (19)-(25) over (10)-(16) lies in the considerable reduction in the amount of algebra required to derive solutions; for example it is easy to see that δ = θ = η = 0 satisfies (23)-(25), leaving a system of four equations (19)-(22) for the four unknowns S, µ, N, M.
The equations (19)-(25) will be used throughout the rest of the paper. In later sections we shall consider the equations for the symmetric growth of homochiral polymer sequences, (19)-(22), separately from the equations describing the chiral purity of the system, (23)-(25).

Steady-state behaviour
For the analysis of the steady-state solutions all the right-hand sides of equations (10)-(16), or equations (19)-(25) are set to zero. For the later stages of the calculation of the steady-states we ignore the small parameter ε (that is, we set ε = 0), since this will allow explicit analytical formulae to be derived and this will not greatly influence the steady-states. The reason for this is that once there are appreciable numbers of polymers of either handedness present in the system, the catalytic breakdown of source will dominate the production of monomers of both handednesses. The O(ε) terms in equations (10)-(16) will be reinstated later where the kinetics of the system starting from zero initial data are investigated.
We solve the equations (10), (12) and (15) to express the solution in terms of L 1 , R 1 .
The microscopic steady-state solution is then found to be In terms of the variables N, M, θ, η, the steady-state solution is then together with The average chain length is given by M/N, which can be written as .
The most important thing to notice about the second term here is that asymmetric (chiral, i.e. δ = 0) solutions permit longer chains to be produced. For the symmetric solution (δ = 0) we have M/N = 2 + 1/χ; whereas if δ approaches ±1 then arbitrarily large chains can be produced. We note that large inhibition rates (χ) reduce the expected sequence length (as one might expect); however, we shall show later on that this effect is compensated for since larger values for χ make it easier for the system to adopt a chiral state. It is more natural to express solutions in terms of µ, δ than L 1 , R 1 . The chiral purity of the polymers can be defined in two ways, one by the mass-weighted purity (η) and the other in terms of the number-weighted purity (θ). At steady-state we have Thus if the chiral purity of the monomer, δ, departs from zero, then the chiral purity of the homochiral sequences (θ and η) also do, and by greater amounts, as illustrated in Figure 3. This can be seen from the linearisations for small δ which give θ = 3δ and η = (5 + 2χ)δ/(1 + 2χ). For large χ the curve for η approaches that for θ. When substituted into equation (20) the above equations yield  (32), and show the accentuated effect of symmetry-breaking in the homochiral sequences compared to the monomers. The solid line corresponds to θ; the dashed lines refer to η. The steepest curve relates to χ = 0, the next to χ = 1 and that for χ = 2 is the closest to the θ curve.
Finally we have to determine the chiral purity of the monomers, δ. This is the most complicated part of the calculation so, for the sake of clarity, we now set ε = 0; δ is then given by either δ = 0, the symmetric (achiral) steady-state solution or, in terms of the fidelity, by We shall discuss this chiral solution in further detail later on (Section 3.2). One question we aim to address in the remainder of this section is how small f could be, and the system still exhibit a bifurcation to a non-symmetric state.

The symmetric, achiral, steady-state solution
It is not immediately obvious that all the solutions given above in equations (27)-(34) for the steady-state exist for all parameter values, or whether any particular solution is unique. Hence we start by considering the symmetric solution, where δ = θ = η = 0. This solution exists for all parameter values and in this case we have with L 1 = R 1 = 1 2 µ, L = R = 1 2 N, P = Q = 1 2 M. We now consider the mathematical stability of this solution, that is we aim to answer the question: 'If the system is close to the steady-state solution, will it be attracted closer to it,  (34), for fidelity in the range 0 < f < 1 and with 0 < χ < 2.
or diverge further away from it ?' To answer this question, we linearise around the solution δ = θ = η = 0. Note that the formulae for µ, N, M, S all have O(δ 2 ) correction terms, and no O(δ) terms. Thus these will be treated as constants. We only need to analyse the evolution of (θ, δ, η) over time, starting from small perturbations away from (0,0,0). Linearising equations (23)-(25) we obtain Note that we have ignored the O(ε) term. Inserting the expressions for µ, N and M from (35)-(36) we obtain which has the characteristic polynomial The Routh-Hourwitz criteria state that all solutions of λ 3 +Aλ 2 +Bλ+C = 0 satisfy ℜ(λ) < 0 if and only if A > 0, C > 0 and AB > C (for further details, see Murray, 1989). If all values of ℜλ are negative then any solution of (40) will have all the quantities η, θ and δ decaying to zero as time increases. In our example, clearly A > 0 whatever values χ and f take, but the other two conditions are less clear. The condition C > 0 gives f > f c (χ) with f c (χ) given by This value for f c agrees with f as given by (34) in the case δ = 0. The condition AB > C implies an instability when f > f c2 (χ); however, f c2 (χ) lies in the region f > 1 for all χ so this instability can be ignored, since only f ≤ 1 is physically realisable. So the symmetric solution is stable for f < f c , and unstable for f > f c . The value of f c lies between 2/9 and 4/5 and depends on χ: at χ = 0; the achiral solution is unstable for f > 4/5, whereas at large χ, the instability of the achiral solution occurs for f > 2/9.

The chiral (asymmetric) steady-state solution
Once f > f c (χ) we have the existence of two chiral steady-states as well as the achiral steadystate. Furthermore, since the Routh-Hourwitz criteria are 'necessary and sufficient' for the existence of a solution of (41) with ℜ(λ) < 0, once f > f c the achiral solution becomes unstable, so will not be observed in any physical system. Thus once the chiral solutions exist, they become the 'preferred' steady-state solutions since they are mathematical stable. Even if a system were artificially put into the symmetric (achiral) state, any small perturbation would cause some chiral imbalance and the natural kinetics of the model would then carry the whole system to one of the two steady chiral states in which one handedness of homochiral polymers dominates the other.
At f = f c a bifurcation occurs, and two mirror-image chiral solutions appear. We now briefly examine the neighbourhood of this point in more detail. These chiral steady-states are governed by equation (34); from this equation we note δ = 0 implies f = f c . Assuming small δ and expanding we find f = f c + βδ 2 with β > 0 so the bifurcation is supercritical, and the new solutions exist in f > f c where the achiral solution is unstable. Thus, as expected, the bifurcation occurs at exactly the same position as that at which the achiral steady-state becomes unstable. This is a standard result in bifurcation theory; see Guckenheimer & Holmes (1983) or Bergé et al. (1984) for more details. Although this is a curve in (f, χ) space, we shall consider f to be the primary bifurcation parameter, with the bifurcation point f c depending on χ.
In figure 4 we show the chiral steady-state solutions given by (34) parametrically, plotting δ against f and χ. The bifurcation occurs at f c = 4/5 if χ = 0, reducing to f c = 2/9 in the limit of large χ. Thus we see the beneficial effect of sequence inhibition if one wants a system which undergoes a symmetry-breaking bifurcation at small values of the fidelity parameter f . Figure 5 shows the chiral purities of the polymers, firstly weighted by number (θ) and then by mass (η). Both show that the chiral purity of the polymers is much greater than that of the monomer (compare Figure 5 with Figure 4). Note also that for the range of fidelities 2/9 < f < 4/5 the bifurcation to asymmetric (chiral) solutions can occur by increasing χ.
Another natural question to pose here is whether only small and moderate values of δ can be accessed; we put ν = 1 − δ with ν ≪ 1 to determine under what conditions δ can approach unity. We find from equation (34) that f = 1 − χν. In this case we also expect θ and η to be near unity. In fact the chiral purity of the polymeric sequences is much enhanced over the  (34). This is to be compared with Figure 4, illustrating that the polymer concentration and mass show more extreme chiral purities than the monomers (described by δ). chiral purity of monomers, since a two-term expansion of equation (32) gives for ν ≪ 1 (in the special case χ = 0 we have η ∼ 1 − 1 16 ν 5 ). The concentrations scale with ν according to together with S = O(ν 3/2 ). In summary we find that chiral solutions do indeed exist provided f > f c and this is easier to satisfy at larger inhibition rates (χ). Thus the presence of stronger cross inhibition aids the manifestation of chiral steady-states. Also, the higher the fidelity (f ), the greater the dominance of one chirality over the other.
We now turn to an analysis of the kinetics of sequence growth, in order to determine at what stage of the reaction the system is likely to manifest a chiral state.

Kinetics and stability of achiral solutions
In this section we reintroduce the O(ε) term omitted in the analysis following (32) and solve the kinetic problem for the achiral solution, in the limit ε ≪ 1. We assume that initially there are no monomers (µ = 0), no polymers (M = N = 0) and no precursor (S = 0), though this source material is added continuously to the system starting at time t = 0. If ε = 0 then since there are no polymers, the precursor cannot be broken down and no polymers will ever form, so we need the O(ε) term to produce some monomer. In section 4.2 we analyse the stability of the solution, that is, whether small, random perturbations to such a solution grow or are damped out, as the total number of polymers and monomers grow from zero concentrations.

Growth of the achiral solution
We assume a set of initial conditions in which η = θ = 0 = δ, and then η = θ = 0 = δ for all subsequent times. This reduces the system (19)-(25) to a system of three ordinary differential equations which we solve asymptotically in the limit ε ≪ 1. Our aim is to develop matched asymptotic expansions for the solution of through a series of timescales.

Timescale
For this timescale the appropriate scalings are t = ε −1/5 t 1 together with thus this is a long induction time over which the leading order equations are where prime denotes d/dt 1 . This system has the solution where µ 1 is given by the solution of which unfortunately is not explicitly available. However, we can see that the timescale ends abruptly with µ 1 , N 1 , M 1 all diverging as t 1 → t 1c , according to for some constant t 1c . These relationships help us determine the scalings relevant in the next timescale. In this timescale we have seen the accumulation of source material, but this is only slowly converted into monomers and chains, so both of these grow very slowly, causing a big build up of source material until, at the end of this timescale, we see the concentration of chains increase to the level where the catalytic mechanism becomes active and accelerates the formation of monomers and chains.

Timescale II
In this timescale all quantities are large and evolve quickly. To be specific we have together with t = ε −1/5 t 1c + ε 1/5 t 2 . Using primes to denote d/dt 2 , the leading order equations are As the chains are present in large enough quantities for the catalytic mechanism to be active, and since there is a large amount of source material present at the start of this timescale, this source material is rapidly converted into monomers and chains so that µ 2 , N 2 and M 2 all increase at the expense of the source species, S, whose concentration now monotonically reduces, so that the only significant simplification in the equations is in the equation for S 2 .
No explicit solution which matches back into Timescale I is available; however, the form of the large-time solution in Timescale II can be determined. Consider new timescale given by 1 µ 2 d dt 2 = d dτ , hence t 2 = 2 a dτ µ 2 (τ ) , then the system (56) can be written as which is a linear system, together with the equation from which Q = 2kS 2 M 2 /aµ 2 is obtained. The solution of equations (57)-(59) is given by a combination of a complimentary function which solves the system with Q(t) = 0 and a particular solution which satisfies the Q(t) input term in equation (57). For general parameter values, the forms of these solutions cannot be explicitly determined, however we shall study the two particular limits of large and small χ in which approximations can be obtained. These correspond respectively to the cases where the solution is dominated by the particular solution and the complimentary function at large times.

The solution in timescale II for χ ≫ 1
If we assume that the input function Q has the form Q 0 e −λτ , then the particular solution has the form In the calculation of Q, we then have S 2 ∼ e −λτ too, and this assumption implies λ = 2k M 2 /a µ 2 , and Q ∼ λ S 2 e −λτ . In the calculation of the prefactors µ 2 , N 2 , M 2 we then have to solve a cubic, and there are constraints that all the prefactors must be positive.
Inserting these into the equations (57)-(59) we obtain and the cubic equation for λ is For asymptotically large χ, this has two roots near λ = χ, and one near λ = 0. The larger roots are at λ = χ ± 2k/aχ, one of which violates the condition λ < χ and the other leads to a solution which rapidly decays in time. The physically relevant solution corresponds to λ ∼ 4k/aχ, hence the solution Since this result has been derived on the basis of large χ, we see a slow exponential decay in τ . The complimentary function also decays exponentially in τ , but with exponents λ 1 = 3, λ 2 = χ − 1 and λ 3 = χ, hence the complimentary function decays much more rapidly than the particular solution.
In terms of the original timescale t 2 the solution (64) leads to This agrees with the numerical observed results, which suggest that for χ ≫ 1 we observe that all quantities decay with 1/t 2 . Figure 6 shows a numerical solution of the system produced by Matlab 6.5.0 Release 13; the right-hand graph shows in detail the decay of the concentrations at the end of the second timescale. When the concentrations S, µ, N, M become O(1) then other terms become relevant in the kinetic equations, and a further timescale is required to describe the final approach to equilibrium.

4.1.4
Timescale III for larger χ: t = ε −1/5 t 1c + O (1) As S 2 , µ 2 , N 2 and M 2 all decay like 1/t 2 at the end of the previous timescale, the new first term to become significant is the S 0 input into the S equation. This becomes significant when t 2 = ε −1/5 , thus the third timescale is t = ε −1/5 t 1c + t 3 . In this timescale all of S, µ, M and N are O(1). Thus the leading order equations are and over this timescale the system approaches its steady-state. Since at leading order the εS terms are neglected, the leading order steady-state solution approached is precisely that described in Section 3. Overall, we see that there is a long induction time, of O(ε −1/5 ), followed by a some rapid kinetics during which the system explores states a long way from its steady-state, and then over a relatively fast (O(1)) timescale the system approaches its steady-state. We now derive the kinetics of the achiral solution for smaller χ before going to address the stability of the growing achiral solution.

The solution in timescale II for χ ≪ 1
Numerical simulations suggest that at the end of TII, S 2 and µ 2 decay exponentially (in t 2 ) to zero, with M 2 and N 2 tending to constant values; eventually reaching their O(1) steady-state values over subsequent and much longer timescales (see later subsections and Figure 7 for example).
The linear system (57)-(59) is solved by the sum of a particular solution and a complimentary function. For smaller χ the large-time solution is dominated by the complimentary solution component of the solution (i.e. that with Q = 0) which is given by a combination of exponentials. In the case χ ≪ 1, the eigenvalues are λ 1 = −χ, λ 2,3 = −1 ± i √ 2χ, thus applying the special result for repeated roots, at leading order, we have the complimentary function The solution in this case corresponds to the scenario in which µ 2 → 0 as τ → τ c < ∞. Let us assume then t 2 = 2 a dτ µ 2 (τ ) implies τ c − τ = e −aµ 2 t 2 /2 so that τ → τ c corresponds to t 2 → ∞. We thus have µ 2 ∼ µ 2 e −aµ 2 t 2 /2 . The solution of equation (60) yields If 1 2 aµ 2 > kM 2 , then this solution predicts that µ 2 decays faster than S 2 ; however, towards the end of the second timescale the dominant terms in the equation for µ 2 are µ ′ 2 = kS 2 M 2 − 1 2 aµ 2 N 2 . This equation implies that µ 2 cannot decay faster than S 2 , and thus we expect 1 2 aµ 2 ≤ kM 2 . If equality holds, then both µ and S reach O(ε 1/5 ) at the same time. In summary, towards the end of this timescale, we find µ 2 and S 2 decaying exponentially in t 2 with S 2 ∼ S 2 e −kM 2 t 2 , µ 2 ∼ µ 2 e −aµ 2 t 2 /2 , together with N 2 → N 2 and M 2 → M 2 . This agrees with our observations of the numerical solution discussed in the opening paragraph of this subsection. The system then passes straight into Timescale III, which becomes relevant when kSM ∼ S 0 . This occurs when S = O(ε 1/5 ) which happens when t 2 ∼ t 2c = (2/(5kM 2 )) log(1/ε). However, if kM 2 > 1 2 aµ 2 then S 2 decays at a faster rate than µ 2 and there is an further timescale between this and timescale III (below) over which S saturates while µ decreases further. This timescale is given by for some −1/5 < γ < 1/5 and t = ε −1/5 t 1c + ε 1/5 log(1/ε)t 2c + ε 1/5 t. The governing equations are then where prime denotes d/d t. Over this timescale, M and N remain constant, S equilibrates to S 0 /k M and µ continues to decrease exponentially, with rate − 1 2 a N. When µ reaches O(ε 1/5 ), timescale III is entered; this occurs when t = O(log(1/ε)), so the relationship between t 3 and t in Timescale III given below remains valid following a redefinition of t 2c to incorporate this extra shift in time.
4.1.6 Timescale III for smaller χ: t = ε −1/5 t 1c + ε 1/5 log(1/ε)t 2c + O(ε 1/5 ) In the new timescale we have together with t = t 1c ε −1/5 + t 2c ε 1/5 log(1/ε) + ε 1/5 t 3 . Using prime to denote time derivative with respect to the new time variable, t 3 , the leading order equations are Over this rapid timescale, M and N do not change from their values at the end of the TII whilst µ and S equilibrate to the values The evolution of M and N occurs over a longer timescale with µ and S constrained to their respective local equilibrium values.
4.1.7 Timescale IV for smaller χ: t = ε −1/5 t 1c + O(ε −1/5 ) Since all concentrations approached constants at the end of TIII, their magnitudes remain unchanged for the fourth timescale, only the scaling for t changes, we now write with t = t 1c ε −1/5 + t 2c ε 1/5 log(1/ε) + ε −1/5 t 4 , and so obtain the equations These imply kS 4 M 4 = S 0 and µ 4 N 4 = 2S 0 /a(1 + χ), thus we have the solution for some constant C. This timescale ends with µ 4 and S 4 increasing hence becoming larger in the next timescale, and M 4 , N 4 decaying hence being smaller in the next timescale.

Timescale V for smaller
The final timescale is given by all of S, µ, N and M being O(1) and varying on an O(1) timescale, thus corresponds to the approach to the global equilibrium solution given by (35)-(36).

Summary
For all values of the parameters, the solution starts with a long induction period during which the concentration of the precursor species S becomes large. There follows a short period of very rapid kinetics where all concentrations become large, and then the precursor and monomer concentrations decay. For larger values of χ the concentration of polymers also decays and the steady-state is reached relatively rapidly. For smaller χ the monomer concentration and that of the source species become very small and the polymer concentrations (mass-weighted and number weighted) both remain high and slowly evolve to their steady-state values over a longer timescale, which is of similar length to the induction timescale. These two distinct behaviours are illustrated in Figures 6 and 7.

Stability of the achiral growing solution
Having determined the form of the kinetic behaviour for the achiral solution (η = θ = δ = 0), we now consider the linear stability of this solution. Assuming there is some small random perturbation during the evolution, we ask whether a perturbation grows or decays as time progresses. We shall use the already determined solution for S(t), µ(t) N(t) and M(t) (from section 4.1), and assume that the perturbation does not make any alteration to these total concentrations of source, monomer and polymer at leading order. This assumption was seen to be valid in the case of the steady-state solution, and since the kinetic equations are symmetric under the transformation (δ, θ, η) → (−δ, −θ, −η), we expect modifications to µ, N, M, S also to be of second order in δ, θ, η. From equations (23)-(25) we have

Stability of achiral solution in Timescale I
We focus our attention on the linear stability of δ(t), θ(t) and η(t) as given by (37)-(39). In the first timescale, using the scalings (50), we find the simplified linear stability problem It is clear that this matrix has a simpler structure than the general case, which has only one zero entry, as given in equation (40). The eigenvalues of the matrix above are all negative if f < 1/2, but if f > 1/2 then one is positive, indicating that a perturbation away from θ, η, δ = 0 would increase in size as time progresses. The temporal evolution of such perturbations is non-trivial since M 1 , N 1 , µ 1 and S 1 are all time-dependent.

Stability of achiral solution in Timescale II, smaller χ
At the end of Timescale II if χ is small then the matrix in equation (83) has the form All elements in the bottom row of this matrix are divergent since µ 2 , S 2 → 0 as t 2 → ∞. When t 2 becomes large the eigenvalues of this matrix are given by solutions of the cubic Applying the Routh-Hourwitz criteria for the signs of the real parts of the eigenvalues, we find stability of the achiral solution requires A > 0 (which always holds), and C > 0 which fails when f > aN 2 µ 2 /2kM 2 S 2 ; and a similar but more stringent inequality from the condition AB > C. So if S 2 decays faster than µ 2 then the symmetric solution is stable, and if S 2 and µ 2 decay at the same rate the stability depends on N 2 and M 2 and requires f > 1/(2(1 − kM 2 /aN 2 )).

Stability of achiral solution in Timescale III, smaller χ
With the scalings of (76) the matrix (83) takes the form (89) only the leading order entries for each term, the cubic governing stability can be written Rewriting this as λ 3 + Aλ 2 + Bλ + C, the Routh-Hourwitz criterion A > 0 is always satisfied. The condition C > 0 implies f < aµ 3 N 3 /2kS 3 M 3 which, at large times, reduces to f < 1, and as such is met in the large-time limit. Finally, AB > C fails if f > N 3 /(N 3 + χM 3 ) indicating an instability; whilst this depends on the unknown number-weighted and massweighted polymer concentrations, an instability is certainly possible for large enough values of the fidelity parameter, f .

Stability of achiral solution in Timescale IV, smaller χ
Since the scalings for the concentrations S, µ, N and M are as in timescale III, the only change from (89) As the solution progresses through timescale IV, the criterion for an instability to exist, namely f > N 4 /(N 4 + χM 4 ), becomes easier to satisfy, since N 4 decays slightly faster than M 4 , as shown by equations (81)-(82).

Summary
In this section we have analysed the kinetics of the concentrations of monomer, source and polymer as they evolve from zero to steady-state following a symmetric (achiral) solution.
For larger values of χ we have found three regimes through which the system evolves. Firstly there is a long induction period (of O(ε −1/5 )) during which the source material builds up (until S becomes O(ε −1/5 )); during this time the monomer concentrations remain small (O(ε 3/5 )) and the concentrations of polymer are extremely small (O(ε)). This timescale ends abruptly as the catalytic feedback of polymer accelerates the breakdown of S into monomer. In the second timescale, which is very brief (O(ε 1/5 )), all monomer and polymer concentrations become large (O(ε −1/5 )) and decay towards steady-state. Finally, over the third timescale all concentrations converge to their steady-state values.
In section 4.2 we analysed the linear stability of the growing achiral state through the sequence of timescales. By linear stability we mean that small external random forces which cause a chiral imbalance will be damped and reduce in amplitude. An instability indicates that such a perturbation will grow and so the system will undergo a symmetry-breaking bifurcation to occur during the evolution.
As the system approaches steady-state in the third timescale we expect to regain the stability criteria f > f c (χ) with f c given by equation (42). However, this criterion for the instability of the achiral solution is not valid in all the timescales; in the first timescale we find the alternative criteria of f > 1 2 for symmetry-breaking to occur. Note that this is independent of χ, since in the first timescale the polymer and monomer concentrations are so low that inhibition of a homochiral polymer by a monomer of the opposite chirality is negligible; this leads to a considerable simplification of the linear stability analysis during the first timescale. In the second timescale we find an instability for f < 2/3.
For f < 2/9 the achiral solution is always stable to such perturbations; whereas for 2/9 < f < 1/2 the system is unstable to such perturbations only in the final stage of the kinetics, when the equilibrium solution is being approached and even then only for some values of χ. For 1/2 < f < 2/3 the kinetics are more complicated since the system is unstable during the first timescale, but linearly stable during the second timescale and unstable during the final approach to equilibrium. Thus for these f -values the chiral purity may oscillate, but the system will eventually approach a chiral state.
For smaller values of χ the kinetics are more complex: there is still a long induction time followed by a period of rapid kinetics. This is more complicated, being split into various timescales; however the end result is always low concentrations of monomer and source species and large concentrations of polymer. There follows a long timescale over which the polymer concentrations reduce towards steady-state and the monomer concentrations and source species increase to steady-state.
For smaller values of χ an instability occurs in Timescale I for f > 1 2 , the instability persists in timescale II, now depending on f > 1/(2(1 − kM 2 /aN 2 )). In timescale III and IV an instability requires f > N/(N + χM), and at steady-state f > 4/5 − 34χ/25. Thus, for smaller χ there is still the possibility of a symmetry-breaking bifurcation occurring, though it requires a larger fidelity parameter.
Although there advantages in χ being large if one is seeking a symmetry-breaking bifurcation, it should be noted that smaller χ has other advantages, in that it allows polymers to form in larger concentrations, and these persist for longer times.

Perfect fidelity
In the extreme case f = 1 the feedback mechanism breaks down the precursor species (S) into chirally pure monomers with unit probability. Instead of a chirally pure homochiral steady state, there is a bifurcation to a state in which δ asymptotically approaches ±1. We refer to such a state as a fully-bifurcated state. From equation (34) we see that such a state cannot arise if f < 1. The large-time asymptotics of the fully-bifurcated state differ significantly from f < 1 since now there is no steady-state solution; instead we find unlimited growth of one set of homochiral polymer sequences and decay to zero for the sequences of opposite handedness. In this case (still ignoring ε) the large time asymptotics are given by Assuming δ, θ and η all decay to zero in the large time limit, we find the leading order equations S 0 = k S M , M = a µ N, These imply Thus in terms of the concentrations of each chirality we have together with S ∼ 1/kt, and for some constants L 1 , L, Q. So in this case no finite steady-state solution is approached. As one might expect from the asymptotic expansions (44)-(45) in the case δ → ±1 we observe the unbounded growth of one type of homochiral sequence and, specifically, unbounded growth in the number of chains (N), the mass of material in polymeric form (M) and the average length (M/N ∼ (aS 0 t 2 ) 1/3 ). Concentrations of the sequences of opposite homochirality decay rapidly, and we expect that the average chain length approaches two, implying Q = 2 L.
All the above analysis has been for the simplified case for which ε = 0; however, with one chain type decaying to arbitrarily small concentrations we may expect that the O(ε) term in (11) is no longer negligible in this limit. Retaining the O(ε) term in this equation yields a slightly different scaling in the large time asymptotics for the high-fidelity case f = 1. We now have together with S 0 ∼ 1/kt and L 1 ∼ ε kχ(3a 2 S 2 0 t 4 ) 1/3 , L ∼ ε 2 k 2 χ 3 S 0 (3a 2 S 2 0 ) 1/3 t 3 , Q ∼ 2ε 2 k 2 χ 2 (3a 2 S 2 0 ) 1/3 t 3 .
Thus, we see the convergence to full chiral purity is more rapid for polymers than for monomers. In all the above analysis there is another chiral solution in which the left-handed homochiral polymer sequences are dominant, and the right-handed homochiral sequences have concentrations which decay to zero asymptotically.

Discussion
After introducing our model in Section 2, we analysed its steady-states, and found that the symmetric steady-state solution exists for all parameter values, but that there are other solutions when the relative inhibition rate χ and the fidelity f are large enough. The critical combination is f > f c = (4 + 2χ)(1 + 2χ) (5 + 6χ)(1 + 3χ) .
At the point f = f c there is a supercritical pitchfork bifurcation, where two unstable steadystate solutions connect to the solution δ = 0 and make δ = 0 unstable for f > f c . When this inequality is satisfied there are two stable steady-state solutions with δ non-zero, that is there are chiral solutions as well as the achiral solution, and the achiral solution is unstable, so any physical system will generically be attracted to one or other of the chiral solutions. An important effect to note from this formula is the role that cross-inhibition plays in making the asymmetric solutions accessible at low values of the fidelity. For small χ, the fidelity has to exceed f c = 0.8 in order to obtain a symmetry-breaking solution; whereas at large χ, this bifurcation point reduces to f c = 0.22 -a dramatic reduction.
In Section 4 we analysed the kinetics of chain growth in a symmetric system, and found that there is a long induction time, during which a large stock of precursor chemical accumulates; an approximate, linear stability calculation shows that during this time, the achiral solution is unstable if f > 1/2. This behaviour is followed by a short timescale over which the precursor species is converted to monomers which are then polymerised. For this short time, monomers and chains are present in large concentrations. The concentrations of chains, monomers and precursor then all decay to their steady-state values, which, if the parameters f, χ are such that an asymmetric steady-state exists, and the system has experienced some external perturbation away from the symmetric state, will be the chiral steady-state discussed earlier (section 3.2). In such a state both monomers and chains have a net chirality or handedness. Even for quite modest values of the chiral purity of monomer (say δ = 0.7), the chiral purity of chains is extremely close to unity (θ = 0.990, η = 0.995 at χ = 2); see Figure 3, and compare Figures 4 and 5.
Finally we have described the large-time asymptotics of the 'fully' bifurcated case which arises when f = 1, wherein chiral purities (δ, θ and η) approach unity in the large time limit; this remains true even when we reintroduce the term which describes the slow spontaneous achiral decay of precursor species into both enantiomeric forms of monomer.

Conclusions
In previous work we showed how qualitatively similar instabilities can lead to the massive amplification of self-replicating RNA polymer sequences over less efficient replicators (Wattis & Coveney, 1999).
In the present paper we have shown that an initially achiral system capable of stepwise polymerisation to homochiral polymer sequences with inhibition from the opposite-handed monomer, is subject to strong instabilities that drive the system overwhelmingly to one or other handedness for all homochiral sequences present. Mechanisms within this class may have played a rôle during the early stages of molecular evolution in determining the chirality of biologically relevant macromolecules, such as nucleotides and proteins, and there is experimental evidence of this behaviour in the literature, for examples see Hitz et al. (2001Hitz et al. ( , 2002Hitz et al. ( , 2003 and Joshi (2000). Although in the system studied by Joshi, addition of the correct enantiomer to a growing polymer chain is more favourable than the wrong one, we have shown that if this cross-inhibition is stronger, then the system is more likely to undergo a symmetry-breaking bifurcation.
Studies of this kind confirm the scope and power of modern methods of theoretical analysis for nonlinear dynamical systems of the kind that abound along the pathway towards the origins of life (Coveney, 1994).