Computing the Parameter Values for the Emergence of Homochirality in Complex Networks

The goal of our research is the development of algorithmic tools for the analysis of chemical reaction networks proposed as models of biological homochirality. We focus on two algorithmic problems: detecting whether or not a chemical mechanism admits mirror symmetry-breaking; and, given one of those networks as input, sampling the set of racemic steady states that can produce mirror symmetry-breaking. Algorithmic solutions to those two problems will allow us to compute the parameter values for the emergence of homochirality. We found a mathematical criterion for the occurrence of mirror symmetry-breaking. This criterion allows us to compute semialgebraic definitions of the sets of racemic steady states that produce homochirality. Although those semialgebraic definitions can be processed algorithmically, the algorithmic analysis of them becomes unfeasible in most cases, given the nonlinear character of those definitions. We use Clarke’s system of convex coordinates to linearize, as much as possible, those semialgebraic definitions. As a result of this work, we get an efficient algorithm that solves both algorithmic problems for networks containing only one enantiomeric pair and a heuristic algorithm that can be used in the general case, with two or more enantiomeric pairs.


Introduction
We study mathematical models of absolute asymmetric synthesis [1,2] that are used to explain the emergence of biological homochirality, which is, according to Frank [3], a natural property of life. Frank proposed, as early as 1953, a chemical mechanism to support his thesis. An important feature of this minimal model is that homochirality is produced by dynamic instability [4]. After Frank many other more sophisticated networks have been proposed (see for instance [5][6][7]), but all of them are based on the same idea: homochirality is the product of chemical instabilities. It is important to remark that the mathematical (stability) analysis [8] of those complex models is a hard piece of work. Fortunately, we found a particular symmetry in the Jacobian matrices of those models [9] that yields semialgebraic definitions of the instability regions where the symmetry-breaking can be observed. Most of those semialgebraic expressions are highly nonlinear and hard to sample. We used Clarke's Stoichiometric Network Analysis (SNA) [10] to reduce the complexity of those expressions.
All those ingredients were put together into an algorithmic tool, and software Listanalchem [11], that can be used to test models proposed to explain the origin of homochirality, and which can also help us to build new and better models. Further thermodynamics constraints must be taken into account, but that point will not be discussed here.
We begin with a mathematical presentation of the method. After that, we use the developed algorithm to analyze three representative models of biological homochirality taken from the available literature.

Network Models of Absolute Asymmetric Synthesis
We use the term absolute asymmetric synthesis to designate all the possible chemical mechanisms that, operating in achiral environments, can transform a racemic mixture into an enantiopure one. We suppose that all those chemical processes reduce to finite sets of chemical reactions acting on finite sets of chemical species. Thus, we suppose that all those processes can be suitably described by chemical reaction networks.

Definition 1.
A chemical reaction over the chemical species X 1 , ..., X n is an expression like c 1 X 1 + · · · + c n X n → d 1 X 1 + · · · + d n X n , where c 1 , ..., c n and d 1 , ..., d n are small integers (some of which could be equal to zero). The latter expression indicates that the mixture of c 1 units of X 1 , ..., and c n units of X n gives place to d 1 units of X 1 , ..., and d n units of X n . A chemical reaction network over the species {X 1 , ..., X n } is a set of chemical reactions, say the set {R 1 , ..., R r } , over this set of species. Notation 1. Given a chemical network Ω = {(X 1 , ..., X n ) , (R 1 , ..., R r )}, we use the expression c 1i X 1 + · · · + c ni X n → d 1i X 1 + · · · + d ni X n (2) to denote the reaction R i , and we use the symbol k i to denote its reaction rate constant. We use variables [X 1 ] , ..., [X n ] to denote the concentrations of the n chemical species.

Remark 1.
It is important to remark that the network Ω F was the first, and it is the most elementary model of absolute asymmetric synthesis proposed in the literature [3].
The dynamics of a network Ω reduces to the temporal evolution of the concentration variables [X 1 ] , ..., [X n ] . We suppose that those temporal evolutions are completely determined by the law of mass action [12], that is: we suppose that the temporal evolution of the variables [X 1 ] , ..., [X n ] is given by the system of Ordinary Differential Equations (ODE) ] c 1j · · · · · [X n ] c nj , i = 1, ..., n, where given j r the symbol k j denotes the reaction rate constant of R j . We get consequently that all those dynamics are deterministic: their evolution in time is entirely determined by their initial states.

Remark 2.
Let s be a (initial) state of network Ω, state s is a vector that encodes the values of all the parameters that participate in the dynamics of Ω. Thus, we have that s can be fully described by a (n + r)-tuple ([X 1 ] , ..., [X n ] , k 1 , ..., k r ) , (4) of nonnegative reals. We must notice that the entries k 1 , ..., k r could remain constant along the dynamics, but in despite this, we choose to include the reaction rate constants in our notion of state.

Definition 2.
A state ([X 1 ] , ..., [X n ] , k 1 , ..., k r ) is said to be a steady state, if and only if, it satisfies the following system of polynomial equations: . . .
Notice that the steady states are the mathematical equilibria of the system. However, if one thinks in chemical equilibrium, it could be possible to consider the states for which only the set of forward reaction rates vanish. Also, it could be possible to consider the existence of complex networks containing loops, where not all forward rates vanish identically at the steady state.
Let us consider the case of network Ω F when the reagent [A] is assumed constant, according to the pool chemical approximation, see [13] chapters 2 and 3. The dynamics of this network is governed by the ODE system hold. Observe that k 1 must be equal to k 2 ; otherwise, the network encodes a chiral environment. Analyzing the dynamics of Ω F corresponds to analyzing the set SS (Ω F ) constituted by all its steady states, as well as the dynamics that can be triggered by arbitrary small perturbations of those states.

The Goal
We want to contribute to the investigation on the emergence of biological homochirality by providing the interested researchers with mathematical and algorithmic tools that can be used in the analysis of any network model of absolute asymmetric synthesis. Definition 3. We say that a chemical reaction network Ω exhibits chiral amplification, if and only if, it has the ability of transforming negligible gaps between the concentrations of the chiral species into larger gaps.
We would like to characterize the set of chemical reaction networks that exhibit chiral amplification.
Let Ω = {(X 1 , ..., X N ) ; (R 1 , ..., R r )} be a chemical reaction network, and suppose that X 1 is the L-form of a chiral biomolecule, and that X 2 is the D-form of the same molecule. Suppose that the system got stuck at a steady state satisfying the equality [X 1 ] = [X 2 ]. There are physical mechanisms that may create small enantiomeric gaps (that could perturb this racemic steady state). Those perturbations trigger dynamics, and those latter dynamics could: 1. Evolve towards the original steady states, if, for instance, those states are stable. 2. Evolve towards different racemic states, if, for instance, those sets of racemic states are attractors of the dynamics. 3. Evolve towards states with large enantiomeric gaps (also called scalemic states).
We are interested in the latter case, and we say that in that case, the system undergoes homochiral dynamics. We want to characterize the steady states of Ω that can undergo homochiral dynamics. We introduce, below, a precise formulation of the algorithmic problem that we study (and solve) in this paper.

Pseudochiral Networks
Let Ω = (X, R) be a chemical reaction network, and suppose that it is a network model of absolute asymmetric synthesis. Then, the set X must be constituted by three disjoint sets of chemical species, the sets {L 1 , ..., L k };{D 1 , ..., D k } and {X 2k+1 , ..., X n }.
The set {L 1 , ..., L k } constitutes the L-side of Ω, the set {D 1 , ..., D k } constitutes the D-side and the set {X 2k+1 , ..., X n } is constituted by all the achiral species occurring in Ω. Moreover, given i ≤ k we have that (L i , D i ) is an enantiomeric pair, that is: L i is the L-form of a chiral molecule, while D i is the D-form of the same molecule.
The enantiomeric gap of s is equal to A state of Ω is said to be racemic, if and only if, its enantiomeric gap is equal to 0. Thus, the racemic condition for Ω corresponds to the following set of polynomial equalities Recall that we are interested in the racemic states of Ω that can produce homochiral dynamics. It is important to take into account the indiscernibility of enantiomeric pairs: it is known that the two species of an enantiomeric pair react with the same achiral species at the same reaction rates. The indiscernibility of enantiomeric pairs implies that any feasible network model of absolute asymmetric synthesis must be a pseudochiral network, as defined below.

Definition 5.
Let Ω = (X, R) be a network model of absolute asymmetric synthesis, and suppose that where {L 1 , ..., L k } is the L-side of Ω, and {D 1 , ..., D k } is its D-side. We say that Ω is a pseudochiral network, if and only if, it is indistinguishable from the network that is obtained when one switches the L-species and the D-species.
Consider the following example of a pseudochiral network. Melvin Calvin, whose work laid the foundations of chemical evolution, proposed an abstract model of biological homochirality [14]. Calvin's mechanism can be suitably described by a chemical reaction network that we denote with the symbol Ω C , and which is defined as follows: The constituent species of Ω C are the species L 1 , L 2 , D 1 and D 2 . The chemical reactions in Ω C are: 1. A pair of autocatalytic reactions as well as the reverse reactions 2. The racemization reactions 3. The enantiomeric conversion reactions as well as their reverse reactions Notice that the reactions in Ω C are organized in pairs: for each reaction involving the species of the L-side there is a dual reaction involving the species of the D-side. We have, for instance, that if the reaction L 1 + L 2 → 2L 2 occurs, then the dual reaction D 1 + D 2 → 2D 2 must also occur; otherwise, the network would distinguish between the species of the L-side and the species of the D-side. The existence of dual reactions is a consequence and is somewhat equivalent to the pseudochirality of network Ω C . Moreover, given a reaction R L and its dual R D , the reaction rate constants of those two reactions are the same. In the particular case of Calvin's mechanism, we have k 0 = k 1 , k 2 = k 3 , k 4 = k 5 , k 6 = k 7 , and k 8 = k 9 .
Notation 2. In the Calvin model, we use the symbol k 0 to denote the reaction rate constant of the dual pair of autocatalytic reactions, the symbol k 2 to denote the reaction rate constant of their reverse reactions, the symbol k 4 to denote the reaction rate constant of the dual pair of racemization reactions, the symbols k 6 to denote the reaction rate constants of the enantiomeric conversions and k 8 to denote the reaction rate constants of their reverse reactions.
We have that Ω C is a pseudochiral network of order 2. Does network Ω C exhibit homochiral dynamics? Calvin claimed that it is the case, but we think that the evidence provided by him is weak, and we would like to consider this question more carefully.

Remark 3.
Let Ω be a network, and suppose that Ω is not pseudochiral. Then, there is an asymmetry in Ω that distinguishes the L-side and the D-side. Thus, if Ω is not pseudochiral, it models a chemical mechanism that works on a chiral environment, and because of this it cannot be regarded as a network model of absolute asymmetric synthesis. We claim that: feasible network models of absolute asymmetric synthesis are pseudochiral networks.

The Algorithmic Problem
Suppose that one wants to introduce a new network model of absolute asymmetric synthesis. One must show that this new model exhibits homochiral dynamics, and one must also show that the model is sound from a thermodynamical point of view. We focus on the first task, which is solved if one exhibits racemic steady states which, after being perturbed, trigger homochiral dynamics. Therefore, we consider the following algorithmic problem: It has been argued that the asymmetric synthesis of chiral biomolecules was a prerequisite for the origin of life [3]. It is supposed that there are chemical processes, which took place in prebiotic earth, and which transformed the initial racemic mixtures, present in prebiotic earth, into the homochiral mixtures that preceded the origin of life. Frank introduced in his seminal work [3] an abstract chemical mechanism, which contains an enantiomeric pair, and which evolves towards enantiopure states independently of the initial state. The chemical mechanism introduced by Frank is well described by the chemical reaction network Ω F . This network was the first network/mathematical model of absolute asymmetric synthesis introduced in the literature. After that, many other network models have been proposed, and we know that some of those proposed models are defective models that cannot support homochiral dynamics (we prove that the Calvin model [14] does not exhibit homochiral dynamics, see below). Therefore, we ask: can we recognize and discard all those defective models? To answer the question, we need an algorithm able to: 1. Recognize the defective network models of absolute asymmetric synthesis that cannot exhibit homochiral dynamics. 2. Recognize the network models that are mathematically sound, and compute samples of their sets of racemic steady states that undergo homochiral dynamics.
Thus, we need an algorithm that solves problem CPVEH.

The MM-Condition
Let Ω be a pseudochiral network, and suppose that we want to show that Ω has racemic steady states that produce homochiral dynamics. How can those states be found? Most authors cope with the latter question using the tools of classical stability analysis: the states that produce homochiral dynamics are unstable (see below).

Remark 4.
From now on we use the symbol Ω to denote a pseudochiral network of order k. Moreover, we suppose that and we say that Ω is a network of size n.

Notation 3.
We use the symbol J Ω to denote the Jacobian of Ω. The Jacobian J Ω is a symbolic matrix whose entries are the partial derivatives [8] of the reaction rates [15]. The entries of J Ω are polynomials over the variables Thus, given a state s, one can evaluate J Ω at s and obtain a numerical matrix J Ω (s). We say that J Ω (s) is the Jacobian of Ω at (state) s.
Let us consider the pseudochiral network Ω C . The steady state equations for the L-species are equal to: In addition, if we assume the racemic condition, those equalities get equal to while the Jacobian matrix gets equal to Remark 5. The reader must observe the symmetrical structure of the above matrix, which we call the racemic Jacobian of Ω C . It happens that the racemic Jacobian of any pseudochiral network Ω reveals the same type of symmetries [9].
The Jacobian matrix J Ω encodes important information related to the dynamics of Ω: We say that a steady state s is unstable, if and only if, the spectrum of J Ω (s) satisfies the following: There exists λ that belongs to the spectrum of J Ω (s) and such that the inequality Re (λ) > 0 holds.
We say that a polynomial p (X) is unstable, if and only if, the roots of p (X) satisfy the same condition imposed on the spectrum of the unstable states. Notice that s is unstable, if and only if, the characteristic polynomial of J Ω (s) is unstable.
Recall that the spectrum of a square matrix is the set constituted by all its eigenvalues.
The unstable states are the states that can get dramatically transformed by the effect of negligible perturbations. We use the term symmetry-breaking states to designate the racemic steady states that can produce homochiral dynamics. Notice that those homochiral dynamics dramatically transform those racemic states. Thus, we get that any symmetry-breaking state is unstable. However, we must observe that there exist unstable states that do not produce homochiral dynamics. We must ask: which are the symmetry-breaking states of network Ω? Which is the mathematical condition that determines the symmetry-breaking status? We observed before that the symmetry-breaking status depends on the eigenvectors of J Ω (s) and not only on its eigenvalues, see [9,16,17].

Definition 7.
We say that a steady state s is symmetry-breaking, if and only if, matrix J Ω (s) satisfies: 1. There exists λ that belongs to the spectrum of J Ω (s), and such that Re (λ) > 0. 2. There exists an eigenvector of λ, say v, and there We get that the symmetry-breaking states of Ω are the unstable states that satisfy a further constraint. The additional constraint refers to the eigenvectors of J Ω (s). In principle, this additional constraint should make harder the search of those states. However, and surprisingly, it is not the case: if one exploits the symmetries of J Ω , the latter problem becomes easier than the former.

Notation 4.
Let Ω be a pseudochiral network of order k (recall that the order k is equal to the number of enantiomeric pairs, k = 2 for the Calvin model), we use the symbol A Ω to denote the submatrix of J Ω that is constituted by its first k rows and its first k columns. We use the symbol B Ω to denote the submatrix that is constituted by the first k rows of J Ω (s) and the columns k + 1, ..., 2k. For instance, if we consider Calvin network Ω C , we get that Theorem 1. (MM-condition) Let s be a racemic steady state of the pseudochiral network Ω, we have that s is symmetry-breaking, if and only if, the characteristic polynomial of A Ω (s) − B Ω (s) is unstable [9].
The above theorem allows us to analyze any pseudochiral network. Let us illustrate the latter claim with the analysis of the network Ω C .
The symmetry-breaking states of Ω C are the 9-tuples that satisfy the following constraints: 1. The racemic condition 2. The steady state condition 3. The non-negativity condition 4. The positivity condition k 0 , k 2 , k 4 , k 6 , k 8 > 0 (30) 5. The (symmetry-breaking) MM-condition asserting that the characteristic polynomial of matrix is unstable.
Notice that the first four constraints are polynomial inequalities (equalities) over the variables 6 , and k 8 . The last constraint can also be expressed in terms of polynomial inequalities. To this end one can use a suitable set of Hurwitz-Routh inequalities [18]. We have, for instance, that the 2 × 2 matrix A Ω C − B Ω C is unstable, if and only if, at least one of the following two inequalities is satisfied: It is easy to prove that there are no steady states satisfying at least one of the above inequalities.

Theorem 2. The Calvin model does not have symmetry-breaking states.
Proof of Theorem 19. First we consider the inequality Let us suppose that [L 2 ] = 0. From the steady state condition, Equation (28), we get that and hence the Inequality (32) cannot be satisfied.
On the other hand, if we suppose [L 2 ] = 0, we get that the symmetry-breaking states of Ω C must be solutions of the system which does not have solutions satisfying the constraints k 0 , k 6 , k 8 > 0 and [L 1 ] ≥ 0. Now, we consider the inequality We get from the steady state condition (28) that and we get consequently that Then, the Inequality (32) cannot be satisfied, and the theorem is proved.
The analysis of the Calvin network, as developed in the previous paragraphs, shows that it is possible to use the MM-condition to analyze small networks thoroughly. What can be done with larger networks? The MM-condition yields an algorithmic solution for the CPVEH problem. However, if the input network Ω (mechanism) is too large, its analysis could be intractable. Consider the problem: The MM-condition allows us to reduce the problem CPVEH to this latter problem. However, this reduction can be useless, given that problem SA is a hard, intractable problem. Thus, we must ask: can we efficiently solve problem CPVEH?
We claimed before that computing the symmetry-breaking states of Ω can be easier than computing its unstable states. The MM-condition strongly reduces the dimensionality of the problem. This mathematical criterion tells us that we must analyze a k × k symbolic matrix instead of the n × n symbolic matrix that we would have to analyze if we were to use classical stability analysis to compute, as in previous approaches, the unstable states of Ω (notice that n ≥ 2k). It seems that it is not possible to further reduce the dimensionality of the problem, and it forces us to consider other completely different reductions. We must observe, at this point, that there are two main sources of complexity related to the instances of problem SA: the dimension of the input matrices, and the degree of their polynomial entries. Thus, we ask: can we also reduce the inherent nonlinearity of the instances of CPVEH? In the following section we introduce some tools of Clarke's Stoichiometric Network Analysis (SNA, see reference [10]), which allow us to achieve the degree reduction we are looking for. We use SNA, and the degree reduction provided by it, to develop an algorithm that can be used in the analysis of pseudochiral networks.

Degree Reduction Using Stoichiometric Network Analysis
Clarke's SNA provides us with tools that can be used in the linear stability analysis of chiral networks, see [19] and the references therein. Here, we use SNA to develop a heuristic algorithm for the CPVEH problem.

A Crash Introduction to SNA
Let Θ be a chemical reaction network and let SS (Θ) be its set of steady states. SNA is based on a change of variables that linearize the definition of SS (Θ): this change of variables maps the latter set onto a polyhedral cone that we denote with the symbol C Θ .

Notation 5.
Along this section we use the symbol Θ to denote a chemical reaction network. Moreover, we suppose that and we suppose that R i is equal to c 1i X 1 + · · · + c ni X n → d 1i X 1 + · · · + d ni X n .
Definition 8. The stoichiometric matrix of Θ that we denote with the symbol S Θ , is a n × r matrix whose entries are called the stoichiometric coefficients of Θ. The stoichiometric coefficients of the network Θ are defined by: A second matrix related to S Θ is the matrix of orders of reaction denoted with the symbol R Θ , and which is the r × n matrix R Θ = c ji i≤n,j≤r .
(42) Definition 9. The velocity function of Θ is the function V Θ : R n → R s that is defined by Clarke observed that the dynamic equation of Θ can be written as and it implies that the function V Θ maps the set SS (Θ) onto the polyhedral cone C Θ that is defined by the linear constraints Thus, the function V Θ allows us to identify the steady states of Θ with the points of C Θ , and it happens that the points of C Θ can be suitably described by the system of convex coordinates provided by its extreme currents [10].

Remark 6.
A set of extreme currents for C Θ is just a minimal set of extreme rays that spans the polyhedral cone C Θ , see [20]. Notation 6. We use the symbol R s + to denote the set of s-dimensional nonnegative vectors, which is the set Definition 10. Let v 1 , ..., v s be a set of extreme currents for C Θ , and let s ∈ SS (Θ); the convex coordinates of s are given by the unique tuple (j 1 , ..., j s ) ∈ R s + that satisfies the equality From now on, we use the symbols j 1 , ..., j s to denote the convex coordinates of the cone C Θ that are determined by v 1 , ..., v s .
If one switches to the system of convex coordinates, then the Jacobian J Θ can be factorized as where ∆ is a n × n scaling matrix, and E Θ is the r × r diagonal matrix defined by

Remark 7.
The term E Θ [i, i] is a linear polynomial over the (convex) variables j 1 , ..., j s .
Clarke noticed that the scaling matrix ∆ has little influence on the stability properties of J Θ , and that one can focus the analysis on the matrix S Θ · E Θ · R Θ . Thus, according to Clarke's theory, the stability analysis of Θ reduces to the stability analysis of the latter matrix. We must ask: what do we gain with this reduction? The entries of S Θ · E Θ · R Θ are linear polynomials over the variables j 1 , ..., j s , and it means that we get an important degree reduction.
We use the symbol V Θ to denote the matrix S Θ · E Θ · R Θ . The SNA-based algorithmic analysis of chemical networks reduces to: 1. Compute a set of extreme currents for C Θ . 2. Compute the symbolic matrix V Θ . 3. Sample the set of parameter values that make the characteristic polynomial of matrix V Θ become unstable. 4. Given α, a s-dimensional vector that belongs to the sample computed in step 3, compute a sample of the set V −1 Θ (α 1 v 1 + · · · + α s v s ) . Here, we use the symbols v 1 , ..., v s to denote the set of extreme currents for C Θ as computed in step 1.

Remark 8.
The latter algorithmic routine can be effectively implemented, see [11].

An Illustrative and Trivial Example
Clarke's velocity function denoted by V Θ 0 is equal to the vector field k 1 The set SS (Θ 0 ) is the subset of R 2 + × R 3 + that is defined by the equation Remark 9. Notice that SS (Θ 0 ) is a 5-dimensional set with a quite complex structure.
The stoichiometric matrix of Θ 0 is the matrix and a set of extreme currents for C Θ 0 is given by {(1, 3, 0) , (0, 1, 1)} . This means that the 5-dimensional set SS (Θ 0 ) is mapped by R Θ 0 onto the infinite triangle spanned by the vectors (1, 3, 0) and (0, 1, 1) . Thus, we get that V Θ 0 (SS (Θ 0 )) has a pleasant conic structure. We can use the extreme currents of C Θ 0 (which are integer vectors located on the borders of C Θ 0 , for instance the vectors (1, 3, 0) and (0, 1, 1)) to define a system of convex coordinates for C Θ 0 . Notice that the points of this triangle are in bijective correspondence with the nonnegative linear combinations of the extreme currents (1, 3, 0) and (0, 1, 1). Thus, if we choose nonnegative values for j 1 , j 2 , we can be sure that is an element of C Θ 0 that represents a nonempty set of steady states. We can use this fact to sample the set SS (Θ 0 ). We proceed in the following way: 1. Pick a tuple of positive convex coordinates; for instance, j 1 = 1 and j 2 = 2. 2. Compute the corresponding convex combination of extreme rays; in our case, we compute 1 · (1, 3, 0) + 2 (0, 1, 1) equal to (1, 5, 2). 3. Compute the solution set of the nonlinear system 4. Pick an element, say (1, 1, 1, 5, 2), of the set computed in the previous step.
This means that given a chemical network Θ, we can use Clarke's velocity function to define a system of convex coordinates for the set SS (Θ). Moreover, we can use this system of coordinates to sample the set SS (Θ). We can also use the factorization of J Θ to analyze the stability of network Θ.
Additionally, it is important to illustrate how this procedure reduces the degree of the polynomials entries of the symbolic Jacobian J Θ 0 . In our illustrative and trivial example, the symbolic Jacobian is equal to We have that R Θ 0 is equal to and we have that E Θ 0 is equal to According to Clarke's factorization the matrix J Θ 0 is equal to and this factorization allows us to focus the stability analysis on the matrix which is equal to Notice that the entries of the matrix V Θ 0 are linear polynomials over the variables j 1 , j 2 , while the entries of J Θ 0 are 3-degree polynomials over the parameters [I] , [A] , k 1 , k 2 and k 3 . Thus, it is true that we gain an important degree reduction if we switch to the system of convex coordinates.

Using SNA in the Analysis of Pseudochiral Networks
We use the tools introduced in the previous paragraphs to develop an algorithm for the stability (symmetry-breaking) analysis of pseudochiral networks.
Let us consider the case of pseudochiral networks of order 1, which we call chiral networks. Thus, let Φ be a chiral network of order 1, and suppose that The MM-condition reduces to the inequality which is a linear inequality over the entries of J Φ . However, the linearity of this algebraic condition is just apparent, given that the terms J Φ [1, 1] and J Φ [1,2] are nonlinear polynomials over the concentration variables Notation 7. We use the symbol SB (Φ) to denote the set of symmetry-breaking states of Φ.
We have that the set SB (Φ) is a highly nonlinear set defined by the conditions: 1. The nonlinear steady state equations.

The nonlinear inequality
3. The non-negativity conditions 4. The positivity conditions 5. The racemic condition stating that the rate constants of the reactions that belong to the same dual pair are equal, and that the initial concentrations of the two enantiomeric species are also equal.
However, if we switch to Clarke's system of convex coordinates, we get that this set is defined by the linear conditions: • j 1 , ..., j s ≥ 0, where j 1 , ..., j s constitute the system of convex coordinates of network Φ.

Remark 10.
In the implementation of the algorithm, we forced condition five (the racemic condition) into the computation of the extreme currents, which are computed from the stoichiometric matrix; for this, we extended the stoichiometric matrix with rows that encode the equalities of the rate constants of the dual pairs; for example, given the stoichiometric matrix of the Calvin model the extended matrix, considering its dual pairs (Reactions (13), (14), (16), (17)), is equal to Forcing this condition into the computations of the extreme currents proved to reduce the number of extreme currents significantly.
The stoichiometric and extended stoichiometric matrices of the models analyzed with Listanalchem can be seen in the output of the computer program. The Supplementary Materials presents those matrices for the models studied in this work.
Sets defined by linear inequalities are easy to sample, see [21]. Then, if we switch to the system of convex coordinates, we will be able to efficiently sample the set of symmetry-breaking states of Φ. Thus, we can use the MM-condition and SNA to efficiently solve the CPVEH problem when it is restricted to pseudochiral networks of order 1.
What about higher orders? Suppose that Ω is a pseudochiral network of order k, suppose that the variables j 1 , ..., j s constitute a system of convex coordinates for C Ω and let be the appropriate submatrices of V Ω . The set SB (Ω) is defined by the algebraic conditions: Notice that we get a strong reduction on the degree of the polynomial expressions defining the set of symmetry-breaking states. The latter is the case given that: • We eliminated all the nonlinearity that was implicit in the definition of the set of steady states. • We reduced the degree of the instability condition. The instability condition is given by the Hurwitz-Routh inequalities [18], and those inequalities involve the determinant of the matrix to be analyzed. The entries of J Ω are nonlinear polynomials, and the determinant of this symbolic matrix is a polynomial whose degree can be larger than 2k. On the other hand, we have that the entries of AV Ω − BV Ω are linear polynomials, and we have that the degree of det (AV Ω − BV Ω ) is upperbounded by k.
However, if k is large, and despite all the reductions achieved so far, the stability analysis of matrix AV Ω − BV Ω can be unfeasible. A full stability analysis of matrix AV Ω − BV Ω presupposes the computation of its determinant, as well as the computation of other large subdeterminants. Take into account that computing symbolic determinants requires exponential time. If we want to avoid the computation of those large (symbolic) determinants, we will have to conform ourselves with a heuristic algorithm.

A Heuristic Algorithm for CPVEH Based on SNA and the MM-Condition
Let Ω be a pseudochiral network of order k, let AV Ω − BV Ω be the symbolic matrix that we want to analyze and let be the characteristic polynomial of this matrix. The coefficient Ω i is equal to the sum of all the diagonal subdeterminants of AV Ω − BV Ω of order i, see [22]. Notice that all those coefficients are polynomials over the set of parameters that we are analyzing, and notice that we are interested in determining the values of those parameters that make the latter (parameterized) polynomial become unstable. We must take into account the following fact: the roots of a polynomial are determined by its coefficients. Let us consider two instances of the latter phenomenon: 1. If the inequality (−1) k Ω k < 0 holds, the polynomial p Ω (λ) has a positive real root. 2. p Ω (λ) can have positive real roots only if it has negative coefficients.
We could focus on the first item and conform ourselves with a sufficient condition for instability, or we could focus on the second item and conform ourselves with a necessary condition for the existence of positive real roots. Notice that Ω k is equal to det (AV Ω − BV Ω ) and recall that we wanted to avoid the computation of large subdeterminants of the matrix AV Ω − BV Ω . This latter observation leads us to focus on the second item and consider the set SB * (Ω), which is the subset of SB (Ω) constituted by the steady states that satisfy the condition: there exists i ≤ k such that (−1) i Ω i < 0. We must ask: what do we gain if we focus on the second criterion? It was observed that in most cases, the chemical instabilities of Ω are determined by small subnetworks, see [23]. Moreover, given i ≤ k, the influence of all the subnetworks of size i is encoded in the coefficient Ω i . We can conclude that in most cases, the chemical instabilities of Ω are encoded in the first coefficients of p Ω (λ). We can use the latter as a further heuristic principle which tells us that: for most pseudochiral networks exhibiting homochiral dynamics, there exist small values of i such that the inequality (−1) i Ω i < 0 gets satisfied for nonnegative values of the parameters. Therefore, we focus on the problem: where Ω is a pseudochiral network. • Problem: compute the minimum i for which the inequality (−1) i Ω i < 0 gets satisfied for nonnegative values of the parameters, and sample this set of nonnegative solutions.
Our heuristic approach for solving CPVEH consists of solving the problem approx-CPVEH. Consider the following algorithm that we denote with the symbol SNA-sampling.
Algorithm SNA-sampling works, on input Ω, as follows: 1. Compute the matrices AV Ω and BV Ω . 2. Given i ≤ 5 determine the set of nonnegative solutions of the inequality (−1) i Ω i < 0. If all those sets are empty reject Ω; otherwise, sample the union of those five sets. 3. Given j = (j 1 , ..., j s ) an element of the computed sample, compute a racemic steady state of Ω, say s j , such that j is the tuple of convex coordinates of s j . Do the same with all the members of the computed sample.
We notice that: • If Ω is a pseudochiral network of order 1, the algorithm correctly and efficiently samples the set SB (Ω).

•
If Ω is a pseudochiral network of order k ≤ 5, then the algorithm samples the set SB * (Ω) . • If Ω is a pseudochiral network of order k > 5, then the algorithm approximates the set SB * (Ω) by a semialgebraic set that is defined by polynomial expressions whose degree is upperbounded by 5.
SNA-sampling is a heuristic algorithm that is as efficient as possible, solves the problem CPVEH for pseudochiral networks of small order, and allows us to compute important information in the case of higher orders. We exemplify, in the next section, the power of this tool with the analysis of three pseudochiral networks that were taken from the literature dedicated to the study of biological homochirality.

Computer Experiments
The algorithm SNA-sampling, as described above, was implemented as a computer program called Listanalchem [11], option six. The algorithm determines if a chemical network Ω, given as input, can produce spontaneous mirror-symmetric breaking (SMSB), and if so, samples the set of parameter values that can produce those dynamics. The computed samples are used to make numerical simulations of the dynamics. For this task, we use Chemkinlator [24], a software tool that can also build bifurcation diagrams. Actually, we used bifurcation diagrams to find the SMSB region for a model which presented problems for the heuristic used in sampling. The construction of the bifurcation diagrams could be a necessary additional step, given that, depending on the model, the instability region can be so small that numerical error could place us out of, but close, to this region. Also, the heuristic used in the algorithm can cause the results to be out of the instability region. In those cases, a fast exploration around the computed values, using bifurcation diagrams, could be enough to find SMSB. Additionally, this procedure can be used to build phase diagrams such as the ones shown in references [25,26].
It is worth remembering that as mentioned in Remark 10, the implementation of the algorithm uses the extended stoichiometric matrix with rows that encode the duality of the pairs of reactions that involve enantiomers. We observed that after adding those constraints the number of extreme currents got reduced. The details of this process, including the extended stoichiometric matrix, and its manipulation, can be seen in the output of Listanalchem option six. These outputs are available in the supplementary material.
Finally, we would like to remark that the units used for simulations are arbitrary. We do not set units to avoid huge or tiny numbers. Instead of that, we use numbers in the interval [0,2]. This particular way to perform the analysis, for each model, does not change the qualitative behavior of the models, and it helps to show clearer images of the relevant facts. Particular units can be obtained using the corresponding factors in concentrations and rate constants, but that fact will be not explored here because we are interested, first of all, in the stability of the models that generate SMSB.

The Replicator Model of Hochberg and Ribo
Hochberg and Ribo [27] have investigated the network described below.
This model has two enantiomeric pairs, ( 1 R D , 1 R L ) and ( 2 R D , 2 R L ). This means, according to the previous definitions, that it is a pseudochiral network of order 2. The analysis of the model, using the developed algorithm, confirms its ability to generate SMSB. However, it is important to remark that the outputs obtained with Listanalchem did not always produce SMSB immediately. In those cases, bifurcation diagrams, around the computed values, allowed us to compute the range of values that produce the symmetric breaking. Those bifurcation diagrams were done using Chemkinlator [24]. Figure 1 presents a typical example of the described situation, including three bifurcation diagrams. The right values are easy to find from the values given by the algorithm; for example, a fast way to find SMSB is by exploring the velocity of the entrance of A (∅ k 12 − − → A), see Figure 1B. In this way, it is easy to tune the region of SMSB.

The APED Model
Plasson et al. [25] have proposed the noncatalytic model presented below.
Our algorithm and software can find values of the rate constants for which the APED model exhibits the breaking of mirror symmetry. We could compute those values without the need for any additional work, as was necessary with the previous model.
The full set of rate constant, for this particular simulation, can be seen in the output of Listanalchem presented in the supplementary material. However, it is important to emphasize that k 0 = k 2 = k 4 = k 6 , k 1 = k 3 = k 5 = k 7 , and k 8 = k 9 = k 10 = k 11 = k 13 . Also, an initial enantiomeric excess, Additionally, we would like to remark that the approach developed here to obtain the instability regions has a solid mathematical background which makes it more efficient and general than the brute-force approach used in the original paper of Plasson et al. [25]: a systematic scan of the rate constants. Figure 2 shows an example of the results given by the algorithm.  It is interesting to remark that Plasson et al. [25] claimed that each one of the three sets of reactions: must have different rate constants, and because of this, they introduce the parameters α, β and γ, which are bigger than zero and different from 1. We found that those constraints are necessary because given α = β = γ = 1, the algorithm finds that the SMSB is not possible. Observe that under the equivalent condition, see the legend of Figure 1, the previous model (Replicator Hochberg-Ribo) exhibited SMSB.
Iwamoto considered to be variables only the species A, L, D, E L and E D . The Iwamoto model under perfect conditions shows SMSB. The instability region is a tiny one as shown by Figure 3. However, the algorithm could sample this region without problem.
On the other hand, the Iwamoto model under imperfect conditions is stable if the stereoselective (R1, R1a, R2, and R2a) rate constants satisfy the equalities k 2 = k 2a = k 4 = k 4a and k 3 = k 3a = k 5 = k 5a , and at the same time the stereospecific (R3, R3a, R4, and R4a) rate constants satisfy k 6 = k 6a = k 8 = k 8a and k 7 = k 7a = k 9 = k 9a . However, if those rate constants are different, as the author assumes, the model presents SMSB, which means: k 2 = k 4 = k 2a = k 4a , k 3 = k 5 = k 3a = k 5a , k 6 = k 8 = k 6a = k 8a and k 7 = k 9 = k 7a = k 9a . This latter condition (different rate constants for particular sets of reactions), seems to be necessary to obtain the desired results in those particular models, but it is not clear that it is a condition that can be presupposed of the prebiotic environment. The results for the Iwamoto model under imperfect conditions can be seen in the supplementary material and in the output of Listanalchem. [L] [D]

Discussion
The developed algorithm and the implemented software are powerful tools, capable of predicting the stability (instability) of models proposed to explain the emergence of homochirality in biological systems. The solid mathematics, behind the algorithm, makes it a robust tool to find the instability regions of chemical reaction mechanisms (networks models) of biological homochirality. This tool can be used to establish the particular rate constant values (intervals), for which a model breaks the initial racemic mixtures. The algorithm and the software can be used to generate phase diagrams of the models proposed to explain the origin of homochirality. Using the previous information, one could study the structure of those models proposed to explain the origin of homochirality in prebiotic earth.
It is important to remark, at this point, that we have focused on the homochiral dynamics that can be triggered by tiny perturbations of the racemic states. Enantiopure states can also be reached by alternate routes, like, for instance, the dynamics that are triggered by large perturbations of those states. We did not study the effect of large perturbations, given that we do not count with the required mathematical tools.