The limited roles of autocatalysis and enantiomeric cross inhibition in achieving homochirality in dilute systems

To understand the effects of fluctuations on achieving homochirality, we employ a Monte-Carlo method where autocatalysis and enantiomeric cross-inhibition, as well as racemization and deracemization reactions are included. The results of earlier work either without autocatalysis or without cross-inhibition are reproduced. Bifurcation diagrams and the dependencies of the number of reaction steps on parameters are studied. In systems with 30,000 molecules, for example, up to a billion reaction steps may be needed to achieve homochirality without autocatalysis.


Introduction
There are many reasons why the chemistry of life is based on carbon.The ability to form complex macromolecules is one of them (Rothery et al., 2008;Longstaff, 2014).Many carbon-bearing molecules also have the property of being chiral, i.e., the three-dimensional structure of such a molecule is different from its mirror image (Avetisov et al., 1991).Even relatively simple amino acids such as alanine, arginine, and proline are such examples.In solution, these molecules rotate polarized light to the left, which is why they are called levorotatory, so we say they are of the l form.Only rarely does life on Earth involve amino acids of the d form (Milner-White, & Russell, 2005).Those molecules are dextrorotatory and would rotate polarized light in the opposite direction.Naturally occurring sugars such as ribose in the backbone of DNA and RNA are also chiral and of the d form.
Abiogenically produced carbon-bearing chiral molecules always occur as a mixture of d and l forms.One calls such mixtures racemic, i.e., the molecules are achiral as opposed to chiral, i.e., the mixture is not homochiral and molecules of both the d and l forms occur.However, synthesizing chiral polymers such as proteins and nucleic acids only succeeds in the predominately isotactic form.This realization goes back to the experiments of Joyce et al. (1984) using template-directed polycondensation.This work demonstrated what was called enantiomeric cross-inhibition, i.e., the termination of polycondensation with enantiomers of the opposite chirality.More recent work, however, showed that enantiomeric cross-inhibition may not always operate (Sczepanski & Joyce, 2014) and may also not be very efficient on the early Earth.This was the reason for Jafarpour et al. (2015Jafarpour et al. ( , 2017) ) to question the long-held belief in the necessity of enantiomeric cross-inhibition for producing homochirality based on the famous proposal of Frank (1953).He suggested that the interplay of autocatalysis combined with what he called mutual antagonism would be a necessary ingredient for reaching complete homochirality.Jafarpour et al. (2015Jafarpour et al. ( , 2017) ) invoked the possibility that local fluctuations could yield results that are different from those expected by solving the kinetic rate equations.Fluctuation effects are generally expected to play a role in dilutes systems, for example in small compartments and when concentrations are low.
In an earlier paper, Sugimori et al. (2008) did already invoke the effects of fluctuations and questioned the concept of autocatalysis.Indeed, the number of known autocatalytic reactions is small -with the reaction of Soai et al. (1995) being basically still the only one.However, the more general case of arbitrary combinations with and without autocatalysis and with and without enantiomeric cross-inhibition has not yet been studied.Such a more general approach is of interest in view of numerical studies such as that of Toxvaerd (2013), where homochirality has been found without apparent autocatalysis or enantiomeric cross-inhibition.

Method
The notion of invoking both autocatalysis and mutual antagonism is based on the governing rate equations for the three reactions where A denotes an achiral molecule that can autocatalyze at a rate k C to a chiral one either of the d or the l form, while D and L can cross-inhibit into an achiral one at a rate k × .However, rate equations no longer provide a suitable description of the relevant kinetics when the system is dilute and reactions are rare (Gillespie, 1977;Toxvaerd, 2014).In that case, a stochastic approach must be adopted, which is the goal of the present paper.
It will be advantageous to regard the reaction rates as probabilities for the respective reactions to occur within a given reaction step.Thus, instead of solving the underlying master equation, as was done by Sugimori et al. (2008Sugimori et al. ( , 2009) ) and Jafarpour et al. (2015Jafarpour et al. ( , 2017)), we adopt here a Monte-Carlo approach in which at each reaction step n, one of several possible reactions occur (Verlet, 1967).The state of the system is governed by the state vector where squared brackets denote the concentrations of the respective molecules, and n denotes the reaction step.At each reaction step n, we select a transition ∆q n out of a set of seven possible transitions such that the new state vector q n+1 is given by For example, for the three reactions (1)-( 3), ∆q n could be one of the three vectors (−1, 1, 0), (−1, 0, 1), or (2, −1, −1).
Let us illustrate the formalism using an example.Initially, at n = 0, the system has, say, 10 molecules of each of the three possible ones (A, D, and L), so we have q 0 = (10, 10, 10).Suppose now that the reaction (1) takes place during the next reaction step, then ∆q 0 = (−1, 1, 0), so in our example we would have in the next step q 1 = (9, 11, 10).This means that one out of the 10 molecules of the d form reacted with an achiral one A to produce two new molecules of the d form.One of the A molecules got removed, so only 9 of them are left.Also one of the D molecules got removed, but since one of them did already participate in the reaction, only one more of them has occurred after the reaction, i.e., we now have 11 molecules of the d form.Note that the model is mass conserving, i.e., the total number of molecules is constant and always equal to the initial value N .No polycondensation reactions are considered here.
In addition to the reactions discussed above, we also allow for spontaneous racemization and deracemization reactions, i.e., Note that the ∆q n of the pair of reactions in Eqs. ( 1) and ( 2) is the same as in Eq. ( 6).The former two reactions, which are autocatalytic, can also be written in a form similar to Eqs. ( 6) and ( 7) by replacing k /N in those two reactions, respectively.Here, N is the total number of all molecules, regardless of whether they are chiral (D or L) or achiral (A).
Likewise, the reaction (3) corresponds to racemization reactions (6) with the coefficients k − in the two parts of that equation being replaced by k × [L]/N and k × [D]/N , respectively.Thus, for the autocatalytic and enantiomeric crossinhibition reactions, we have The curly bracket to the left on the last two reactions indicates that those two reactions happen at the same time.
We recall that in this paper we regard the coefficients k + , k − , k C , and k × not as rates, but as probabilities.The sum of all probabilities must then not exceed unity, i.e., the four coefficients must obey the constraint We summarize all the reactions modelled in this paper in Table 1.The coefficients k + , k − , k C , and k × are used to determine seven intervals, w 1 , w 2 , ..., w 7 , with w 1 + w 2 + ... + w 7 = 1.The values of w i are equal to the respective probabilities of the seven possible reactions ( 6)-( 10).The two reactions in (6) occur with equal probability, so we set w 1 = w 2 = k + /2.Likewise, the two reactions in (10) occur with equal probability, so we set is, in general, less than unity, it may be possible that no reaction occurs during a particular step.
We denote the interval boundaries by r i , so we have 0 ≤ r 1 ≤ r 2 ≤ ... ≤ r 7 ≡ 1.At reach reaction step, we generate a new random number r with 0 ≤ r ≤ 1, where the value of r determines which of the various reactions occurs.Reaction i is performed when r i−1 ≤ r < r i (for reaction i with i = 0, 1, ..., 7). ( Table 1 Summary of all seven reactions (6)-(10).
The widths of the seven intervals, r i − r i−1 = w i , are listed in Table 1.
To make statistically meaningful statements, we have to perform sufficiently many experiments, and then consider averages over all experiments.In practice, we perform many experiments at the same time and refer to them as populations.The populations are completely independent of each other.We then compute normalized averages over all populations, denoted by angle brackets, so A , D , and L are the fractional concentrations of each of the three types of molecules.Therefore, we have Within each population, the enantiomeric is defined as It can be between −1 and +1 for homochiral states and is close to zero for a racemic mixture.From a statistical point of view, however, the two homochiral states, η = ±1, are equivalent, so only the average of the modulus of η is of interest.To determine for which parameters the bifurcation occurs, we monitor the mean enantiomeric excess, |η| , which is close to zero if most of the members of a population are and remain achiral.Altogether, our models have five parameters: k + , k − , k C , k × , and the population size N .The number of populations is an additional parameter that is chosen to be 512 in all cases.The models are completely symmetric this respect to D and L, i.e., there is no preference with respect to D and L.
We report here the results of numerical experiments performed with the Pencil Code1 , a publicly available time stepping code that is designed to perform computations on massively parallel computers.It uses the third order Runge-Kutta time stepping scheme of Williamson (1980).Normally, the code is used for solving partial differential equations using meshpoints.Here, however, no spatial extent will be considered and each mesh point can be regarded as an independent population.This is how many populations can then be solved for in parallel.The time step is chosen to be unity to reproduce a discrete reaction step.Furthermore, the derivative module noderiv is in the code, which means that no extra ghost zones are used, which would only be needed in a model with spatial extent.

Strategy
We begin by studying the familiar case of autocatalysis and enantiomeric crossinhibition, i.e., we choose the value of k × and, since k + = k − = 0 in this case, we set k C = 1−k × , so the bound given by Eq. ( 11) is saturated.We refer to this as model I, which corresponds to that of Frank (1953).As already explained above, each reaction step can correspond to a different time interval, which does not need to be specified for our present purpose.Next, we consider the case proposed by Jafarpour et al. (2015Jafarpour et al. ( , 2017)), where we set k + = k × = 0 and vary k C such that k C + k − = 1.This is referred to as model II.The case proposed by Sugimori et al. (2008) corresponds to k − = k C = 0, so we vary k + and adjust k × = 1 − k + .This is model III.Finally, we consider the case k C = k × = 0, vary k + , and adjust k − = 1 − k + .This is our model IV.We also consider intermediate cases that we refer to as I/III and II/IV, where we have considered three non-vanishing probabilities; see Table 2 for an overview of all the six models.For models I-IV, we vary p with 0 < p < 1, while for models I/III and II/IV, we vary q, but now in the range 0 < q < p, keeping p = 0.4 fixed.When q = 0, model I/III is identical to model I, while for q = p, model I/III is identical to model III.Likewise, model II/IV is identical to model II for q = 0 and identical to model IV for q = p.

Models I-IV
In models I-III, there is a bifurcation of solutions within each population to either +1 or −1, depending on chance.To determine for which parameters a bifurcation occurs, we monitor the mean enantiomeric excess, |η| , which is  close to zero if most of the members of a population are and remain achiral.Thus, we expect A two vary complimentary to |η| .This is indeed the case.For model IV, however, we see that A gradually changes from −1 for k + ≤ 0.2 to A = 0 for k + ≥ 0.8.Thus, for large values of k + , there are hardly any achiral molecules left, but this does not mean that the final stage is chiral.Instead, there are large fluctuations among different populations, where some of the molecules can be close to fully chiral, with many of them being of the d form in some populations, and many of them being of the l form in other populations.This is particularly the case for smaller populations with, e.g., 300 or 3000 members.For 30, 000 members, on the other hand, some strongly chiral states are still possible when k + is not too large.Indeed, we see that the black dashed line in Figure 2 shows a maximum for k + = 1 − k − = 0.5.However, this only happens after a significant number of reaction steps.
To characterize the necessary number of reaction steps for different parameter combinations, we define the parameters n η and n A as the reaction step after which |η| and A , respectively, have reached their final values within 1%.In Figure 3, we plot these times as a function of parameters for different population sizes.
We see that for model IV, n η can be extremely large -of the order of 10 8 -especially when the population size is large.For n = 30, 000, for example, even n η = 10 9 is found when k + → 1 and k − → 0. An increasing trend of n η with k + is also found for model II, i.e., the one of Sugimori et al. (2008), again for large population sizes.In those cases, there is also a dramatic drop in n A when k + → 1, which shows that most of the molecules are either of the d or of the l form.
For models I and II, on the other hand, both n η and n A are significantly smaller -typically between 10 3 and 10 6 .Nevertheless, we see again an increasing trend of n η and a decreasing trend of n A when k C is increased, which is analogous to the increase with k + in models III and IV.

Models I/III and II/IV
The models I/III and II/IV were designed to assess the possibility of cooperative effects between autocatalysis (k C ) and spontaneous deracemization (k + ).In other words, can we trade some fraction of autocatalytic reactions for deracemization and still have the same, of an even stronger, effect on achieving homochirality?
Looking at Figure 4, we see that for the model I/III, the values of |η| and A are unchanged, but both n η and n A increase as k + is increased and thus k C decreased.This suggests that there is no cooperative effect of spontaneous deracemization on autocatalysis, because it now takes more steps to achieve the same result.For model II/IV, we see that A is again unchanged, but now |η| decreases.Thus, replacing some of the autocatalytic reactions by deracemization has a detrimental effect.The values of n η increase with increasing values of k + , while n A remains of the order of 10 4 , except for some departures at k + = 0.5.

Conclusions
In this work, we have presented a unified approach to homochirality by considering all possible combinations discussed above: with and without autocatalysis, and with and without enantiomeric cross-inhibition.In most of the cases, we have allowed for population sizes between 30 and 30,000 members.Our approach allows us to understand the earlier work of Sugimori et al. (2008Sugimori et al. ( , 2009) ) and Jafarpour et al. (2015Jafarpour et al. ( , 2017) ) in the broader context of models that include these two models as special cases.In fact, the close relation between these approaches does not seem to have been broadly recognized yet.The only paper that mentions both of them is the recent review of Walker (2017).
The unified approach to stochastic effects in chemical systems discussed in the present paper is conceptually simple and can easily be generalized to other systems, for example those with additional spatial extent.Such an approach is particularly important to astrobiology and the origin of life, given that independent geneses may have occurred at different locations on the early Earth (Davies & Lineweaver, 2005) and that concentrations may have been low.Independent geneses may yield opposite chiralities at different locations on the early Earth (Brandenburg and Multamäki, 2004).Another possible extension of the present approach is to include the more general case of polymerization or polycondensation reactions of nucleotides (Sandars, 2003;Brandenburg et al., 2005) and of peptides (Plasson et al., 2004;Brandenburg et al., 2007).The polycondensation reactions can easily be included and would simply increase the number of possible reactions from seven in the present work to any arbitrary number.Particularly useful would be the study of network catalysis (Plasson & Brandenburg, 2010), which may be strongly affected by fluctuations having either an enhancing or diminishing effect on achieving homochirality.

Fig. 2
Fig. 2 Bifurcation diagrams of |η| (black) and A (red) for N = 3000 (solid lines) and N = 300 (dotted lines) as a function of parameters for models I-IV.

Fig. 3
Fig.3Typical reaction numbers nη (black lines) and n A (red lines) as a function of parameters for models I-IV for N = 3000 (solid lines) and N = 300 (dotted lines).

Fig. 4
Fig. 4 Bifurcation diagrams of |η| (black) and A (red) (upper row) and time scales nη (black lines) and n A (red lines) (lower row) for N = 3000 (solid lines) and N = 300 (dotted lines) for the mixed cases I/III (left, for k × = 0.6) and II/IV (right, for k − = 0.6) as a function of k + , keeping k C = 0.4 − k + in each case.

Table 2
Summary of the six sets of experiments discussed in this paper.The columns "auto" and "inhib" indicate where autocatalysis and enantiomeric cross-inhibition are possible.