Reaction-noise induced homochirality David Hochberga, , Mar´ıa-Paz Zorzanoa, ∗ aCentro de Astrobiolog´ıa (CSIC-INTA), Ctra. Ajalvir Km. 4, 28850 Torrejo´n de 7 0 Ardoz, Madrid, Spain 0 2 n a Abstract J 3 Starting from the chemical master equation, we employ field theoretic techniques ] to derive Langevin-type equations that exactly describe the stochastic dynamics E of the Frank chiral amplification model with spatial diffusion. The intrinsic multi- P plicative noise properties are completely and rigorously derived by this procedure. . o We carry out numerical simulations in two spatial dimensions. When the inherent i b spatio-temporal fluctuations are properly included, then complete chiral amplifica- - q tion results from a purely racemic initial configuration. Phase separation can also [ arise in which the enantiomers coexist in spatially segregated domains separated by 1 a sharp racemic interface or boundary. v 5 0 0 1 0 7 1 Introduction 0 / o i b Questions about the origin of life have always held great fascination. One - q especially importantproblemistoexplainchiralsymmetry breakinginnature, : v why for example, it came to be that the nucleotide links of RNA and DNA i X incorporate exclusively dextro-rotary (D) ribose and D-deoxyribose while the r enzymes involve only laevo-rotary (L) enantiomers of amino acids. Mirror a symmetry is broken in the bioorganic world and life as we know it is invariably linked with homochirality. For recent reviews that survey existing hypotheses concerning this phenomenon, specific experimental realizations and additional references, see [1,2,3]. A simple chemical scheme to explain spontaneous generation of chiral asym- metry was proposed by Frank in 1953 [4]. In it original version, the Frank model leads to an amplified production of a chiral enantiomer from a racemic Corresponding author. ∗ Email addresses: [email protected] (David Hochberg), [email protected] (Mar´ıa-Paz Zorzano). Preprint submitted to Elsevier Science 9 February 2008 mixture seeded with an initial bias. If however the initial state of the system is strictly racemic, it remains racemic for all time. The problem is to explain the origin of the initial chiral seed. The purpose of this Letter is to demonstrate that when the inherent spatio-temporal fluctuations are properly included, then complete chiral amplification results from a purely racemic initial con- figuration. Phase separation can also arise in which the enantiomers coexist in spatially segregated domains separated by a racemic interface or bound- ary. The internal fluctuations responsible for these two a-priori unpredictable outcomes are a natural consequence of the fact that chemical reactions are generally diffusion-limited, and take place in imperfectly mixed systems [5]. ThegeneralmethodweemployinthisLetterisbestappreciatedbythoroughly working through a simple example based on an extension of the Frank model, which we now introduce. Let L and D denote a pair of enantiomers. Then the modified Frank model we will study is described by the following reaction scheme [6]: Autocatalytic production: ⇋k1 ⇋k1 L+A L+L, D+A D+D (1) k k 3 3 Mutual destruction, or dimerization, in a second order reaction: L+D k2 P. (2) −→ The k denote the rate constants and we take the achiral substance A as a i uniform constant background. The difference between this and the original Frank model1 lies in the open-flow reactor nature of the process and the fact that the reaction (1) is allowed to be reversible (k 0) [7]. The system is 3 ≥ fed by an input of the achiral substrate A, whereas the output consists of the inactive product P, and the excess of the enantiomers. We further assume that each enantiomer diffuses with the same diffusion constant D and include this feature in the master equation description of this process. 1 For the direct formation of chiral matter from achiral substrate, we would have had to include the steps A L(k ) and A D(k ). According to [8], k must be 0 0 0 → → kept sufficiently small. Otherwise, these chirally unspecific reactions would proceed too rapidly thus generating large amounts of racemic matter that swamps the am- plification process driven by the autocatalytic (1) and mutual inhibition steps (2). To take this important observation into account, we simply set k = 0 from the 0 outset. 2 2 Effective Bosonic Field Theory Our starting point for a systematic treatment of the above reaction scheme is an appropriate master equation. On a microscopic level, this comprises an exact description of the dynamics. From this equation, it is then straightfor- ward to derive an effective stochastic field theory. The salient steps are given in detail below. The chemical master equation for the kinetic scheme in Eqs.(1,2) is mapped to a “second-quantized” description following Doi [9]. Briefly, we introduce annihilation and creation operators ai and a†i for L, bi and b†i for D at each latticesite i,obeying thecommutation relations[ai,a†j] = δij,[bi,b†j] = δij.The vacuum state (corresponding to the configuration containing zero particles) satisfies ai|0i = bi|0i = 0. Furthermore, ai|mii = mi|mi−1i,a†i|mii = |mi+1i and a†iai|mii = mi|mii (and similarly for the b-sector). We then define the time-dependent state vector |Ψ(t)i = P({m},{n},t) (a†i)mi(b†i)ni|0i, (3) mX, n Yi { }{ } where P( m , n ,t) is the probability distribution to find m ,n particles of i i { } { } type L,D, respectively, at each site i, and satisfies the master equation. Then by means of Eq.(3), the master equation can be written as a “Schr¨odinger” equation ∂ Ψ(t) | i = H Ψ(t) , Ψ(t) = exp( Ht) Ψ(0) , (4) − ∂t | i ⇒ | i − | i where the lattice hamiltonian H is calculated to be D D H = l2 (a†i −a†j)(ai −aj)+ l2 (b†i −b†j)(bi −bj) X X (i,j) (i,j) −Xi nk1((a†i)2ai −a†iai)+k1((b†i)2bi −b†ibi) +k3(a†ia2i −(a†i)2a2i)+k3(b†ib2i −(b†i)2b2i) + k2(aibi −a†iaib†ibi)o, (5) and (i,j) denotes the sum over all lattice sites i (with lattice spacing l) and their nearest neighbors j in a d-dimensional space. Now take the continuum limit (l 0) [10] in order to obtain the path integral → representation of the underlying stochastic dynamics. Thus we can write the 3 time-evolution operator ¯ ¯ exp( Ht) = a a¯ b bexp( S[a,a¯,b,b]) (6) − Z D D D D − intermsofcontinuousfieldsa(x,t),a¯(x,t),b(x,t),¯b(x,t)withweightexp( S[a,a¯,b,¯b]). ¯ − Next, perform the shift a¯ = 1 +a , b = 1 + b on S. This yields the shifted ∗ ∗ action for any space dimension d, which is found to be: S[a,a ,b,b ]= ddx dt a ∂ a D 2a+k ab k a+k a2 ∗ ∗ ∗ t 2 1 3 Z Z n (cid:16) − ∇ − (cid:17) +b ∂ b D 2b+k ab k b+k b2 +a 2(k a2 k a) ∗ t 2 1 3 ∗ 3 1 (cid:16) − ∇ − (cid:17) − + b 2(k b2 k b)+k a b ab . (7) ∗ 3 1 2 ∗ ∗ − o The action S contains all the dynamics of the reaction scheme (1,2) including diffusion, and apart from taking the continuum limit, is exact. In particular, no assumptions regarding the form of the noise are required. At this stage, there are at least two independent ways to extract detailed information from S: by (i) employing field theoretic renormalization group techniques [11] or (ii) by exploiting the exact equivalence of this action to coupled Langevin equations. The former is suited for understanding chiral symmetry breaking as a problem in dynamic critical phenomena, the latter is ideal for understanding chiral amplification, evolution of enantiomeric excesses and net reactor yield as a problem in spatial pattern formation. 3 The Exact Langevin Equations For the final step we use the gaussian transformation [12] which allows us to integrateexactly over theconjugatefields a ,b ,appearing inthepathintegral ∗ ∗ a a b b e S[a,a∗,b,b∗]. This final step yields a product of delta-functional ∗ ∗ − D D D D Rconstraintswhichinturn,immediatelyimplythefollowingsetofexact coupled stochastic partial differential equations: ∂ a=D 2a+k a k ab k a2 +η (x,t), (8) t 1 2 3 a ∇ − − ∂ b=D 2b+k b k ab k b2 +η (x,t), (9) t 1 2 3 b ∇ − − where theintrinsic multiplicative reactionnoise iscompletely and exactly spec- ified as follows, namely η = η = 0 and a b h i h i η (x,t)η (x,t) =(k k a(x,t))a(x,t)δd(x x)δ(t t), (10) a a ′ ′ 1 3 ′ ′ h i − − − 4 η (x,t)η (x,t) =(k k b(x,t))b(x,t)δd(x x)δ(t t), (11) b b ′ ′ 1 3 ′ ′ h i − − − 1 η (x,t)η (x,t) = k a(x,t)b(x,t)δd(x x)δ(t t). (12) a b ′ ′ 2 ′ ′ h i −2 − − Intheabsence ofnoise(mean fieldlimit) thefields a(x,t),b(x,t) correspond to the coarse-grained local densities of the L,D particles. With noise, these fields are generally complex, and do not represent the physical densities. However, the averages a(x,t) , b(x,t) are real and do correspond to particle densities h i h i [13]. Forconvenience, we cast this system interms of nondimensional couplings and fields. Introducing L and T for spatial and temporal dimensions, dimensional analysis implies that [a] = [b] = 1/Ld, [D] = L2/T, [k ] = T 1, [k ] = [k ] = 1 − 2 3 Ld/T, [η] = 1/TLd. So, define dimensionless fields a˜ = (k /k )a, ˜b = (k /k )b, 2 1 2 1 time τ = k t, coordinates xˆ = k /Dx and dimensionless noises η˜ = 1 j 1 j a q (k /k2)η , η˜ = (k /k2)η . Then the equations Eqs.(8,9) can be written as 2 1 a b 2 1 b ∂ a˜= ˆ2a˜+a˜ a˜˜b ga˜2 +η˜ , (13) τ a ∇ − − ∂ ˜b= ˆ2˜b+˜b a˜˜b g˜b2 +η˜ , (14) τ b ∇ − − and the correlation matrix of the dimensionless noises is given by a˜(1 ga˜) 1a˜˜b  − −2  B = ǫ , (15)        1a˜˜b ˜b(1 g˜b)  −2 −    d/2 where g = k3 and ǫ = k2 k1 is the spatial dimension-dependent noise k2 k1(cid:16)D(cid:17) amplitude. Note that in d = 2 this reduces to ǫ = k2, indicating that the noise D is controlled via a competition between the dimerization reaction and spatial diffusion. We apply the Cholesky decomposition B = MMT to extract the square root of the matrix in the form of a lower triangular matrix M [14]. We will use ~ this decomposition to relate the noise to a new real white Gaussian noise ξ = (ξ ,ξ ), ξ (xˆ,τ)ξ (xˆ ,τ ) = δ δ(xˆ xˆ )δ(τ τ ), such that (η˜ ,η˜ ) =~η˜= Mξ~ 1 2 i j ′ ′ ij ′ ′ a b h i − − (and thus ~η˜T = ξ~TMT). Notice that by doing so the condition < ~η˜~η˜′T >=< Mξ~ξ~TMT >= M < ξ~ξ~T > MT = MMT = B is satisfied. This decomposition ′ ′ will allow us to separate the real and imaginary parts of the noise, a useful feature to have for setting up a numerical integration of the (complex-valued) 5 stochastic reaction-diffusion system [14]. We verify that √B 0 11   M = (16)      B21 B B221   √B11 r 22 − B11  √a˜(1 ga˜)1/2 0  −  = √ǫ (17)       −21˜b√a˜(1−ga˜)−1/2 (1−ga˜)−1/2q˜bq(1−g˜b)(1−ga˜)− 41a˜˜b yields the desired decomposition of B. 4 Numerical results The non-trivial mean-field solutions of Eqs. (13,14), in the absence of fluctu- ations (deterministic case, ǫ = 0), have been worked out in [6]: (I) if g < 1: ˜b = 0, a˜ = 1, homochiral absorbing state solution, • g (II) if g < 1: a˜ = 0, ˜b = 1, homochiral absorbing state solution, • g (III) if g > 1: a˜ =˜b = 1 racemic active state solution, • 1+g where g = k3. In solutions (I) and (II), one of the two enantiomers is amplified k2 to its maximal value and the other one disappears, this corresponds to pure chiral amplification. On the contrary, in solution (III) the two enantiomers coexist with the same concentration, this corresponds to the racemic solution. This scenario isonlyvalidinthemean-field limit.Ifthefullproblemistreated, then the noise due to spatial fluctuations and the diffusion of reactants must properly be taken into account. We will solve the full stochastic two-dimensional version of Eqs. (13,14) nu- merically, using reflecting boundary conditions and a finite difference scheme with ∆τ = 0.05, ∆xˆ = ∆yˆ= 0.51, and a grid of size L L = 154 154 2. × × We first considered the evolution of a system with g = 0.09 = 0.643 < 1. Two 0.14 representative results are shown in Figs. 1 and 2. The initial condition was ˜ set to a homogeneous, racemic situation (a˜ = 0.5,b = 0.5), and the noise 0 0 2 Inournumericalstudies,thesystem sizedidnothave any effect on thestationary solutions, only on the time required to reach equilibrium. 6 strength ǫ = k2 = 10 3 was fixed. A sequence of simulations were indepen- D − dently executed: every time the same mathematical problem was solved, the only difference being the random sequence of applied noises, and the result- ing evolution of each realization. In the two sets of rows, we plot the time evolution of the spatial distribution of the real part of a˜(xˆ,yˆ,τ) (upper row) ˜ and b(xˆ,yˆ,τ) (lower row). From left to right, the spatial distributions are shown at τ = 0,83,166,1111,8333 and 140000, respectively. When displayed in color, red represents a concentration of 0, purple represents a concentration of roughly 0.5, dark blue roughly 0.6 and green is the maximal concentration, 1.555 (=1). In both simulations, after the transient time (τ > 100, i.e. after g the first two frames of the sequences in Figs. 1 and 2), we found that even if the initial state was homogeneous and racemic (see purple squares in Fig- ures 1 and 2), the spatial fluctuations rapidly drive the system to a pattern with marked spatial compartamentalization, where, locally, the system is pure in one of the enantiomers 3. This means that we can find coexistence of lo- cal environments with the homochiral solutions (I) or (II). In the boundaries between regions of different chiralities there are racemic fronts of type (III), which in the mean-field limit would have only been expected for g > 1. Note the fully complementary nature between the spatial patterns formed of the ˜ a˜ and b-sequences during their evolution. Once the system has reached this state, the different regions compete and merge and the system evolves to one of the following final states: A) A homochiral, homogenous solution where one of the enantiomers has • entirely vanished, see Fig. 1. B) A mixed state with strong spatial segregation of two immiscible phases, • where locally in each phase, one of the homochiral solutions is realized, separated by planar racemic interface, see Fig. 2. Notice that when the solution is homochiral (solutions (I) and (II)) the noise vanishes since the correlation matrix, Eq. 15, vanishes as well. Thus this final state is absorbing. The racemic fronts are fluctuating and active. We verified that the complex solutions of Eqs.(13,14), once averaged over the fluctuations, do indeed yield real results, in accord with the theoretical expectations [13]. Next, we evaluate the time evolution of the spatially av- 3 Thus from the third square onward from left to right in both Figures 1 and 2, we see the formation of red and green zones. Top row Figure 1, red and green signal the absence of and maximal concentration of a˜, respectively. Bottom row Figure 1: red and green signal absence of and maximal concentration of ˜b, respectively. Final state in this sequence is zero concentration a˜ (red square top row) and maximal concentration for ˜b (green square bottom row). In Figure 2, the final state (right- hand most squares in top and bottom rows) represents half the domain filled by maximal a˜ and the complementary half square filled by maximal ˜b concentrations, respectively, each half domain separated by a racemic boundary line. 7 a˜ ˜ b Fig. 1. Simulation 1. Starting from a homogeneous racemic condition, the local random fluctuations induced by imperfect mixing drive the system to a homochiral solution. Time runs from left to right. See the text for details. a˜ ˜ b Fig. 2. Simulation 2. The system is solved again with the same conditions as in Fig. 1. The new random fluctuations induced by imperfect mixing lead the system to a different final state with spatial segregation. Time runs from left to right. eraged (over xˆ,yˆ) densities for the two cases shown above and for two other simulations representative of the inverse situation (i.e., where the system is dominated by the other enantiomer). We represent in Fig. 3 the time evolu- tion of the averaged enantiomeric excess ee(τ) = <a˜> <˜b>. When the solution <a˜>−+<˜b> is spatially homogeneous with a˜ at its maximal value, then ee = +1, for the ˜ inverse situation where b is amplified maximally, then ee = 1, whereas for − theintermediate situation withspatial compartmentalization 0 < ee < 1if the ˜ region with a˜ is greater than the one with b, and 1 < ee < 0 for the opposite − situation. All these outcomes were found with roughly equal probability over the total of the 20 numerical simulations performed. As for the evolution of a system with g = 0.14 = 1.555 > 1 the system con- 0.09 ˜ verges to the homogeneous solution a˜ = b = 1/(1+g) with small fluctuations about this value, i.e. the racemic state (III). All the results shown appear for a wide range of fluctuation intensities (ǫ = 10 2 10 6). − − − 8 1 Simulation 4 0.5 Simulation 2 ) τ ( 0 e e Simulation 3 -0.5 Simulation 1 -1 0 50000 100000 150000 τ Fig. 3. Time evolution of the spatially averaged enantiomeric excess ee(τ) = <a˜> <˜b> of four representative simulations of the samestudycase, starting <a˜>+−<˜b> fromracemicconditions(notresolvableinthefigure).Duetotheinherentsymmetry of the reaction, either enantiomer a˜ or ˜b can dominate. The final result varies from one simulation to another. One can have total, homogeneous, chiral amplification (Simulations 1 and4) or local chiral amplification with aracemic front (Simulations 2 and 3). 5 Conclusions Internal particle density fluctuations automatically and inevitably lead to chi- ral symmetry breaking in an extension [7,6] of the Frank model. The inclusion of this noise is carried out through a first-principles technique which allows us to map the chemical master equation to the corresponding coarse-grained Langevin equations. This multiplicative noise is rigorously determined by this procedureandisnotanad-hocadditiontothemeanfieldequations. Amarked advantage of the fully stochastic model is that no initial chiral seed nor ex- ternal chiral field [15] are needed: full homochirality is achieved starting from a strictly racemic initial condition. In the limit as ǫ 0, we would likely → be sensitive to round-off errors, which also seem to induce chiral symmetry breaking, as reported in [16]. In two dimensions, there is a second class of solutions (see Simulation 2 in Fig. 2) for which the final state consists of two enantiomeric pure regions separated by a racemic interface or boundary. Re- garding the solutions displayed in Figs 1 and 2, it is amusing to point out that Frank, in conceiving of what may be expected to happen in imperfectly mixed systems, predicted the formationof separate “colonies” of thetwo kinds of enantiomers bounded by racemic surfaces, with chiral homogeneity main- tained within each such colony. He also stated that any initial curvature in these boundaries makes the boundary move towards the side from which it is concave, a feature which we have also observed in our simulations. The con- 9 sequence of this is that any enclave of the one species will eventually shrink, and the other survives; this is what we observe in simulations of type 1 (see Fig 1). An exception to this rule occurs “when the two seas are connected by a strait. A different species could survive in each sea, with the (racemic) boundary running stably from cape to cape” [4]. This is indeed what happens in simulations of type 2 (see Fig.2) An entirely different stochastic version of the Frank model was presented in [17]. Chiral amplification is also achieved starting from a racemic initial state. There the noise is incorporated in an ad-hoc fashion by assuming fluctuating kinetic constants k . i We thank Dr. Vladik Avetisov for his interest in this work, as well as Prof. Albert Moyano, and Prof. Josep M. Rib´o for useful discussions and comments. M.-P. Z. is supported by the INTA and the research of D.H. is supported in part by funds from CSIC and INTA. References [1] V. Avetisov and V. Goldanskii, Proc. Natl. Acad. Sci. USA 93 (1996)11435. [2] M. Avalos, R. Babiano, P. Cintas, J.L. Jim´enez and J.C. Palacios, Tetrahedron Asymm. 11 (2000) 2845. [3] D.K. Kondepudi and K. Asakura, Acc. Chem. Res. 34 (2001) 946. [4] F.C. Frank, Biochim. et Biophys. Acta 11 (1953) 459. [5] I.R. Epstein, Nature 374 (1995) 321. [6] I. Gutman, D. Todorovi´c and M. Vuˇckovi´c, Chem. Phys. Lett. 216 (1993) 447. [7] S.F. Mason, Nature 314 (1985) 400. [8] T. Buhse, J. Mex. Chem. Soc. 49(4) (2005) 328. [9] M. Doi, J. Phys. A: Math. Gen. 9 (1976) 1479. [10] L. Peliti, J. Physique 46 (1985) 1469. [11] For a recent review of these methods, see e.g., U.C. Ta¨uber, M. Howard and B.P. Volmayr-Lee, J. Phys. A: Math. Gen. 38 (2005) R79. [12] D.J. Amit, Field Theory, the Renormalization Group,and Critical Phenomena, McGraw-Hill, New York, 1978. [13] M.J. Howard and U.C. Ta¨uber, J. Phys. A:Math. Gen. 30 (1997) 7721. [14] D. Hochberg, M.-P. Zorzano and F, Mor´an, Phys. Rev. E 73 (2006) 066109. 10