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 and 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 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 and Joyce 2014) and may also not be very efficient on the early Earth. This was the reason for Jafarpour et al. (2015, 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). In that work, Frank (1953) 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. (2015, 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 small compartments, or 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.

Both the model of Frank (1953) that is based on rate equations and the alternative stochastic models discussed above describe what is known as spontaneous chiral symmetry breaking. These mechanisms do not require an external chiral influence, even though a very small chiral influence is always present. Possible candidates for causing a systematic chiral influence include polarized light from the UV radiation in star-forming regions (Bailey et al. 1998; Bailey 2001) or from nearby neutron stars (Bonner 1999) or the weak force in fundamental physics (Hegstrom et al. 1980; Hegstrom 1984; Mason and Tranter 1984). The expected effect is likely to be very small (Bada 1995; Lente 2006). This can be quantified by studying the interplay between spontaneous chiral symmetry breaking and the effect of an external chiral influence (Kondepudi and Nelson 1983, 1985). We expect those conclusions to carry over to the stochastic models as well. To verify this, we allow in some of our models for the presence of a chiral influence in the catalytic effect, as was done in the context of a polycondensation model using non-stochastic rate equations (Brandenburg et al. 2005). They found that, as the fidelity of the autocatalytic reactions is increased, a typical bifurcation diagram emerges. Owing to the effect of an external chiral influence, the diagram becomes slightly asymmetric with respect to positive and negative enantiomeric excess. It is therefore referred to as an imperfect bifurcation, just as it was anticipated by (Kondepudi and Nelson 1983, 1985).

Method

The notion of invoking both autocatalysis and mutual antagonism is based on the governing rate equations for the three reactions

$$ \begin{array}{@{}rcl@{}} A&+&D \stackrel{k_{\mathrm{C}}~}{\longrightarrow} 2D, \end{array} $$
(1)
$$ \begin{array}{@{}rcl@{}} A&+&L \stackrel{k_{\mathrm{C}}~}{\longrightarrow} 2L, \end{array} $$
(2)
$$ \begin{array}{@{}rcl@{}} D&+&L \stackrel{k_{\times}~}{\longrightarrow} 2A, \end{array} $$
(3)

where A denotes an achiral molecule that can autocatalyze at a rate kC 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. Developing an intuitive and simple Monte Carlo method will be 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. (2008, 2009) and Jafarpour et al. (2015, 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

$$ {q}_{n}=([A], [D], [L]), $$
(4)

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 Δqn out of a set of seven possible transitions such that the new state vector qn+ 1 is given by

$$ {q}_{n+1}={q}_{n}+{\Delta} {q}_{n}. $$
(5)

For example, for the three reactions (1)–(3), Δqn 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 q0 = (10,10,10). Suppose now that reaction (1) takes place during the next reaction step, then Δq0 = (− 1,1,0), so in our example we would have in the next step q1 = (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.,

$$ \begin{array}{@{}rcl@{}} A \stackrel{k_{+}}{\longrightarrow} D, \quad &&\quad A \stackrel{k_{+}~}{\longrightarrow} L, \end{array} $$
(6)
$$ \begin{array}{@{}rcl@{}} D \stackrel{k_{-}}{\longrightarrow} A, \quad &&\quad L \stackrel{k_{-}~}{\longrightarrow} A. \end{array} $$
(7)

Spontaneous racemization can be the result of natural degradation; see, e.g., Bada et al. (1970). Spontaneous deracemization was introduced by Sugimori et al. (2008). Note that the Δqn 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+ by kC[D]/N for the first reaction and kC[L]/N for the second one. Here, N = [D] + [L] + [A] 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 cross-inhibition reactions, we have

$$ \begin{array}{@{}rcl@{}} &&A \stackrel{\underrightarrow{~k_{\mathrm{C}}[D]/N~}}{} D, \end{array} $$
(8)
$$ \begin{array}{@{}rcl@{}} &&A \stackrel{\underrightarrow{~k_{\mathrm{C}}[L]/N~}}{} L, \end{array} $$
(9)
$$ \quad\left\{\begin{array}{ll} D \stackrel{\underrightarrow{~k_{\times}([D]+[L])/N~}}{} A,\\ L \stackrel{\underrightarrow{~k_{\times}([D]+[L])/N~}}{} A. \end{array}\right. $$
(10)

The curly bracket to the left on the last two reactions indicates that those two reactions happen at the same time.

To study the effects of an imperfect fidelity and of an external chiral influence on the system, we extend our model in a way that is analogous to the formulation employed by Brandenburg et al. (2005). This is accomplished by changing Eqs. (8) and (9) to

$$ \begin{array}{@{}rcl@{}} &&A \stackrel{\underrightarrow{~k_{\mathrm{C}}(f_{+}[D]+f_{-}[L]+\beta_{D}[A])/N~}}{} D, \end{array} $$
(11)
$$ \begin{array}{@{}rcl@{}} &&A \stackrel{\underrightarrow{~k_{\mathrm{C}}(f_{+}[L]+f_{-}[D]+\beta_{L}[A])/N~}}{} L, \end{array} $$
(12)

where f± = (1 ± f)/2 with 0 ≤ f ≤ 1 being the fidelity (f = 1 for perfect fidelity) and βD/L quantifies the bias of the system toward D or L, respectively (βD = βL = 0 in the absence of an external chiral influence on the system).

We recall that in this paper we regard the coefficients k+, k, kC, 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

$$ k_{+} +k_{-} +k_{\mathrm{C}}+k_{\times}\leq1. $$
(13)

We summarize all the reactions modelled in this paper in Table 1.

Table 1 Summary of all seven reactions (6)–(10)

The coefficients k+, k, kC, and k× are used to determine seven intervals, w1, w2, ..., w7, with w1 + w2 + ... + w7 = 1. The values of wi 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 w1 = w2 = k+/2. Likewise, the two reactions in (10) occur with equal probability, so we set w3 = w4 = k/2. Finally, we have w5 = kC[D]/N, w6 = kC[L]/N, and w7 = k×([D] + [L])/N. Since ([D] + [L])/N is, in general, less than unity, it may be possible that no reaction occurs during a particular step. This is true even if βD≠ 0 and βL≠ 0. (Both βD and βL are below unity.)

We denote the interval boundaries by ri, so we have 0 ≤ r1r2 ≤ ... ≤ r7 ≡ 1. At each 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} \leq r < r_{i}\quad (\text{for reaction } i \text{ with } i=1, 2, ..., 7). $$
(14)

The widths of the seven intervals, riri− 1 = wi, 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

$$ \langle A\rangle+\langle D\rangle+\langle L\rangle=1. $$
(15)

Within each population, the enantiomeric excess is defined as

$$ \eta=\frac{[D]-[L]}{[D]+[L]}. $$
(16)

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 primary interest. To determine for which parameters the bifurcation occurs, we monitor the mean enantiomeric excess 〈|η|〉, which is close to zero if there are about equally many molecules of the d and of the l forms.

Altogether, our models have five parameters: k+, k, kC, 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 with 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 CodeFootnote 1, a publicly available time stepping code that is designed to perform computations on massively parallel computers. By default, 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. A first order time step of length unity is chosen to reproduce a discrete reaction step. Furthermore, the derivative module noderiv is selected in the code, which means that no extra ghost zones are used, which would only be needed in a model with spatial extent.

Results

Strategy

We begin by studying the familiar case of autocatalysis and enantiomeric cross-inhibition, i.e., we choose the value of k× and, since k+ = k = 0 in this case, we set kC = 1 − k×, so the bound given by Eq. (13) 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. (2015, 2017), where we set k+ = k× = 0 and vary kC such that kC + k = 1. This is referred to as model II. The case proposed by Sugimori et al. (2008) corresponds to k = kC = 0, so we vary k+ and adjust k× = 1 − k+. This is model III. Finally, we consider the case kC = k× = 0, vary k+, and adjust k = 1 − k+. This is our model IV. Examples of evolutionary tracks are displayed in Fig. 1 for all four models. We also consider intermediate cases that we refer to as I/III and II/IV, where we have considered three non-vanishing probabilities, and a model V, where all four probabilities are non-vanishing; see Table 2 for an overview of all the seven models. For models I, II, III, and 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. For model V, we only consider one case where k+ = k = k× = 0.2 and kC = 0.4.

Fig. 1
figure 1

Examples of the evolution of the enantiomeric excess (e.e.) for cases I, II, III, and IV. Upper left: model I for kC = 1 − k× = 0.1; upper right: model II for kC = 1 − k = 0.9; lower left: model III for k+ = 1 − k× = 0.3; lower right: model II for k+ = 1 − k = 0.3

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

Models I, II, III, and IV

In models I, II, and 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 unsigned enantiomeric excess, 〈|η|〉. It will be close to unity if most of the achiral A molecules of the substrate are turned into D or L. Thus, we expect 〈A〉 to vary inversely with 〈|η|〉. 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 Fig. 2 shows a maximum for k+ = 1 − k = 0.5. However, this only happens after a significant number of reaction steps.

Fig. 2
figure 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, II, III, and IV

To characterize the necessary number of reaction steps for different parameter combinations, we define the parameters nη and nA as the reaction step after which 〈|η|〉 and 〈A〉, respectively, have reached their final values within 1%. In Fig. 3, we plot these numbers as a function of parameters for different population sizes.

For models III and IV, nη can be extremely large – of the order of 108 – especially when the population size is large. For model IV with N = 30,000, for example, even nη = 109 is reached when k+ → 1 and k→ 0. An increasing trend of nη with k+ is also found for model III, i.e., the model of Sugimori et al. (2008), again for large population sizes. In those cases, there is also a dramatic drop in nA when k+ → 1, which shows that most of the molecules are either of the d or of the l form.

Fig. 3
figure 3

Typical reaction numbers nη (black lines) and nA (red lines) as a function of parameters for models I, II, III, and IV for N = 3000 (solid lines) and N = 300 (dotted lines)

For models I and II, on the other hand, both nη and nA are significantly smaller – typically between 103 and 106. Nevertheless, we see again an increasing trend of nη and a decreasing trend of nA when kC 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 (kC) and spontaneous deracemization (k+). In other words, can we trade some fraction of autocatalytic reactions for deracemization and still have the same or an even stronger effect on achieving homochirality?

Looking at Fig. 4, we see that for model I/III, the values of 〈|η|〉 and 〈A〉 are unchanged as k+ = 0.4 − kC is changed, but both nη and nA increase as k+ is increased and thus kC 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 k+, while nA remains of the order of 104, except for some departures at k+ = 0.5.

Fig. 4
figure 4

Bifurcation diagrams of 〈|η|〉 (black) and 〈A〉 (red) (upper row) and time scales nη (black lines) and nA (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 kC = 0.4 − k+ in each case

Model V with Finite Fidelity and an External Chiral Influence

It may have appeared strange that in Figs. 2 and 4, 〈|η|〉 was always close to unity – even for rather small values of kC. The reason is that, until now, autocatalysis was modelled with perfect fidelity (f = 1), i.e., the presence of molecules of the d form always favors the production of more molecules of the d form and not of the l form, and vice versa. This is different when f < 1 in Eqs. (11) and (12).

Models with f < 1 were first studied by Sandars (2003) in a fairly detailed polycondensation model. Later, Brandenburg et al. (2005) also allowed for the presence of a small bias (βD≠ 0 or βL≠ 0) in such a model. In those cases, a bifurcation diagram emerges, where 〈η〉≠ 0 for f = 0. This is referred to as an imperfect bifurcation, because it is asymmetric with respect to positive and negative values of 〈η〉.

We employ what we call here model V, where, in addition to autocatalysis and enantiomeric cross-inhibition, also spontaneous racemization and deracemization are included. Adding these two effects has the advantage that they cause the system to be well “mixed”. Without this, even the case of f = 0 can in some cases lead to a gradual loss of all achiral molecules, which is unrealistic. This problem is alleviated by having k+≠ 0 and k≠ 0.

The result for model V is shown in Fig. 5, where we plot not only 〈|η|〉, but now also 〈η〉, which is the average of the signed enantiomeric excess over all populations. Unlike 〈|η|〉, this quantity can approach + 1 or − 1 only if a large number of independent populations or geneses produce the same chirality. We recall that the number of populations is still 512. The number of members of each population is here taken to be 300. This model shows an imperfect bifurcation at f ≈ 0.2. There is a small neighborhood around this point where for small enough perturbations only positive values of 〈η〉 are possible.

Realistic values of βD or βL could be in the range 10− 12 to 10− 17; see Kondepudi and Nelson (1985), who demonstrated that even such a small chiral influence can tip the balance systematically in one direction. This ignores, however, the effects of fluctuations associated with small system sizes. In practice, this would always dominate over thermal fluctuations. It would therefore remain an exciting prospect to look for biomolecules of unconventional chirality, for example, through metabolic experiments on Mars, as proposed by Sun et al. (2009).

Fig. 5
figure 5

(a) 〈|η|〉 (red) together with 〈A〉 (black) versus f and (b) 〈η〉 versus f for model V with N = 300, k+ = k = k× = 0.2, kC = 0.4, βD = 10− 3, and βL = 0. Note the slight asymmetry for positive and negative values of 〈η

Conclusions

For many decades, the paradigm of Frank (1953) of producing homochirality by a combination of both autocatalysis and mutual antagonism or enantiomeric cross-inhibition has shaped much of our thinking for many decades. It is now clear that the implied necessity of the combined presence of both of these ingredients may not strictly hold. Earlier numerical simulations of Toxvaerd (2013) have already hinted at such a possibility, which may well be a viable one in small enough systems where fluctuations can be important. Two separate aspects of this were already studied in some detail by Sugimori et al. (2008, 2009) and Jafarpour et al. (2015, 2017).

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. (2008, 2009) and Jafarpour et al. (2015, 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 and Lineweaver 2005) and that concentrations may have been low. Although some of the processes reported above may involve altogether up to a billion reaction steps, this may not be long on geochemical time scales and would correspond to only about 40,000 years if we assumed a reaction time of 20 minutes. 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 and Brandenburg 2010; Hochberg et al. 2017), which may be strongly affected by fluctuations having either an enhancing or diminishing effect on achieving homochirality.