(1) We consider a set of monomeric molecules (monomers) of different nature labelled from 1 to n.

(2) These monomers are ‘nutrients’ that are present outside a ‘reaction chamber’ or ‘cell’.

(3) Different numbers PD (for Polymerisation Degree) of these monomers can be assembled in linear (non-branched) polymeric molecules (polymers) to give different lengths. The symbol P_j({k}) represents a polymer containing PD=j monomeric units and {k} is an ordered set of j symbols, each symbol being an element of the set of n monomeric units defined above. For example, if the polymer P_j({k}) is the string 23112, then j=5 and {k}={2,3,1,1,2}. Polymers containing PD=j monomeric units therefore exist in a maximum of n^j different permutations.

(4) A polymer containing PD=p monomeric units may be cleaved or added to other polymers (PD=q) or monomers (PD=1) such that the order and total number of monomeric units is conserved in the reaction. Reactions between molecules are of the form, P_p({k})⊕P_q({l})↔P_p+q({k}{l}) where {k} and {l} are ordered sets of the monomeric units and ({k}{l}) is the result of the addition of these sets in the order k to l. It is therefore the equivalent of reactions of addition (left to right) or cleavage (right to left). This reaction is reversible but is not commutative, i.e., P_p({k})⊕P_q({l}) is not the same as P_q({l})⊕P_p({k}) and P_q({l})⊕P_p({k})↔P_p+q({k}{l}) is not a valid reaction.

(5) Reactions are catalysed by a particular polymer P_p({k}) that may itself be a reactant (autocatalysis) and that we term ‘enzyme’ for convenience. More than one variety of enzymes, e.g., P_p({k}) and P_q({l}), may separately catalyse the same reaction. A single variety of enzyme may catalyse more than one reaction. No monomer may catalyse a reaction. Some polymers do not catalyse reactions.

(6) Reactions of the above type are only studied in the confined volume of a cell or reaction chamber that, in its initial form, we regard as created by the self-association of a random number of each monomer and a random number of a random selection of polymers made outside the cell (by ‘abiotic’ mechanisms that may be different from those in (5) and that we do not study).

(7) Nutrients are then supplied to the cell at regular or intermittent intervals. The cell may also be supplied with polymers, but at a much lower rate than that of the supply of nutrients.

(8) The dynamics of the system is described by representations of its state at discrete time steps. At each time step, a nutrient may or may not be incorporated into the cell, depending on the availability of the nutrients outside the cell. In principle, there could be an efflux of material through the cell, since the possibility also exists of losing monomers and polymers from the cell (for example, the probability of a monomer or polymer being lost could be inversely proportional to the number of reactions in which it is involved). This possibility is not implemented in the version described here.

(9) At each time step, two lists that describe the system are updated. The first list, the MoleculeList, contains a reference to each molecule present in the cell. Each molecule has a description comprising its label (name), the number of copies, whether it is a monomer or a polymer and in the latter case, whether it is an enzyme. The second list, the EnzymeList, contains a reference to each enzyme present in the cell. Each enzyme has a description comprising its label, the number of copies, its activity status (the possibility of using the parameter ‘active or inactive’ is built into the model but has not been used to generate the data shown here) and the reaction(s) it catalyses. A reaction is defined by three molecules, corresponding to two substrates, P_p({k}) and P_q({l}), and one product, P_p+q({k}{l}), and by the k_f and k_r for this reaction where k_f and k_r are the equivalent of the rate constants for the forward and reverse reactions, respectively. An enzyme can catalyse more than one reaction. Initially, the number of varieties of enzyme is chosen at random. Then the cell is fed in accord with (7).

(10) At each time step, the system is updated by calculating whether each variety of enzyme catalyses its cognate reactions. Each variety of enzymes is examined. The forward reaction can take place if the enzyme P_r({m}) is present and active and if k_f×N({k})×N({l})>k_r×N({k}{l})×N_t, where N_t is the total number of molecules in the cell and k_f and k_r are the rate constants for the particular reaction as catalysed by the enzyme. If the inequality is reversed, the reverse reaction occurs. Nothing happens at equilibrium.

(11) Iteration during the time step. A variety of enzyme is chosen from the EnzymeList either according to its concentration (as presented here) or at random (in this case, the choice is weighted by the number of copies, N_m, of the enzyme of a given variety, P_r({m}), such that the probability of choosing this variety is proportional to N_m/N_Et, where N_Et is the total number of enzymes).

The total increase and decrease in the number of copies of each molecule involved in a single catalysed reaction during the time step is obtained from:

(12) At the start of each time step, the cell is tested to see whether it has grown to a threshold at which cell division could occur. This test is based on the total number of polymers in the cell but alternatives include the number of monomers in the form of polymers as well as the number of copies of a specific polymer. In the present model, all reactions cease at this critical cell size and the program ends.

Here, we present a typical run of the program using only three types of monomer. The number of copies of these monomers present in the cell is given in C₀, for concentration at origin (column 4, Table 1). The open circle symbol in column 1 (Table 1) indicates which of the polymers in ID, for identity (column 3), were present in the cell at the start of the experiment. C₀ also gives the number of copies of polymers at the start – note that this cell did not initially contain polymers of types ID 233 to ID 31321131212313211121. C_f, for final concentration (column 5), gives the number of copies of the polymer at the end of the experiment – note that if a polymer has neither an entry in C₀ nor C_f, it is because it was formed during the experiment but disappeared before the end. N_s for node class of substrate (column 8) indicates the number of different reactions in which a molecule was involved as a substrate – in other words, it is the node class of the molecule with respect to the reactions in which it has been consumed; N_p for node class of products (column 8) indicates the number of different reactions in which a molecule was involved as a product; N_t for total number (column 7) is the sum of N_s and N_p. We can therefore consider each molecule as a node belonging to a network of either substrates or products or both; a molecule that can be made into 3 different molecules (i.e. is a substrate) and that is made from two others (i.e. is a product) has a node class of N_s=3, N_p=2 and N_t=5. In this analysis, the cell can be regarded as a network of molecules that are connected by catalysed reactions. Reactions (column 9) is the number of times the reaction occurred in the course of the experiment. The right arrow (column 10) is the number of times the molecule was consumed, the left arrow (column 11) is the number of times the molecule was produced and C_n is the difference between consumption and production.

The asterisks in ez (column 2, Table 1) indicate which of the polymers were enzymes. Polymers either started as enzymes or acquired activity later. This latter is simply a device to model the entry of a new enzyme by ascribing at random an activity to a polymer present in the original cell (column 1, Table 1); an alternative possibility not explored here is that acquisition of activity is due to a cofactor entering the cell.

The kinetics of the numbers of monomers and polymers (Fig. 1) show that an event occurred at time step 360 when enzyme 11 ‘entered’ the cell. The reaction catalysed by 11 is 2+3↔23 and this led to a rapid increase in the numbers of 23 (Fig. 2). There were no enzymes catalysing the reaction 1+1↔11, so catalytic closure did not occur [1]. Nevertheless, the system grew using its initial complement of polymer 11 until time step 688 when it had doubled the number of its polymers (from 68 to 136). 23 was also an enzyme and catalysed the reaction leading to another enzyme 33+21↔3321; 3321 catalysed 2+33↔233; 233 catalysed 11+31↔1131, a reaction that consumes or generates 11 (Table 2). Hence there is a cycle. (Note that the reactions of all enzymes are easily displayed, although we do not show them here; many other reactions are indicated in Table 1.) It might be expected that enzymes that catalyse the addition of monomers to other molecules would be the most abundant. This is not the case in this example, where 11, present in six copies at the end, is not as abundant as 23, present in a hundred copies at the end. Note that another enzyme in the cycle, 3321, also catalysed addition of the monomer, 2 (Table 2), and that this enzyme was present neither at the start nor the end of the experiment. Such fluctuations in potentially key enzymes were common in the program.

The numbers of members of a particular node class (where the class is the sum of the reactions in which the molecule is either a substrate or product) reveals a distribution with a long tail (Fig. 3). The molecules with the highest connectivities were 21 and 1121 and those with the next highest were 2, 313 and 31321 (Table 1). We have not analysed connectivity in terms of catalysts in a context where an enzyme could catalyse several reactions, although this would be possible. In summary, these results show that even in this simple model, autocatalytic networks form readily. It should be noted that, reassuringly, the system does not grow if it does not use the monomeric food set in its reaction network!

Artificial chemistry offers a powerful possibility for both testing existing biological concepts and for deriving new ones [2,10]. To contribute to an ambitious, integrative modelling of biological cells, I-cell [11], we have used a variant that we term artificial microbiology to study how autocatalytic networks develop when they are confined to cells that are allowed to grow to a fixed size. The results from a limited range of simulations, with a restricted set of molecules, are compatible with the idea that the most abundant enzymes are those that catalyse the addition of monomers. This appears obvious – as a network evolves, the first limiting step in the growth of the cell is the flow of monomers into the cell and the second limiting step is the addition of these monomers onto polymers or onto other monomers. The enzymes that add monomers may therefore correspond to the ancestral precursors of modern polymerases and ribosomes, which also add monomers onto polymers; such enzymes can constitute a major proportion of cell mass. That said, the proportion of cell mass in the form of ribosomes varies with the nutrient supply in bacteria such as Escherichia coli and it should be noted that the most abundant enzyme in the example given in this paper does not catalyse monomer addition.

Abstract networks or graphs consist of nodes and the connections or links between them. A particular node belongs to a node class that corresponds to the number of connections per node. By allocating a metabolite to a node class on the basis of its connections to other metabolites, we and – independently – others found that real metabolic networks have a power law distribution of node classes characterized by a long tail [7,8,12]. Similar distributions have been found for proteins [13]. Artificial microbiology permits the origin and properties of real metabolic distributions to be investigated. In many individual runs of the system presented here, the distribution of the connectivities did indeed show a long tail. Hence, the system may be useful in answering such questions as: What network structure confers robustness to fluctuations in nutrition? Where should the highest node classes be located so as to assure a maximum use of inputs and a steady output? What allows the network to be regenerated after starvation?

In the version presented here, the program ends when the cell attains a size at which division might occur. The importance of cell division in the evolution of autocatalytic networks has recently been described [3,5]. Our program offers the possibility of exploring how competing autocatalytic networks within the mother cell may be separated by cell division into individual ones. By selecting the faster growing daughter cell, the role of cell division in the evolution of networks can then be studied. This role may be facilitated by the co-localisation of many different macromolecules into a non-equilibrium hyperstructure in response to the cellular need to perform a function [14,15]. Cell division in bacteria, for example, would be performed by a hyperstructure comprising division genes, their mRNA and enzymes, together with particular lipids and ions such as calcium [5]. Future development of the model will therefore entail attributing increased probabilities of reactions to polymers that are colocalised so as to allow evaluation of the consequences of hyperstructure formation. Finally, artificial microbiology may be extended to study what happens when other types of reactions are introduced [16], when the cell shifts from a growing to a non-growing state, when inactive enzymes are preferentially inhibited, and when a coding polymer (DNA/RNA) is present.