The origin of the central dogma through conflicting multilevel selection

The central dogma of molecular biology rests on two kinds of asymmetry between genomes and enzymes: informatic asymmetry, where information flows from genomes to enzymes but not from enzymes to genomes; and catalytic asymmetry, where enzymes provide chemical catalysis but genomes do not. How did these asymmetries originate? Here, we show that these asymmetries can spontaneously arise from conflict between selection at the molecular level and selection at the cellular level. We developed a model consisting of a population of protocells, each containing a population of replicating catalytic molecules. The molecules are assumed to face a trade-off between serving as catalysts and serving as templates. This trade-off causes conflicting multilevel selection: serving as catalysts is favoured by selection between protocells, whereas serving as templates is favoured by selection between molecules within protocells. This conflict induces informatic and catalytic symmetry breaking, whereby the molecules differentiate into genomes and enzymes, establishing the central dogma. We show mathematically that the symmetry breaking is caused by a positive feedback between Fisher’s reproductive values and the relative impact of selection at different levels. This feedback induces a division of labour between genomes and enzymes, provided variation at the molecular level is sufficiently large relative to variation at the cellular level, a condition that is expected to hinder the evolution of altruism. Taken together, our results suggest that the central dogma is a logical consequence of conflicting multilevel selection.


Alternative agent-based models in which the mutation of k c pt is modelled differently
In this section, we describe alternative models for the mutation of k c pt . In the agentbased model described in the main text, the mutation of k c pt is modelled as unbiased random walks in a half-open interval p´8, k max q with a reflecting boundary at k c pt " k max . To ascertain that this specific model of mutation does not critically affect results, we additionally examined two alternative models of mutation. The first alternative model is nearly the same as the above, except that the reflecting boundary condition is set at k c pt " 0. In the second alternative model, each k c pt value is mutated by multiplying expp q, where is a number randomly drawn from a uniform distribution on the interval p´δ mut , δ mut q, with a reflecting boundary at k c pt " k max . Both models of mutation produce essentially the same result as described in the main text (Figs. S3 and S4), indicating that the results do not depend on the specific models of mutation.

The derivation of equation (1)
In this section, we describe the derivation of equations (1) that is outlined in Methods.
To derive equations (1), we simplified the agent-based model in two ways. First, we assumed that k c pt is independent of p and t. Under this assumption, a catalyst does not distinguish the replicator types of templates (i.e., k c pt " k c pt 1 for t ‰ t 1 ) and products (i.e., k c pt " k c p 1 t for p ‰ p 1 ). This assumption excludes the possibility of numerical symmetry breaking, but still allows catalytic and informatic symmetry breaking as described in the main text (see Results).
Second, we abstracted away chemical reactions by defining ω t ij as the probability that replicator j of type t in protocell i is replicated or transcribed per unit time. Let n t ij pτ q be the population size of this replicator at time τ . Then, the dynamics of n t ij pτ q can be mathematically described as The fitness of the replicator can be defined as the dominant eigenvalue λ ij of the 2ˆ2 matrix on the right-hand side of equation (S1). The equilibrium frequencies of P and Q are given by the right eigenvector v ij associated with λ ij . Fisher's reproductive values of P and Q are given by the corresponding left eigenvector u ij . These eigenvalue and eigenvectors are calculated as follows: Based on the above simplification, we now derive equations (1). For concreteness, we focus on the evolution of the average catalytic activity of P (denoted byk P in the main text). However, the same method of derivation is applicable to that of Q if P and Q are swapped.
Let κ P ij be the catalytic activity of replicator j of type P in protocell i (we use κ instead of k to distinguish κ P ij from k P pt ). Price's equation [1,2] states that where xx ij y, xx˜ijy, and E˜irxs are x averaged over the indices marked with tildes, σ 2 i rx, ys is the covariance between x and y over protocells, and σ 2 ij rx, ys is the covariance between x and y over the replicators in protocell i (one replicator is always counted as one sample in calculating all moments). Below, we show that equation (S3) is approximated by equations (1) up to the second moments of xκ P ij y and κ P ij , namely, σ 2 i rxκ P ij y, xκ P ij ys and E˜irσ 2 ij rκ P ij , κ P ij ss. To approximate the first term on the right-hand side of equation (S3), we assume that xλ ij y is a function of xκ P ij y and xκ Q ij y that can be expanded as a Taylor series around xκ P ij y and xκ Q ij y. Substituting this series into σ 2 i rxλ ij y, xκ P ij ys, we obtain where Opσ 3 i q consists of terms involving the third or higher (mixed) central moments of xκ P ij y and xκ Q ij y over protocells [3]. To approximate the second term on the right-hand side of equation (S3), we likewise assume that λ ij is a function of κ P ij and κ Q ij that can be expanded as a Taylor series around xκ P ij y and xκ Q ij y. Substituting this series into σ 2 ij rλ ij , κ P ij s, we obtain where Opσ 3 ij q consists of terms involving the third or higher (mixed) central moments of κ P ij and κ Q ij over the replicators in protocell i [3]. Applying E˜i to both sides of the above equation and assuming that Bλ ij {Bκ c ij is independent of σ 2 ij rκ P ij , κ c ij s, we obtain Substituting equations (S4) and (S5) into equation (S3), we obtain ∆xκ P ij y " where O 1 " Opσ 3 i q`E˜irOpσ 3 ij qs. Next, we assume that covariances σ 2 i rxκ P ij y, xκ Q ij ys and E˜i because the mutation of κ P ij and that of κ Q ij are uncorrelated in the simulation model (this assumption is alternatively justified in Supplementary Material Text 1.6). Under this assumption, equation (S6) is transformed into Moreover, it can be shown that Using the above equation, we can transform equation (S8) into ij qs`E˜irOpσ 2 ij qsE˜irOpσ 2 ij qs. We adopt the following notation: xω t ij y xλ˜ijy , σ 2 cel " σ 2 i rxκ P ij y, xκ P ij ys, σ 2 mol " E˜irσ 2 ij rκ P ij , κ P ij ss, whereω t is the normalised average reproductive value of type-t replicators, σ 2 cel , σ 2 mol , and k P are the simplification of the notation, γ P P is an average decrease in the replication rate of a type-P replicator due to an increase in its own catalytic activity, and β t P is an increase in the average replication rate of type-t replicators in a protocell due to an increase in the average catalytic activity of type-P replicators in that protocell.
We assume that V is so large that xκ P ij y and κ P ij can be regarded as mathematically independent of each other, provided i and j are fixed (if i and j are varied, xκ P ij y and κ P ij may be statistically correlated). Under this assumption, increasing κ P ij does not increase xκ P ij y, so that γ P P reflects only the cost of providing catalysis at the molecular level. Likewise, increasing xκ P ij y does not increase κ P ij , so that β t P reflects only the benefit of receiving catalysis at the cellular level. Moreover, the independence of xκ P ij y from κ P ij implies that Bω Q ij {Bκ P ij " 0, which permits the following interpretation: if a replicator of type P provides more catalysis, its transcripts, which is of type Q, pay no extra cost (i.e., γ Q P " 0). Using the above notation and the fact that where O 2 is omitted. Equation (S10) is identical to equations (1). Finally, to derive the equation for ∆k Q (i.e., ∆xκ Q ij y), we swap P and Q in the above derivation. Moreover, we assume that ij rκ P ij , κ P ij ss because no difference is a priori assumed between P and Q.

The mathematical analysis of numerical symmetry breaking
In this section, we show that numerical symmetry breaking occurs because while it is neither favoured nor disfavoured by molecular-level selection, it is favoured by cellularlevel selection if catalytic and informatic symmetry breaking has occurred. To this end, we will again simplify the agent-based model into mathematical equations in a mannar analogous to that used to derive equations (1). Before describing the mathematical analysis, we first need to note that the proximateas opposed to ultimate-cause of numerical symmetry breaking is the self-replication of catalysts (i.e., k c cc ą 0, where c is the replicator type of catalysts) in the absence of the reverse transcription of catalysts (i.e., k c tc " 0, where t is the replicator type of templates). This fact can be inferred from the following two results. First, when catalytic, informatic, and numerical symmetry breaking occurs, the replication and transcription of templates are catalysed at about the same rate, i.e., k c tt « k c ct (Fig. 2b). Therefore, the replication and transcription of templates cannot cause numerical asymmetry. Second, when catalytic and informatic symmetry breaking occurs without numerical symmetry breaking, the selfreplication of catalysts is absent (Fig. S5). Taken together, these results indicate that the proximate cause of numerical symmetry breaking is the self-replication of catalysts in the absence of the reverse transcription of catalysts. Therefore, to understand why numerical symmetry breaking occurs, we need to understand why the self-replication of catalysts evolves if catalytic and informatic symmetry breaking has occurred.
To address the above question, we assume that replicators have already undergone catalytic and informatic symmetry breaking and consider how the fitness of those replicators depends on the self-replication of catalysts. The population dynamics of replicators with catalytic and informatic asymmetry can be described as follows. Let n t ij pτ q be the population size of replicator j of type t in protocell i at time τ . Let catalysts and templates be P and Q, respectively. Then, the dynamics of n t ij pτ q is mathematically described as follows: where w PP ij is the self-replication probability of catalysts, and ω Q ij is the replication and transcription probabilities of templates, which are assumed to be identical to each other. The fitness of replicators can be defined as the dominant eigenvalue (denoted by λ ij ) of the 2ˆ2 matrix on the right-hand side of equation (S11): (S12) The associated right eigenvector, which determines the stationary frequencies of P and Q, is (S13) Equation (S13) shows that we must assume ω Q ij ą w PP ij in order for P and Q to coexist. Equation (S13) also shows that the frequency of catalysts (i.e., 1{p2´w PP ij {ω Q ij q) increases with the self-replication of catalysts (i.e., w PP ij ), as stated in the beginning of this section. We first examine whether the self-replication of catalysts is favoured by molecularlevel selection. To this end, we consider how the fitness of replicators (i.e., λ ij ) depends on the self-replication of catalysts (i.e., w PP ij ). According to equation (S12), λ ij does not directly depend on w PP ij . However, λ ij can indirectly depend on w PP ij because λ ij increases with the frequency of catalysts in a protocell (i.e., E ij r1{p2´w PP ij {ω Q ij qs). This frequency increases with w PP ij if V is so small that a particular replicator can influence the frequency of catalysts in the protocell. However, if λ ij increases with w PP ij , the average fitness of replicators in the protocell (i.e., xλ ij y) must also increase. Therefore, we need to consider the relative fitness (i.e., λ ij {xλ ij y). The relative fitness is independent of w PP ij because catalysis is equally shared among templates within a protocell. Therefore, the self-replication of catalysts is neither favoured not disfavoured by molecular-level selection. We next examine whether the self-replication of catalysts is favoured by cellular-level selection. To this end, we consider how the fitness of a protocell depends on the average self-replication of catalysts in that protocell (i.e., xw PP ij y). The fitness of a protocell can be defined as the average fitness of the replicators in that protocell (i.e., xλ ij y). According to equation (S12), xλ ij y does not directly depend xw PP ij y. However, xλ ij y indirectly depends on xw PP ij y because xλ ij y increases with the frequency of catalysts in a protocell (i.e., E ij r1{p2ẃ PP ij {ω Q ij qs). This frequency increases with xw PP ij y, so that xλ ij y must also increase with xw PP ij y. Therefore, the self-replication of catalysts is favoured by cellular-level selection. Taken together, the above considerations indicate that the self-replication of catalysts is neutral with respect to molecular-level selection, but advantageous with respect to cellular-level selection. Therefore, numerical symmetry breaking results from the maximisation of fitness at the cellular level in the presence of catalytic and informatic asymmetry.
Finally, we mention an important consequence of numerical symmetry breaking. Numerical symmetry breaking causes a bottleneck effect on the population of replicators within a protocell. This bottleneck effect increases among-cell variance relative to withincell variance (i.e., σ 2 cel {σ 2 mol ); therefore, it has a stabilising effect on protocells [4,5]. In this regard, numerical symmetry breaking can be compared to life-cycle bottlenecks displayed by multicellular organisms and eusocial colonies (i.e., an organism or colony develops from only one or a few propagules), which are considered to reduce within-group conflict [6][7][8].

The hierarchical Wright-Fisher model
In this section, we describe a model that stochastically simulates the population dynamics described by equations (1), in which σ 2 mol and σ 2 cel are treated as dynamic variables dependent on m and V .
The simplifications involved in the derivation of equations (1), while illuminating, make the comparison between equations (1) and the agent-based model indirect. Specifically, equations (1) cannot be compared with the agent-based model in terms of the same parameters, because the equations treat σ 2 mol and σ 2 cel as parameters, which are actually dynamic variables dependent on m and V in the agent-based model. To fill this gap, we constructed a model that stochastically simulates the population dynamics described by equations (1) and treats σ 2 mol and σ 2 cel as dynamic variables dependent on m and V . This model is formulated as a hierarchical Wright-Fisher process. Replicators are partitioned into a number of groups (hereafter, protocells). Each replicator is individually assigned replicator type c P tP, Qu and two k c values. The fitness of a replicator is calculated according to equation (S14). In each generation, replicators are replicated or transcribed with probabilities proportional to ω c ij , so that the population dynamics matches equation (S1) on average. After the replication-transcription step, the protocells containing greater than V replicators are divided with their replicators randomly distributed between the two daughter cells. The protocells containing no replicators are discarded.
The mutation of k c is modelled as unbiased random walks with reflecting boundaries. That is, each k c value of a replicator is mutated with a probability m per replication or transcription by adding a number randomly drawn from a uniform distribution on the interval p´δ mut , δ mut q (δ mut " 0.1). The values of k c are bounded in r0, 1s with reflecting boundaries at both bounds.
To determine the condition for symmetry breaking, we simulated the above Wright-Fisher model for various values of V and m. The simulations show that symmetry breaking occurs only if V and m are sufficiently large (Fig. S8), a result that is consistent with the outcomes of the original agent-based model (Fig. 2). Given that the Wright-Fisher model involves many of the simplifications involved in equations (1), the above consistency supports the validity of the symmetry breaking mechanism described by equations (1).

The phase-plane analysis
In this section, we describe the phase-plane analysis outlined in Methods.
To perform the phase-plane analysis depicted in Fig. 3, we adapted equations (1) by defining ω t ij as a specific function of κ t ij (see the previous section for the meaning of ω t ij and κ t ij ). The following definition was employed: where the factor e xκ P ij y`xκ Q ij y represents the cellular-level benefit of catalysis provided by the replicators in protocell i, the numerator e´s κ t ij represents the molecular-level cost of catalysis provided by the focal replicator, the denominator 1{pxe´s κ P ij y`xe´s κ Q ij yq normalises the cost, and s is the cost-benefit ratio. The above definition of ω t ij was chosen to satisfy the requirement that a replicator faces the trade-off between providing catalysis and serving as a template, so that γ t t and β t c are positive; for example, if the cost γ t t were negative, it would actually be a benefit, so that there would be no trade-off. This requirement is satisfied if Bω t ij {Bκ t ij ă 0 and Bxω t ij y{Bxκ c ij y ą 0 for c " t and c ‰ t. Apart from this requirement, the definition was arbitrarily chosen for simplicity.
Under the definition of ω t ij in equation (S14), we obtain equations describing the evolution of xκ c ij y (denoted ask c in the main text) as follows. Since the evolution of xκ c ij y is described by equation (S6), we substitute equation (S14) into equation (S6). For this substitution, we need to calculate the derivatives of fitness. According to equation (S2), the fitness of a replicator is λ Moreover, the average fitness of replicators in a protocell is xλ ij y " e xκ P ij y`xκ Q ij y , so We substitute these derivatives into equation (S6) and use the fact that where c 1 ‰ c, ρ cel is the correlation coefficient between xκ P ij y and xκ Q ij y (i.e., ρ cel " σ 2 i rxκ P ij y, xκ Q ij ys{σ 2 cel ), and ρ mol is the average correlation coefficient between κ P ij and κ Q ij (i.e., ρ mol " E˜irσ 2 ij rκ P ij , κ Q ij ss{σ 2 mol ). To derive equation (S15), we have assumed that the variances of xκ c ij y and κ c ij are independent of c; i.e., σ 2 cel " σ 2 i rxκ c ij y, xκ c ij ys and σ 2 mol " E˜irσ 2 ij rκ c ij , κ c ij ss for c " P and c " Q. Equation (S15) can be expressed in a compact form as follows: where ∇ is a nabla operator (i.e., ∇ " rB{Bxκ P ij y, B{Bxκ Q ij ys T , where T denotes transpose), σ 2 tot " σ 2 mol`σ 2 cel , R " σ 2 cel {pσ 2 cel`σ 2 mol q, B " p1`ρ cel qpκ P ij`κ Q ij q, and C " pρ mol1 q lnpe´s κ P ij`e´s κ Q ij q`ρ mol spκ P ij`κ Q ij q. R can be interpreted as the regression coefficient of xκ c ij y on κ c ij [9] and, therefore, the coefficient of genetic relatedness [10]. The potential function RB´p1´RqC can then be interpreted as inclusive fitness.
Next, we set ρ mol " 0 and ρ cel " 0 in equations (S15) and let xκ c ij y be denoted byk c , obtaining where c 1 ‰ c. Comparing equations (S16) and (S10), we infer that which are identical to equations (5). Next, we omit O 2 in equation (S16) and replace ∆ with time derivative d{dτ , obtaining d dτk c " σ 2 cel´s e´sk c e´sk P`e´skQ σ 2 mol . (S17) Finally, to allow for the restriction on the range ofk c (i.e.,k c P r0, k max s), we multiply the right-hand side of equation (S17) with a function, denoted by Θpk c q, that is 1 if 0 ăk c ă k max and 0 ifk c " 0 ork c " k max . Multiplying Θpk c q with the right-hand side of equation (S17), we obtain The above equation was numerically integrated for s " 1 to obtain the phase-plane portrait depicted in Fig. 3. Equation (S15) allows for statistical correlations between κ P ij and κ Q ij at the molecular and cellular levels, i.e., ρ mol and ρ cel . Therefore, it can be used to examine the consequence of ignoring these correlations, which is one of the simplifications made in the derivation of equations (1) described in Supplementary Material Text 1.3. For this sake, we calculate the nullcline of ∆xκ c ij y. Setting ∆xκ c ij y " 0 in equation (S15) and omitting O 2 , we obtain xκ c 1 ij y « xκ c ij y`s´1 ln ρ mol sσ 2 mol´p 1`ρ cel qσ 2 cel p1`ρ cel qσ 2 cel´s σ 2 mol .
This equation shows that all parameters only appear in the intercept of the nullcline with the xκ c 1 ij y-axis. Let us denote this intercept as s´1 ln I. The way I qualitatively depends on σ 2 cel and sσ 2 mol is independent of ρ cel because´1 ă ρ cel ă 1. Therefore, we can assume that ρ cel " 0 without loss of generality. Next, to see how ρ mol influences I, we focus on the singularity of I by setting p1`ρ cel qσ 2 cel " sσ 2 mol` , where ą 0. Then, I " p1´ρ mol qsσ 2 mol { ´ρ mol . The way I qualitatively depends on sσ 2 mol { is independent of ρ mol because´1 ă ρ mol ă 1. Therefore, we can assume that ρ mol " 0 without loss of generality. Taken together, these calculations show that ignoring correlations between κ P ij and κ Q ij does not qualitatively affect the results, supporting the validity of equations (1). . The rate constants of complex formation were defined in such a way that coexistence between P and Q is neither favoured nor disfavoured by cellular-level selection. a, Phase diagram with a symmetric initial condition: k c pt " 1 for all combinations of c, p, and t, with both P and Q present at the beginning of each simulation. The symbols are the same as in Fig. 2a, except that the circles include cases in which one replicator type goes extinct. b, Dynamics of k c pt averaged over all replicators for m " 0.01 and V " 10000 in a. c, Phase diagram with an asymmetric initial condition:    Figure Figure S5: The absence of numerical symmetry breaking for small m and large V (see Supplementary Material Text 1.4). a, b, The dynamics of k c pt averaged over all replicators is shown for V " 10000 and m " 0.001 with two different initial conditions: a symmetric initial condition, where k c pt " 1 (a); an asymmetric initial condition, where k P PP " 0.95, k P PQ " 0.1, k P QP " 1, k P QQ " 1, and k Q pt " 0.1 (b). The self-replication of catalysts does not evolve for the symmetric initial condition, whereas it is maintained for the asymmetric initial condition (t min ą 1.2ˆ10 7 ). The dependence of the results on the initial conditions suggests the presence of bistability for V " 10000 and m " 0.001. c, d, The frequencies of P (catalysts) and Q (templates) are plotted as the functions of time. Numerical symmetry breaking does not occur for the symmetric initial condition, whereas it occurs for the asymmetric initial condition. The results indicate that numerical asymmetry depends on the self-replication of catalysts. e, f, Catalytic activities evolved for the symmetric initial condition (e) and for the asymmetric initial condition (f).  Figure S6: The effect of symmetry breaking on catalytic activities. The fraction of replicators 1´N S {N tot , which is a proxy for the overall catalytic activity of replicators, is shown as a function of m and V , where N S is the total number of S molecules in the system, and N tot " N P`NQ`NS . a, The original model, which allows symmetry breaking (i.e., Fig. 1). b, The model that excludes the possibility of symmetry breaking; specifically, it allows only one type of replicator (either P or Q). Black squares indicate extinction (i.e. N tot " N S ). t min ą 1.5ˆ10 7 .   Figure S7: Result for large m and V values. The dynamics of the agent-based model is shown for m " 0.1 and V " 10 5 , parameters outside the range examined in Fig. 2a and Fig. S6a. a, The dynamics of k c pt averaged over all replicators. b, The dynamics of the fraction of replicators 1´N S {N tot , where N tot and N S are the total numbers of particles and S molecules in the system, respectively. t min ą 1.8ˆ10 6 .  (1), treating σ 2 mol and σ 2 cel as variables dependent on m and V (see Supplementary Material Text 1.5). a, Phase diagram. Circles indicate no symmetry breaking (i.e.,k P «k Q « 1); diamonds, symmetry breaking (i.e.,k c « 0 andk c 1 « 1 for c ‰ c 1 ); stars, extinction (i.e.,k P «k Q « 0). s " 1 (cost-benefit ratio). The total number of replicators was 50V (approximately 130 protocells throughout simulations). The initial condition was k P " k Q " 1 for all replicators. Each simulation was run for 4ˆ10 5 generations. The extinction (i.e.,k P «k Q « 0) for large m and V is consistent with the phase-plane analysis of equations (1), which also shows extinction (i.e., k P «k Q « 0) for sufficiently large σ 2 mol {σ 2 cel (parameters outside the range examined in Fig. 3). The discrepancy between Fig. S8a and Fig. 2a is due the simplifying assumption made in equations (1) that k c pt is independent of p and t. If k c pt is allowed to depend on p and t, the flow of information from templates to catalysts can become completely unidirectional. Such unidirectional flow of information can resolve the dilemma between catalysing and templating and leads to the maintenance of high catalytic activities as described in Results. b, The dynamics ofk c for m " 0.001 and V " 1000 (no symmetry breaking). c, m " 0.01 and V " 1000 (symmetry breaking). d, m " 0.1 and V " 1000 (extinction).