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

TL;DR

This study uses Monte Carlo simulations to investigate how autocatalysis and enantiomeric cross-inhibition influence the emergence of homochirality in dilute systems, revealing their limited roles and the extensive steps needed without autocatalysis.

Contribution

The paper demonstrates through simulations that autocatalysis and cross-inhibition have limited effectiveness in achieving homochirality, especially in dilute systems, requiring many reaction steps.

Findings

Homochirality is hard to achieve without autocatalysis in dilute systems.

Up to a billion reaction steps may be needed without autocatalysis.

Bifurcation diagrams reveal parameter dependencies in the process.

Abstract

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.

Equations

$$A+D$$

$$A+L$$

$$D+L$$

$${\mathbf{q}}_n=([A],[D],[L]),$$

$${\mathbf{q}}_{n+1}={\mathbf{q}}_n+\Delta {\mathbf{q}}_n.$$

$$A\stackrel{k_{+}}{⟶}D,$$

$$D\stackrel{k_{-}}{⟶}A,$$

$$A\stackrel{\underrightarrow{k_{\mathrm{C}}[D]/N}}{}D,$$

$$A\stackrel{\underrightarrow{k_{\mathrm{C}}[L]/N}}{}L,$$

$$\{\begin{array}{l}D\stackrel{\underrightarrow{k_{\times }([D]+[L])/N}}{}A,\\ L\stackrel{\underrightarrow{k_{\times }([D]+[L])/N}}{}A.\end{array}$$

$$A\stackrel{\underrightarrow{k_{\mathrm{C}}(f_{+}[D]+f_{-}[L]+{\beta }_D[A])/N}}{}D,$$

$$A\stackrel{\underrightarrow{k_{\mathrm{C}}(f_{+}[L]+f_{-}[D]+{\beta }_L[A])/N}}{}L,$$

$$k_{+}+k_{-}+k_{\mathrm{C}}+k_{\times }\le 1.$$

$$r_{i-1}\leq r<r_{i}\quad\mbox{(for reaction $i$ with $i=1$, $2$, ..., $7$)}.$$

$$⟨A⟩+⟨D⟩+⟨L⟩=1.$$

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

Full text

1 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. (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 & 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 & 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 & Nelson (1983, 1985).

2 Method

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

[TABLE]

where $A$ denotes an achiral molecule that can autocatalyze at a rate $k_{\rm 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_{\times}$. 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

[TABLE]

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 $\Delta\bm{q}_{n}$ out of a set of seven possible transitions such that the new state vector $\bm{q}_{n+1}$ is given by

[TABLE]

For example, for the three reactions (1)–(3), $\Delta\bm{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 $\bm{q}_{0}=(10,\,10,\,10)$. Suppose now that reaction (1) takes place during the next reaction step, then $\Delta\bm{q}_{0}=(-1,\,1,\,0)$, so in our example we would have in the next step $\bm{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.,

[TABLE]

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 $\Delta\bm{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_{+}$ by $k_{\rm C}[D]/N$ for the first reaction and $k_{\rm C}[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_{\times}[L]/N$ and $k_{\times}[D]/N$, respectively. Thus, for the autocatalytic and enantiomeric cross-inhibition reactions, we have

[TABLE]

[TABLE]

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

[TABLE]

where $f_{\pm}=(1\pm f)/2$ with $0\leq f\leq 1$ being the fidelity ($f=1$ for perfect fidelity) and $\beta_{D/L}$ quantifies the bias of the system toward $D$ or $L$, respectively ($\beta_{D}=\beta_{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_{-}$, $k_{\rm C}$, and $k_{\times}$ 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

[TABLE]

We summarize all the reactions modelled in this paper in Table LABEL:Tab.

The coefficients $k_{+}$, $k_{-}$, $k_{\rm C}$, and $k_{\times}$ 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 $w_{3}=w_{4}=k_{-}/2$. Finally, we have $w_{5}=k_{\rm C}[D]/N$, $w_{6}=k_{\rm C}[L]/N$, and $w_{7}=k_{\times}([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 $\beta_{D}\neq 0$ and $\beta_{L}\neq 0$. (Both $\beta_{D}$ and $\beta_{L}$ are below unity.)

We denote the interval boundaries by $r_{i}$, so we have $0\leq r_{1}\leq r_{2}\leq...\leq r_{7}\equiv 1$. At reach reaction step, we generate a new random number $r$ with $0\leq r\leq 1$, where the value of $r$ determines which of the various reactions occurs. Reaction $i$ is performed when

[TABLE]

The widths of the seven intervals, $r_{i}-r_{i-1}=w_{i}$, are listed in Table LABEL:Tab.

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 $\langle A\rangle$, $\langle D\rangle$, and $\langle L\rangle$ are the fractional concentrations of each of the three types of molecules. Therefore, we have

[TABLE]

Within each population, the enantiomeric excess is defined as

[TABLE]

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, $\eta=\pm 1$, are equivalent, so only the average of the modulus of $\eta$ is of primary interest. To determine for which parameters the bifurcation occurs, we monitor the mean enantiomeric excess $\langle|\eta|\rangle$, 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_{-}$, $k_{\rm C}$, $k_{\times}$, 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 Code111https://github.com/pencil-code, DOI:10.5281/zenodo.2315093, 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.

3 Results

3.1 Strategy

We begin by studying the familiar case of autocatalysis and enantiomeric cross-inhibition, i.e., we choose the value of $k_{\times}$ and, since $k_{+}=k_{-}=0$ in this case, we set $k_{\rm C}=1-k_{\times}$, 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_{\times}=0$ and vary $k_{\rm C}$ such that $k_{\rm C}+k_{-}=1$. This is referred to as model II. The case proposed by Sugimori et al. (2008) corresponds to $k_{-}=k_{\rm C}=0$, so we vary $k_{+}$ and adjust $k_{\times}=1-k_{+}$. This is model III. Finally, we consider the case $k_{\rm C}=k_{\times}=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, and a model V, where all four probabilities are non-vanishing; see Table LABEL:Tab2 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_{\times}=0.2$ and $k_{\rm C}=0.4$.

3.2 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, $\langle|\eta|\rangle$. 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 $\langle A\rangle$ two vary inversely with $\langle|\eta|\rangle$. This is indeed the case. For model IV, however, we see that $\langle A\rangle$ gradually changes from $1$ for $k_{+}\leq 0.2$ to $\langle A\rangle=0$ for $k_{+}\geq 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_{\eta}$ and $n_{A}$ as the reaction step after which $\langle|\eta|\rangle$ and $\langle A\rangle$, respectively, have reached their final values within $1\%$. In Figure 3, we plot these numbers as a function of parameters for different population sizes.

For models III and IV, $n_{\eta}$ can be extremely large – of the order of $10^{8}$ – especially when the population size is large. For model IV with $N=30,000$, for example, even $n_{\eta}=10^{9}$ is reached when $k_{+}\to 1$ and $k_{-}\to 0$. An increasing trend of $n_{\eta}$ 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 $n_{A}$ when $k_{+}\to 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_{\eta}$ and $n_{A}$ are significantly smaller – typically between $10^{3}$ and $10^{6}$. Nevertheless, we see again an increasing trend of $n_{\eta}$ and a decreasing trend of $n_{A}$ when $k_{\rm C}$ is increased, which is analogous to the increase with $k_{+}$ in models III and IV.

3.3 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_{\rm C}$) 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 Figure 4, we see that for model I/III, the values of $\langle|\eta|\rangle$ and $\langle A\rangle$ are unchanged as $k_{+}=0.4-k_{\rm C}$ is changed, but both $n_{\eta}$ and $n_{A}$ increase as $k_{+}$ is increased and thus $k_{\rm 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 $\langle A\rangle$ is again unchanged, but now $\langle|\eta|\rangle$ decreases. Thus, replacing some of the autocatalytic reactions by deracemization has a detrimental effect. The values of $n_{\eta}$ increase with increasing $k_{+}$, while $n_{A}$ remains of the order of $10^{4}$, except for some departures at $k_{+}=0.5$.

3.4 Model V with finite fidelity and an external chiral influence

It may have appeared strange that in Figs. 2 and 4, $\langle|\eta|\rangle$ was always close to unity – even for rather small values of $k_{\rm C}$. 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 ($\beta_{D}\neq 0$ or $\beta_{L}\neq 0$) in such a model. In those cases, a bifurcation diagram emerges, where $\langle\eta\rangle\neq 0$ for $f=0$. This is referred to as an imperfect bifurcation, because it is asymmetric with respect to positive and negative values of $\langle\eta\rangle$.

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_{+}\neq 0$ and $k_{-}\neq 0$.

The result for model V is shown in Figure 5, where we plot not only $\langle|\eta|\rangle$, but now also $\langle\eta\rangle$, which is the average of the signed enantiomeric excess over all populations. Unlike $\langle|\eta|\rangle$, 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\approx 0.2$. There is a small neighborhood around this point where for small enough perturbations only positive values of $\langle\eta\rangle$ are possible.

Realistic values of $\beta_{\rm D}$ or $\beta_{\rm L}$ could be in the range $10^{-12}$ to $10^{-17}$; see Kondepudi & Nelson (1985), who demonstrated that even such a small chiral influence good 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).

4 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 & 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 & Brandenburg, 2010; Hochberg et al., 2017), which may be strongly affected by fluctuations having either an enhancing or diminishing effect on achieving homochirality.

Acknowledgements

The author thanks Hannah Zona for comments and encouragement. I am indebted to the two reviewers for their constructive remarks and suggestions. This work was supported through the National Science Foundation, grant AAG-1615100, the University of Colorado through its support of the George Ellery Hale visiting faculty appointment, and the grant “Bottlenecks for particle growth in turbulent aerosols” from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048. The simulations were performed using resources provided by the Swedish National Infrastructure for Computing (SNIC) at the Royal Institute of Technology in Stockholm and Chalmers Centre for Computational Science and Engineering (C3SE).