Spontaneous chiral symmetry breaking in early molecular networks

Background An important facet of early biological evolution is the selection of chiral enantiomers for molecules such as amino acids and sugars. The origin of this symmetry breaking is a long-standing question in molecular evolution. Previous models addressing this question include particular kinetic properties such as autocatalysis or negative cross catalysis. Results We propose here a more general kinetic formalism for early enantioselection, based on our previously described Graded Autocatalysis Replication Domain (GARD) model for prebiotic evolution in molecular assemblies. This model is adapted here to the case of chiral molecules by applying symmetry constraints to mutual molecular recognition within the assembly. The ensuing dynamics shows spontaneous chiral symmetry breaking, with transitions towards stationary compositional states (composomes) enriched with one of the two enantiomers for some of the constituent molecule types. Furthermore, one or the other of the two antipodal compositional states of the assembly also shows time-dependent selection. Conclusion It follows that chiral selection may be an emergent consequence of early catalytic molecular networks rather than a prerequisite for the initiation of primeval life processes. Elaborations of this model could help explain the prevalent chiral homogeneity in present-day living cells. Reviewers This article was reviewed by Boris Rubinstein (nominated by Arcady Mushegian), Arcady Mushegian, Meir Lahav (nominated by Yitzhak Pilpel) and Sergei Maslov.


Background
The derivation of chemical reactions that spontaneously generate an excess of one enantiomeric form (i.e. one of two stereo-isomers of an asymmetric molecule endowed with the property of handedness or chirality, and mutually related by mirror symmetry) has been a central ambition of numerous theoretical and experimental studies [1][2][3][4][5][6]. The challenge is to depart from a racemic mixtures (having equal amounts of both isomers), and reach enantiomeric excess without the aid of external chiral selectors. Thus (reviewed in [5]), some authors have proposed that a catastrophic symmetry breaking event was necessary to explain why in a class of biomolecules (e.g. amino acids) all members have the same chiral configuration. Energy imbalance of enantiomers due to a lack of antimatter parity, or enantioselective breakdown by circularly polarized light from space was invoked. It was argued, however, that a viable statistical model could replace these cosmic explanations, a model merely based on evolutionary properties such as propagation and competition.
Indeed, several studies invoked relatively simple kinetic models, in which initial racemates with fluctuations undergo reactions that lead to chiral purity, thus demonstrating the plausibility of symmetry breaking in a nonequilibrium regimen [5,[7][8][9]. Many such treatises assume that chiral selection has occurred under abiotic conditions, and preceded (or even served as a prerequisite for) life's origin. Among these are models that involve bifurcation in small molecules [8][9][10] In parallel, systems were reported that involve polymerization [11,12] as well as interactions within crystals (reviewed in [12]). The basic principles that guide such papers include the notion of statistical fluctuations, namely that in an ensemble of asymmetric molecules of a given type, there will always be an excess of one enantiomer, particularly apparent in small ensembles, and that such fortuitous excess may be greatly amplified by catalytic or replicative reactions [8]. The present paper rests on such view, and attempt to provide a novel concrete and quantitative framework for its realization.
Life is believed to have emerged by self organization processes occurring within a random and highly heterogeneous chemical environment [10,13]. One of the hallmarks of some other prebiotic evolution studies is the assumption that homochirality (the prevalence of only one of the two chiral isomers), widespread in present-day life, has emerged as part of the processes that led to cellular life [14][15][16][17][18][19]. For example, it has been argued [14] that information theory conclusions can explain why chiral building blocks, as well as sets thereof, are necessary in living systems, and that simplest forms of life likely constituted autocatalytic reactions such as the Soai reaction [4] where a chiral product acts as a chiral catalyst for its own production.
It is thus crucial to ask how chiral symmetry breaking could become possible under the conditions that prevailed at the early emergence life (see for example [19], [20] and references thereof ). By one school of thought, the origin of life is proposed to have occurred through kinetically self organizing processes controlled by defined chemical interaction networks [21][22][23][24][25][26]. In this respect, models accounting for life's origin could be helpful for the understanding the generation of chiral purity. A case in point is the Graded Autocatalysis Replication Domain (GARD) model we have developed [23,[27][28][29][30][31][32][33][34]. The model entails a chemically diverse set of mutually catalytic amphiphiles that spontaneously aggregate to form molecular assemblies (e.g. micelles, a "Lipid World" scenario [29], see also [35]). It was shown that such assemblies often self-organize into kinetically stable mutually catalytic networks, termed composomes, display homeostatic growth and reproduction-like processes, as well as a limited capacity for selection [30,[36][37][38]. One unique aspect of GARD, with respect to other proposed models, is that the life-like chemical networks are not a-priori engineered for this specific purpose, but rather spontaneously emerge.
GARD nominally belongs to a set of models that assume that early prebiotic evolution took place completely in aqueous solution, in distinction from scenarios that point to a liquid-solid interface as indispensable for life's origin. Mineral interfaces are presumed to have provided catalysis, compartmentalization and sometimes also a free energy source ( [39][40][41] and refs thereof ). However, all three such aspects are provided by in GARD and similar models via lipid-water interfaces. Furthermore, lipid-based models are not contradictory to the involve-ment of solid interfaces, which lipid assemblies could interact.
We ask here whether a collection of simple organic chiral amphiphilic molecules could assemble into a reaction network that, in the absence of any external chiral perturbation, would catalyze preferential enantiomer selection. In contrast to previous theoretical models, we do not engineer a reaction network that is designed to generate symmetry breaking, but rather apply random chemistry by generating a large diversity of molecules and chemical reactions. Our results suggest that random catalytic networks may self organize into chirally selected compositions, which in addition portray homeostatic growth and dynamic properties akin to self reproduction.

Chiral compositional dynamics shows enantiomeric selection
To investigate the network behavior of collection of chiral amphiphilic molecules, we employ the GARD model (Methods) in ways that accommodate chirality (Chiral-GARD or C-GARD). We thus consider an environment with a population of N G types of asymmetric molecules in a racemic mixture that contains equal amounts of the D and L optical isomers of each molecule type. For sufficiently complex molecular structures it is justified to assume that essentially all molecules are chiral [5,[51][52][53][54][55] ( Figure 1). All 2 × N G molecule types are treated as different compounds with different kinetic parameters, keeping in mind that they actually constitute 100 enantiomer pairs (Methods and Figure 2).
We asked whether in C-GARD, composomes may display enantioselection, in analogy to the chemical selection seen in the dynamics of the GARD model [28][29][30]32,34]. For this, we employed a measure termed here "weak enantioselection", denoted W W (Methods, Eq. 7). The top panel of Figure 3 shows a correlation diagram for all 4000 time points in one such simulation. Importantly, appreciable non zero values of W W occur at most time points, suggesting that many of the compositions are symmetry-broken. We note that at high W W values, each of the compounds has a high enantiomeric excess, although not necessarily in the same handedness for all compounds ("Strong enantioselection" W S ≈0, Eq. 8). The simulation further displays the appearance of distinct composomes along the time axis, and it is apparent that different composomes have different average W W values. Figure 4 shows a typical variation of W W as a function of time for several values of σ ε , a parameter determining the distribution of the enantiodiscrimination α (Eqs. 4 and 6). Larger values of σ ε result in increasing W W to the point of reaching nearly complete weak enantioselection (W W T1), as also seen at the top apex of Figure 5. The observed W W fluctuations at low σ ε result from the assembly size variation with assembly growth and split. Such fluctuations are significantly reduced at higher σ ε , indicating compositional stability through the action of the catalytic network [21][22][23]29,30,56]. The variation on longer time scale seen for σ ε = 3 represent transitions among different composomes with varying degrees of chiral inhomogeneity.
Interestingly, significant non-zero values of W W are obtained even in the absence of any mutually catalytic effect, i.e. when both enantiodiscrimination and cata- Red represents a case where one carbon is replaced by a hetero atom, and blue denotes a case of no hetero atom. Data is taken from [54]. This figure demonstrates that for sufficiently complex molecular structures it is a good approximation to assume that all molecules are chiral [51][52][53][54][55]. A B C lytic-potency-related parameter are low (σ ε and σ, respectively. See Eqs. 6 and 3) ( Figure 5). This enantiomer excess is distinct from the presently described dynamic enantioselection, and relates to previously published predictions regarding statistical fluctuations at low molecular copy numbers [5,[57][58][59]. Intriguingly, allowing high values of the rate enhancements suppress this statistical effect due to compositional bias whereby only a few molecules are present at low copy number (lower left apex of Figure 5). Figure 6 shows a global analysis of the dependence of enantioselection on molecular enantiodiscrimination, integrating the results of 6000 different simulations. A probability distribution for the average W W values of assemblies is plotted for three different σ ε values. Increasing σ ε enhances the probability of assemblies to have high W W , yet even for the highest σ ε studied here (σ ε = 6) there is an almost even chance for assemblies to show high or low W W . This may indicates a stochastic effect: high enantiodiscrimination is necessary, but not sufficient to lead to symmetry breaking.

Antipodal composomes
We analyzed the symmetry properties of composomes that emerged in assemblies of chiral compounds. Figure 7 shows a correlation diagram in which different types of It is seen that the typical W W increases with σ ε , because higher enantiomeric discrimination (α) values are allowed. The saw-tooth patter arises from growth-fission cycles of the C-GARD assembly.   composomes are apparent, with symmetry properties highlighted by specific colors (inset). It may be seen that at certain time intervals highly asymmetric composomes appear while at different times the composomes are not symmetry-broken (C1). The figure further demonstrates the existence of antipodal composomes (C2 and C3, see compositional bar charts in Figure 8), each with broken symmetry, and showing mutual compositional mirror relationships as indicated by the off-diagonal blue areas (See Methods, Eq. 5). Interestingly, in certain cases network dynamics leads to abrupt transitions among composomes with different symmetry properties, including between antipodal composomes (e.g. near time point 2400, where a transition occurs directly between composomes C2 and C3).

The dynamics of composome populations
We examined the capacity of the C-GARD model to portray evolution-like processes with selection of particular enantioselected molecular assemblies. To this end, we explored the competitive coexistence of numerous assemblies with different degrees of chiral symmetry. We simulate the time course of a single specific assembly and reconstructing from it an approximation for a time dependent population behavior, similar to a previously reported procedure [30]. The approximated population size of a given composome C m was computed as , where t m is a cumulative elapsed time period of assembly homeostatic growth/fission while being in the composomal state C m and ? m is the average time between fission events characterizing that composome. The rational for such a computation is that if a bona fide population of C-GARD assemblies were observed for a time period , the assemblies in composomal state C m will have undergone ~ fission events.

Figures 9 and 10A
shows results for a particular simulation. A competition-like behavior between the composomes is observed in which the relative proportions composomal states in the emulated population changes over the time course. This reflects the differences in composome "fecundity", represented by their time to fission τ m , manifested by the exponential nature of the population growth. We note that these simulations do not include events of assembly "death", reflecting an assumption that the simulated population size is sufficiently small to justify the buffered environment assumption [27,30].
Among all the composomes in Figure 9, only composomes 3 and 7 are appreciably enantioselected and the  rest are racemic ( Figure 10B). Other simulations show a majority of symmetry broken composomes (not shown). Interestingly, while composome 7 (red solid line) and composome 3 (red dashed line) are nearly antipodal, composome 7 which fortuitously first emerges later than composome 3, ends up with a simulated population size more than 100 fold larger towards the end of the examined time period. This is rationalized by the notion that each composomal state constitutes an ensemble of disparate, though similar compositions and small chance fluctuations may lead to noticeable differences in long-term dynamic behavior (cf. [60]).

Multi-component kinetic enantioselection
The breaking of chiral symmetry in multi-molecular assemblies presented here constitutes a distinct class of plausible stereo-selective processes. An advantage over other models is by considering many pairs of enantiomers, thus offering systems-related enantioselective mechanisms. The same concept is manifested in a reported kinetic simulation of a simple network of replicating peptides [61], as well as in mutual interaction within molecular assemblies such as monolayers and 3-dimensional crystals ( [12] and references thereof). It should be also noted that heterogeneous multi-component chemistry is more appropriate for describing early symmetry breaking processes, since prebiotic environments have likely been highly chemically diverse [27,33,62]. The C-GARD model presented here specifically assumes that symmetry breaking has occurred within assemblies of amphiphilic, lipid-like molecules. This, and similar concepts involving lipid molecular assemblies [63], stochastic aggregates [64] and crystalline arrays [12,16,45,50] have been previously explored by others. Enantioselection is often portrayed as a non equilibrium kinetic process [3,9,61,65]. Many relevant kinetic formalisms [3,[66][67][68] are based on the original model of Frank for spontaneous mirror symmetry breaking [9], which assumes a "chemical substance which is a catalyst for its own production and an anti-catalyst for the production of its optical antimer". Another set of models derives from the mechanism proposed by Kondepudi [8] with two autocatalytic achiral precursors whose dynamics result with homochirality. The C-GARD model presented here is also based on defining a set of kinetic equations for the different reaction paths. However, these kinetic equations are not designed a-priori to produce symmetry breaking. Rather, enantioselection spontaneously emerges in some of the molecular assemblies and is propagated by compositional homeostasis. The detailed mechanisms responsible for enantioselection may vary from one assembly to another and could conceivably include autocatalysis and mutual pairwise catalysis. Another advantage of our model is its ability to provide an estimate of the propensity of assemblies with different levels of symmetry breaking based on the kinetic equations derived from a statistical formalism based on molecular interaction [32,[69][70][71].
We assume in the present analysis that all molecules in the C-GARD simulation are asymmetric (Figure 1). An interesting question is what would be the dynamic fate of  Figure 9, with color scheme as in Figure 7; B, the compositional bar charts for all 9 composomes.

A B
symmetric compounds intermixed with a majority of asymmetric ones. Consider Eq. 4 for the case shown in Figure 2, for i = 17 and j = 11. Making molecule 11 symmetric, by the loss of a chiral center, which is assumed to have a minimal chemical effect otherwise, necessitates , where X denotes the symmetrized molecule. It may be rather safely assumed that (on average of many such cases) the value of will be somewhere between and , perhaps their geometric mean. Thus, symmetric molecules will not have an appreciable kinetic advantage or disadvantage.

Symmetry breaking due to statistical constraints
An interesting aspect of the C-GARD model is a capacity to show a distinction between kinetically-controlled chiral selection, and apparent enantioselection in the absence of stereospecific molecular recognition. The latter arises due to statistical fluctuations relating to assembly size and chemical heterogeneity, as has been explored previously for polymer systems [5,[57][58][59] and for generally diverse chemical systems [5]. In contrast, the chiral constitution of assemblies of kinetically-interacting molecules displays fluctuations between high and low values of W W , in agreement with the general characteristics predicted for non equilibrium symmetry breaking systems [61,65,72,73]. A mechanism for the symmetry breaking in non equilibrium systems was previously proposed [9] and additionally revised [9,65,73]. There is a relationship between this mechanism and the one depicted by C-GARD. Basically Frank's model constitutes a special case of a two dimensional C-GARD, in which the autocatalytic values in the β matrix are greater than 0 and the cross catalytic values are smaller than 0. An advantage of GARD is its generality, having less restrictive assumptions. Still, an appreciable symmetry breaking ensues, both for single and for multiple assembly simulations.

Chirality and the origin of life
The C-GARD model predicts a considerable degree of chiral selection, but this happens at relatively high values of the enantiodiscrimination parameter α, as compared to the typical values obtained by a statistical analysis [74] of values reported in CHIRBASE [75]. This may be taken as an indication that assembly-based enantioselection could not have acted at the earliest stages of life's origin. CHIR-BASE is, however, restricted mainly to small, relatively simple molecules. A study quantifying the capacity of random molecular structures to manifest enantioselection has predicted, through the use of random diffusion limited aggregates, a linear relationship between the size of the asymmetric structure and its enantiomeric dis-crimination [55]. This is also in line with analyses using string complementarily models, that generally show growing recognition capacity with increasing molecular sizes [69,70,76,77]. As more information is accumulated on the nature size and complexity of prebiotically-available molecules, better fine-tuning of the prediction of C-GARD will become possible. Furthermore, based on such relations between molecular size and enantiomeric discrimination it may become possible to predict the minimal size of prebiotic molecules that would generate sufficient symmetry breaking as inferred from the C-GARD model. Such analysis cannot however be presently performed based on CHIRBASE information, as this database does not have explicit display of molecular sizes. A considerable number of publications perceive chirality not merely as a central characteristic of life but as a prerequisite for its emergence ( [78] and references therein). This stems from the widely accepted notion that self-replication of biopolymers is essential for life's inception. Indeed, previous experiments [78,79] as well as in theoretical models [73,80] has indicated that a very high degree of chiral purity is required for successful polymerbased information transfer. Consequently, many works assumed that at some point in early earth history physical and/or chemical processes have led to pronounced symmetry breaking, which occurred in an inanimate environment and allowed the initiation of life processes. The alternative scenario presented here involves rudimentarily replicating compositional entities, such as GARD assemblies, independent of biopolymer templating. These afford a potential origin of chiral selection as part of the mutually-catalyzed accretion dynamics of noncovalent molecular assemblies. Once a sufficient degree of chiral selection is achieved, more elaborate informational biopolymers may become possible. Thus, the C-GARD model highlights the possibility that chiral selection is a result of, rather than a prerequisite for early lifelike processes.

The GARD formalisms
The C-GARD model is built upon the GARD kinetic model [30,73]. A GARD molecular assembly, typically assumed to consist of amphiphilic molecules, grows by accretion within an buffered environment containing N G different molecule types, and undergoes a stochastic fission process designed to produce two (potentially similar) daughter assemblies. The assembly is represented by a compositional vector n, such that the component n i depicts the number of molecules of type i within the assembly.
Assembly growth rate is governed by the following set of kinetic equations b L X 17 11 , b L L 17 11 , b L D 17 11 , where k f and k b are the forward and backward reaction rates, ρ is buffered extraneous concentration of all molecule types, N G is the number of different molecule types, the assembly size is , and when the assembly size reaches the value N max we impose a stochastic split generating two progenies of equal size [30]. Equation 1 has an obvious steady state fulfilling at , where * indicates a specific steady state value. We note that this is a stable uniform equilibrium steady state (with all n i equal), which is different from the dynamic quasi-stationary states [21,43] constituting the composomes. The dynamics involving periodic fission events averts the attainment of equilibrium, and induces continuous transitions among quasistationary states typical of GARD dynamics. Such behavior is in fact the result of stochastic small perturbations of the concentrations and rates, corresponding to GARD's life-like characteristics. For evaluation of compositional similarity among different assemblies (e.g. at two different time points) we use the normalized dot product of the corresponding composition vectors (i.e. cosine of an angle between the two composition vectors) [29,30,32]: A similarity threshold of H ≥ 0.95 is used in an iterative procedure to classify an assembly within one of several predefined composomes or, if necessary, to define a new composome. This procedure is akin to that which previously referred to as clustering of multiple instances of composomes into compotypes [81]. Mutual rate enhancement exerted by molecule type j on molecule type i is represented by the non-negative element β ij in an N G × N G matrix (Eq. 1). The choice of rate enhancement distribution characteristics in GARD is guided by an embodiment of the Receptor Affinity Distribution (RAD) formalism for catalytic activities [32,[69][70][71], which is supported experimentally by analyses of immunoglobulin and phage display libraries [69]. The extension of RAD from affinities to catalytic rate parameters [32] derives from the relation between binding and catalysis governed by transition state theory and implies a lognormal distribution for the catalytic intensities β ij : Where μ and σ are respectively the mean and standard deviation of the distribution, and γ is a constant related to the subsite binding energy in the RAD model [69].
In the present embodiment we use a Poisson approximation with a single statistical parameter, λ, interpretable as the average number of successful intermolecular subsite recognition events in the RAD model (μ = λ and σ = λ 1/2 , [32] appendix). Except where otherwise indicated, we use λ = 6 which has been proven appropriate in a study that addresses GARD heritability properties [32].

C-GARD and its symmetry properties
In a C-GARD simulation, half of the entries in the compositional vector represent the D isomers and the other half -the antipodal L isomers (Figure 2), thus leading to the definition of the compositional vector where and are the counts of the two enantiomers of the molecule type i within the assembly. In C-GARD the assembly prefission size is .
We employ a parity principle of space inversion equivalence by requiring that the catalytic interaction coefficient for a given pair of molecules would be equal to that corresponding to their respective enantiomers, resulting in a chiral 2N G × 2N G β matrix (Figure 2). Parity-violating energy difference between enantiomers is excluded from the analysis as it is generally regarded too minute to account for macroscopic behavior [82][83][84][85].
A key property intrinsically associated with chirality is the ability of a chiral compound (e.g. L i ) to differentiate between two encountered enantiomers, L j and D j , ( Figure  2). Such capacity is represented by the value of the enantiodiscrimination factor α.
To assess the degree to which two compositional vectors n 1 and n 2 are antipodal to each other, we define the antipodicity M as the similarity H (Eq. 2) between a composition n 1 and, the antipodal composition of n 2 : where and vice-versa. Note that since M is defined through the measure H, it is also normalized and

Distribution of enantiomeric discrimination
It is necessary to utilize values of that would satisfy a lognormal distribution of α ij , yet would not alter the overall distribution of values within the matrix β. This requirement is fulfilled here by the use of a dummy variable, ε, obeying a Gaussian distribution with a mean of zero and with σ ε <σ (σ 2 being the variance of the values and σ ε 2 the variance for the ε ij values). This variable is used to generate the interactions between L and D isomers according to: whereby one quadrant of the chiral beta matrix (labeled in Figure 2) is generated based on the RAD lognormal distribution with its parameters μ and σ, and a second quadrant (labeled in Figure 2) is obtained via Eq. 6. μ and σ are mean and standard deviation of the lognormal distribution of β ij values. Plotting a probability distribution of α values (Eq. 6) using pairs of β values obtained by this procedure, a best fit to earlier experimental data [74,75] is obtained for σ ε = 0.8 ( Figure 11).

Weak and strong enantiomeric selection
Weak enantioselection is defined as the average absolute value of the enantiomeric excess over a given set of N G molecular types: W W measures the extent of symmetry breaking for a given molecular population, regardless of the direction (L or D) it assumes for each individual molecular types.
On the other hand, strong enantioselection is defined as the average value of the signed enantiomeric excesses for the same N G molecular types: We note that W S is usually used when addressing homochirality in a group of similar compounds, such as amino acids, where it indicates the tendency of such molecular repertoire to have enantiomeric excess in the same direction (say L) for all compounds of a given class. W W is more suitable for analyzing early enantioselection in diverse molecular repertoires, as done in the present work. Figure 12 illustrates the difference between W W and W S for several specific compositions.   The authors say nothing about sign of the elements β ij , but from the Equation 5 (Authors comment: Equation 6 in the final version) it follows that all these elements should be positive. Consider the steady state solution of (1). As the components n i of the compositional vector n are all non-negative, it is easy to see that the expression the square brackets is always larger than 1, so that the steady state compositional vector n* elements n* i can be found from the conditions producing the uniform distribution with Employing the parameter values the authors used for numerical simulations: k f = 5*10 -2 , k b = 5*10 -3 , ρ = 10 -2 , N G = 100 we find that in steady state regime n* I = 0.1 and N = N G n* = 10. This value is much smaller than fission level N max = 200 so that when the steady state is reached the system arrives at the uniform distribution and there is no way to observe the chiral inhomogeneity.

Reviewers' comments
Minor comments: 1. Equation 2 is just a cosine of an angle between two vectors in multidimensional space.  (2) is the stable one. The authors do not consider the stability of the numerical solutions of Equation 1 to small perturbations of the parameters. Because the stability is not addressed one cannot be sure that the produced solutions correspond in any way to the real life situation. Author's response. In the version reviewed by these referees, Eq. 1 was missing the factor N multiplying k f due to a typing error. With this, (2) comes to be and taking for simplicity N* > Nmax/2 we at steady state N = 1000, which is much higher than the fission size used in this work. In addition, text was introduced in the Methods addressing equilibrium steady state vs. quasistationary states away from equilibrium, as well as relating to the effect of small perturbations. We have also made additional corrections to the text as suggested by these referees, including a ~2-fold extension of the background section, and a detailed paragraph thereof addressing the question of liquid-solid interface as a site of prebiotic evolution.

Reviewer report 2
Meir Lahav, Weizmann Institute of Science, Rehovot, Israel. This manuscript suggests a new stochastic model for spontaneous "Mirror Symmetry Breaking" of possible in a pre-biotic environment. Whereas, I am not competent to evaluate the technicality of the mathematical aspects of the model, I can identify some of the advantages that such model proposes.
The message that the model conveys is clearly presented. The "mirror symmetry breaking " phenomenon is an outcome from a more general model of the authors on the origin of life. In variance to some other models, which consider the role played by a single racemic component only, the present one considers the involvement of a complex mixture of racemic α-amino acids, which makes it more realistic to a primeval environment. Such model may enjoy an additional reward provided it can be supported by pertinent experiments.
I suggest to add some additional references, to the ones mentioned in the manuscripts, which deal with other theoretical proposed stochastic models of "Mirror Symmetry Breaking" Decker P. The Origin of molecular Asymmetry through the Amplification of "Stochastic Information" (noise) in Bioids, Open Systems which can exist in Several Steady States. J.Mol. Evol. 1974, 4, 49.

Reviewer report 3
Sergei Maslov, Brookhaven National Laboratory, New-York, USA. The manuscript describes an interesting extension of authors' earlier Graded Autocatalysis Replication Domain (GARD) model. The new model called Chiral-GARD (or C-GARD) is introduced to explain the symmetry breaking between left-or right-handed biomolecules in the modern biosphere. The main conclusion is that a strong symmetry breaking can result from a relatively modest asymmetry in mutual (auto)catalytic activity: that is D-enantiomers that are more likely to catalyze other D-enantiomers than L-enantiomers, while Lenantiomers preferentially catalyze their L-brothers and sisters. BTW, the very term enantiomer should be explained early on in the text for the benefit of uninitiated. This catalytic asymmetry is quantified by the parameter alpha. Authors have some idea about the range and distribution of alpha from the CHIRBASE database (see Fig. 12 of the manuscript (Authors comment: Figure  11 in the final version)). The minimal value of alpha required for symmetry breaking in their model is somewhat larger than the typical value of entries in the CHIR-BASE database. Since this database is marketed to pharmaceutical industry it is mostly limited to relatively small molecules. The fact that enantiomeric discrimination of such small molecules is too weak indicates that the L-D symmetry breaking during the prebiotic evolution must have occurred at a later stage when prebiotic chiral molecules were already larger than entries in the CHIR-BASE. Authors cite a paper (Ref. [48]) (Authors comment: ref. [55] in final version) reporting a linear relationship between molecule's length and the strength of enantiomeric discrimination (value of alpha). Based on this trend is it possible to predict the minimal size of prebiotic molecules that would generate sufficient symmetry breaking in the C-GARD model? Is this linear correlation also present in CHIRBASE data? From Fig. 1 it follows that even for molecules of very modest length chiral isomers outnumber non-chiral ones. This allowed authors to disregard non-chiral isomers in their C-GARD model. I wonder, would their conclusion be qualitatively different if non-chiral isomers were added to the model. To rephrase it, do non-chiral isomers that are spared the mutually exclusive fight between their L-and D-forms get a competitive advantage over their chiral counterparts? When authors introduce their Receptor Affinity Distribution (RAD) formalism it is very easy to miss that it is *the logarithm* of affinity that follows the Poisson distribution. Only my previous interest in theories explaining lognormal distributions of dissociation constants spared me from this confusion. Authors mention that beta has a lognormal distribution in only one inconspicuous place on this page. I suggest authors explicitly mention it when introducing their GARD model and maybe even write a Poisson distribution formula for P(log(beta)). On a similar note, when introducing the Eq.(5) (Authors comment: Eq. 6 in the final version) authors describe all the variables except for mu which has to be traced back to their verbal discussion of the Poisson distribution. In Figure 4 (Authors comment: Figure 5 in the final version) the parameter sigma epsilon goes as high as 10 for lambda = 10. How it can be reconciled with the earlier requirement that sigma_epsilon < = sigma = square root of lambda? Perhaps authors mislabeled the X-axis in this figure which should read sigma? What is the functional form of the distribution of alphas from the CHIRBASE database in Figure 12? Is it indeed lognormal as stated in the beginning of the section "Distribution of enantiomeric discrimination"? Perhaps, in Fig. 12 authors can change axes to log-log (or show a log-log insert) which would let readers verify this fact? Author's response. Text was added in the introduction to clarify the terminology used, including "enantiomers". Additional text is now in place in the discussion to address the correlation between molecular size and enantiodiscrimination. A new paragraph in the first section of the discussion addresses the intriguing question of asymmetric and symmetric molecular mixtures. An explicit formula for the lognormal distribution of the rate enhancement parameters β has been added (present Eq. 3), and μ and σ are defined at this earlier instance. We have indeed mislabeled the X-axis in figure 5 (thanks for seeing this!), and it is now corrected to read σ. A double logarithmic transformation of the data presented in Fig.  11 was added as an insert.