Introduction

A central effort in theoretical biology and ecology is to provide an effective description of the intimate but often subtle relationship between a given environment and the evolution of its ecosystem. In the case of simple environments without geographical isolation and physical barriers, as the one here considered, the emergence of species and phenotyping clustering Davis (1943); Rundle and Nosil (2005); Keymer et al. (2012); Gupta et al. (2021) is generally considered an outcome of the competition for the limited available resources Dieckmann and Doebeli (1999); Pigolotti et al. (2007); de Aguiar et al. (2009); Anceschi et al. (2018). The coexistence of stable species is ecologically understood in terms of “niches”, indicating the unique role and position that a particular species occupies within an ecosystem. According to the “niche hypothesis” Chase and Leibold (2009); Peterson (2011); Anceschi et al. (2018), biodiversity is limited by the number of “niches” (or types of resources) that are present since no two species can occupy the same niche indefinitely, as one would eventually out-compete the other Levin (1970); Gupta et al. (2021).

Fitness, which quantifies species reproductive success and thus also their relationship with the environment De Visser and Krug (2014), cannot generally be defined in a predictive way even in the most idealized systems Thurner et al. (2010); Wiser and Lenski (2015). Indeed, despite the efforts to identify simple case study systems, the evolution of populations formed by a variety of species remains difficult to model and to quantitatively characterize because of the inherent complexity of ecosystems and living beings, of the large number of potentially relevant variables of difficult access, and of the role of stochasticity Catalán et al. (2017).

In this scenario, introducing new tools to explore and test eco-evolutionary models, concepts and interpretations appear as the best strategy toward new understanding Solé (2016). Along this line, various synthetic biological platforms have been proposed in the last years Ichihashi et al. (2013); Tizei et al. (2016); Parrilla-Gutierrez et al. (2017); Kauffman et al. (2018); Adamala and Szostak (2013); Katla et al. (2023) that exploit different principles and mechanisms, and focused on in-vivo, in-vitro, ex-vivo, or in-silico approaches. For example, in Ichihashi et al. (2013) a “cell-like” model system is proposed, in which the evolution of one long genomic RNA (> 2,000 nt-long) was investigated in detail under the action of a selective pressure ultimately provided by its biological meaning. Despite the tremendous simplification provided by this approach, the constructed artificial cell is still a complex system that includes ribosomes, lipids, translation factors, accessory proteins, tRNAs and amino acids. While most of the synthetic biological platforms are based on the establishment of cell-like, compartmentalized systems involving complex biomolecular milieus and processes Tizei et al. (2016), a few are devoted to investigate, with distinct strategies, the competition between two coexisting vescicle-based species competing for resources (vescicle-forming molecules) in a non-biochemical context Adamala and Szostak (2013); Katla et al. (2023). These articles demonstrate that typical evolution concepts, such as competition or resources and niche exclusion principle, can be applied to “non-biological” systems.

We here propose a synthetic eco-evolutionary scheme, with no cell-like compartmentalization and no connection to the molecular biology of the cell: no coding sequences, no translation, no proteins involved. We consider an ecosystem formed by a large crowd (∼ 10¹⁵) of distinct molecular individuals interacting with each other and competing for survival in an environment with fixed resources, and we focus on the process by which selection and competition drive the emergence of dominating species. Specifically, we study solutions of 50-base-long single strand DNA oligomers with random sequences, thus initially all distinct, whose mechanisms of survival and mutual interactions are based on DNA hybridization. By exploiting the wealth of tools and knowledge around DNA code selective pairing and DNA synthesis, amplification and sequencing, we create a condition in which a limited number of accessible and computable variables controls the destiny of an ecosystem formed by biological molecules that evolve in a non-biological way, i.e. with no reference to the biological meaning of their sequence.

The results of our work, demonstrating the intimate connection between fitness and ecological interactions, belong to the growing body of studies Fussmann et al. (2007); Vetsigian (2017); Camacho Mateu et al. (2021) showing that ecological and evolutionary processes are strictly connected in the emergence and maintenance of species.

Results

Affinity-based DNA Synthetic Evolution

We introduce here a variant of SELEX for in-vitro synthetic evolution of oligonucleotides to develop protein-binding aptamers Ellington and Szostak (1990); Tuerk and Gold (1990). In standard SELEX protocol, the evolving oligonucleotides are selected at each cycle by their interactions of with the target protein. In the experiment here reported, we implemented a selective mechanism based on the affinity capture provided by magnetic beads carrying single-stranded DNA (ssDNA) filaments of fixed length L = 20 and sequence, that act as targets (or resources, Figure1a). Selection is thus primarily based on the sequence of the DNA individuals and its level of complementarity to the targets. This marks a significant difference with SELEX, in which the aptamer-protein interaction depends instead on higher order factors such as the secondary structure of the oligonucleotides and the variety of binding sites on the folded protein. Being these hard to model, predict and control, SELEX has never been considered, to the best of our knowledge, as a useful experimental test bench to understand evolution.

In our Affinity-based DNA Synthetic Evolution (ADSE) protocol, evolution start from an initial pool of DNA Individuals (DNAi), chosen to be of fixed length L = 50 and random-sequence. Since the potential molecular variety is 4⁵⁰ ∼ 10³⁰ while our experiments use about 10¹⁵ initial molecules, each sequence is represented one time only. The following evolution process is sketched in Figure1b, Figure 1-figure supplement 2, and is given by three steps.

Figure 1.

1. Selection: the seed population is mixed with a given amount of dispersed capture beads. After a suitable incubation time, the beads are extracted from the solution and the bound oligomers released and saved. The rest of the original solution is discarded;

2. Amplification: the pool of “survived” oligomers is PCR amplified about 1000 times to recover the initial molarity;

3. Sequencing: a small portion of the amplified sample (1-3 10⁶ molecules) is analyzed by massive parallel sequencing.

These steps constitute a cycle - one generation of evolution - that we repeated up to 24 times in two independent evolution histories.

In actuality, to enable amplification and sequencing, DNAi are built by flanking the 50mer with two 25 base long fixed sequences that enable primer binding for a total length of 100 bases. Such two terminal segments of the DNAi are made inactive during selection by hybridization with oligomers of perfect complementarity, as sketched in Figure 1a, Figure 1-figure supplement 1.

A key feature of ADSE is that DNAi can interact not only with the resources, but also with itself and with each other, and also potentially form complexes binding to the affinity beads as summarized in Figure 1c. Indeed, the choice of a length of 50 for the DNAi interaction was thought to enable - in principle - simultaneous resource and mutual binding.

PCR amplification can be of high fidelity or error-prone, the latter choice enabling genetic drift and, potentially, speciation. Since our primary goal was to first investigate fitness in a pool of competing species within a given niche - defined in our case by the 20 base long resources, we opted for a high fidelity amplification, and left the investigation of high mutations regimes to a follow-up experiment. Because of the many PCR cycles required in the ADSE scheme, the amplification inevitably leads to the formation of artifacts in the form of longer sequences Tolle et al. (2014), a phenomenon that intrinsically sets a limit to the number of generations that can be explored (see Appendix 1).

As detailed below, the ADSE protocol enables observing a non-trivial evolution of the DNAi ecosystem across generations, the emergence of dominating species and nonmonotonic population evolution. It also enables to appreciate the role of inter-specific interactions and their contribution to fitness.

Evolution of the DNAi ecosystem

The main output of the synthetic evolution that we are proposing is the dataset {DNAi}_j obtained by sequencing DNAi at the various cycles (1 ≤ j ≤ 24 being the index of the cycle). We find that the initial random ecosystem markedly changes with ADSE generations as shown in Figure 2a, in which we have plotted the evolution of: (i) the fraction F_D of the total population formed by distinct nucleotide sequences (red dots). F_D drops from 1 to nearly 0, indicating that initially DNAi are all different from each other, while after 24 generations the number of distinct sequences is much smaller than the number of DNAi. This indicates that most of the initial sequences become progressively “extinct”; (ii) The fraction F₁₀ of the total population formed by the 10 most abundant DNA sequences (black dots). F₁₀ is initially close to 0 - being each sequence represented by a single DNAi, and grows to about 25% at cycle 12 and 55% at cycle 24 - indicating that the offspring of 10 initial “mother” sequences becomes the majority of the system. This loss of diversity can also be quantified by the zipped file size of the list of sequences Benedetto et al. (2002) or by computing its Shannon entropy Hill (1973), both markedly decreasing during evolution (see Appendix 1).

Figure 2.

Since in ADSE survival depends on bead capture, we expect the DNAi-resource binding strength to be an essential ingredient of fitness. To this aim we computed the mean DNAi-resources binding free energy ΔG_DR across generations, plotted in Figure2b (blue squares, right-hand side y-axis). Specifically, ΔG_DR has been computed based on DNAi and resource sequences by using the standard “nearest-neighbor approximation” for the thermodynamics of DNA hybridization SantaLucia Jr (1998); Ghosh et al. (2019); Plata et al. (2021), as implemented in the NUPACK tool in Python Zadeh et al. (2011). Being this computation of some complexity, we could perform it only on batches of 1000 DNAi randomly chosen within {DNA_i}_j, the error bars expressing the uncertainty introduced by such downsampling (see Material and Methods section).

To perform a faster - and thus complete - analysis of the resource-binding energy in the evolving population, we introduced a simpler, albeit more approximate, quantifier: the length ω of the longest consecutive number of basis within each DNAi that is complementary to the resource sequence (considering all possible relative positions of the two), in which we allow for single pairing errors and 1-base bulges Mambretti et al. (2022a). While ⟨ω⟩, missing key ingredients such as the fraction of CG vs. AT base pairs and their sequence, is certainly inadequate to evaluate the binding strength of a specific DNAi, we find that its average value ⟨ω⟩computed for the whole population in each generation (Figure2b, full red dots, left-hand side y-axis) grows almost identically to ⟨ΔG_DR⟩. The ⟨ω⟩ axis has been scaled so to make its value computed on a pool of ideal randomly sequence DNAi (dashed red line) match, on the plot, ⟨ΔG_DR⟩ computed on the same population (dashed blue line). A more detailed comparison between ω and ΔG_DR is in Figure2c, showing box plots with ΔG_DR on the y axis and ω expressed through color code, built on a random selection of 1000 distinct DNAi in the initial (left hand side box) and the final populations (central box), and on top 10 most frequent sequences in the final generation (right hand side box). The result further strengthens the validity, in our statistical context, of ω as a quantifier of interaction strength to the resources Mambretti et al. (2022a).

Figure2b reveals that, during ADSE, ω tends to saturate at a value ω_sat ≈ 10, indicating that no relevant selective advantage is gained when ω > 10, in agreement with the notion that the residence time of hybridized oligomers becomes larger than typical experimental times when ω > 12 Di Leo et al. (2022).

Figure 2d-f describes the evolution of the ecosystem {DNAi_j by showing P (ω), the fraction of DNAi having overlap ω with the resource evaluated in the initial pool (panel d) and at generation 12 (e) and 24 (f). P (ω) clearly evolves, its small ω components being progressively lost, while individuals with large ω grow in number, as they are more successful in being selected and amplified. It might be worth pointing out that the appearance of non-zero p(ω > 13) at generation 12, while p(ω > 13) = 0 in the initial distribution, is a pure effect of the undersampling involved in the sequencing procedure.

The different destiny of DNAi with distinct ω values is shown in Figure 2g-i, where we show the temporal (i.e. across generation) evolution of the fraction of the total population whose ω is in the following ranges: 3 ≤ ω ≤ 5 (panel g), 6 ≤ ω ≤ 10 (panel h) and 11 ≤ ω (panel i). As expected, DNAi with weak affinity to the resources decrease while those with large affinity increase. Interestingly, DNAi with intermediate affinity exhibit a non-monotonic behavior, indicating that the conditions for survival change during the evolution, reflecting the evolution of the ecosystem.

Null Model and Eco-evolutionary Algorithm

In order to better understand the experimental outcomes, we first build a null random model without evolution and then an Individual-Based Eco-Evolutionary (IBEE) model enabling predictions for P (ω).

The null model describes the interaction between individuals in the initial pool {DNAi_j=0} of random sequences and the resources within a purely combinatorial framework (see Methods and Appendix 3 for mathematical details). In this model, we attach a random string of length 50 to the resource string in a random position, and we compute the maximum consecutive overlap, accepting the binding only if it is at least formed by 3 complementary basis. For analytical tractability we do not allow the random string to shift over the resource, but we calculate the ω constrained to a specific binding position.

The resulting analytical p₀(ω), plotted in Figure 2d (orange crosses) closely matches the data obtained from sequencing, confirming the random-sequence nature of {DNAi_j=0}. The analytical P₀(ω) extends to a ω range where no data are available because of the limited size of the sequenced pool. The comparison between P₀(ω) and P (ω) in panels e and f shows that the exponential decay of P₀(ω) at large ω is maintained even in generation 12 and partly in generation 24, further supporting the notion that the selective advantage of ω saturates at large ω.

In the IBEE algorithm we consider N_p individuals. Each has a fitness f_i(ω) (i = 1, 2, .., N_p) that depends on its affinity ω with the resource. At time t = 0 we assign ω to each individual based on P₀(ω) from the null model. Then, for each evolutive cycle we model competition and selection so that out of N_p individuals, only N_r survives. This is attained as a combination of two processes: a fraction x of the N_r sequences results from sampling individuals from the N_p population, each with a survival probability f (ω) ∈ (0, 1); the remaining fraction 1 − x is formed by a straight random selection from the N_p individuals. The latter group is meant to mimic “neutral drift” (as it would be expressed in evolutionary language) provided by non-specific binding to the beads. The N_r survivers are amplified by identical copying back to N_p, the starting population of the next evolutionary cycle (the introduction of very rare mutations do not affect the results). We perform 24 evolutionary cycles.

As expected, the outcome of the eco-evolutionary dynamics strongly depends on the shape of the fitness function. We model it as:

where ω_max is the (cycle-dependent) largest value of ω within the actual population and, ω^* = ω if ω < ω_sat, while ω^* = ω_sat otherwise. ω_sat and γ are parameters to be tuned in the comparison with the observed⟨ω⟩. ω_sat expresses the loss of fitness gain for ω > ω_sat yielding the saturation of⟨ω⟩. γ represents the strength with which f (ω) depends on ω, and thus the rate at which low ω individuals are discarded. γ may of course depend on time on account of the evolving ecosystem.

After grid search we find x = 0.9, indicating that the random drift contributes to about 10% of survival at each cycle, and ω_th = 10, in agreement with the observations prompted by the saturating behavior of ⟨ω⟩. We also find that the data cannot be approximated with a single value of γ, but they can be very well matched assuming γ = 3 for the first 5 cycles and γ = 1 for the remaining cycles, as shown in Figure 2b, green line. With the same choice of parameters, the IBEE fitness based model captures the decrease in the diversity of the DNAi ecosystem (Figure 2a, green and grey lines) and the temporal (i.e. across generations) evolution of the relative population abundance of DNAi in the three ω intervals in Figure 2, panels g, h and i.

The effect of the system size and initial conditions on the IBEE results, discussed in the Appendix 2, do not qualitatively change the models results. Error-bars on simulations have been obtained by averaging 20 independent runs, starting from the same initial conditions. As it can be observed, the variability among simulations is negligible.

The change in γ, and thus in fitness, necessary to reproduce the observed ⟨ω⟩ with the IBEE model is remarkable because it contrasts with the generally simple (fixed resource and low mutation rate) conditions of ADSE environment and indicates that survival is controlled by more than affinity to the resources, as discussed in the analysis below.

Self and mutual DNAi interactions are evolutionary drivers

While ω is certainly a key driver of the observed evolution, it is clearly not the only one. The fact that the top 10 most represented sequences come to be, in the last cycle, about 55% of the total population (Figure 2a) implies that, even among sequences with large ω, only a tiny minority come to dominate, while the largest part of them eventually disappear. Moreover, the 10 most represented sequences do not stand out for their particularly large ΔG_DR or ω, as noticeable in Figure. 2c. The same data are plotted in Figure 3a as a P (ΔG_DR) distribution to enable comparing the free energy distribution in the initial population (purple shading), in the final population (blue columns) and in the top 10 (black columns). These elements support the notion that survival and dominance must be also due to factors other than ω, which we thus explored.

Figure 3.

Figure 3b shows the value, across the 24 generations, of two other contributions to the total free energy, normalized to their value computed in random sequences, so to enable comparing their relative variations: ⟨⟨ ΔGself⟩⟩ (green symbols), the average unimolecular free energy, expressing the average strength of the internal folding of each DNAi; ⟨ΔG_DD⟩ (red dots), the total bimolecular DNAi-DNAi free energy, comprising both self and mutual interactions. For comparison, we also plot the normalized value of ⟨ΔG_DR⟩ (blue dots). As for the case of ΔG_DR⟩, ⟨⟨ ΔGself⟩⟩ and ⟨ΔG_DD⟩ are computed by randomly selecting 1000 individuals or pairs, respectively, from the population at cycle j and using the NUPACK tool to compute the values.

Figure 3b reveals that ⟨⟨ ΔGself⟩⟩ grows in time even more than ⟨ ΔG_DR⟩, but with a different progression: ΔG_DR grows faster in the first cycles to later saturate, while the growth of ⟨ΔGself⟩is more uniform.

Since ΔG_DR is computed as the free energy of the whole DNAi-resource structure, it includes contributions of self energy associated to hairpins in the DNAi. Thus the growth of ΔG_DR could actually depend on the growth in ⟨ΔGself⟩(but not the contrary). To investigate this possibility, Figure 3c displays the scatter plot between ΔG_DR and ⟨ ⟨ ΔGself⟩⟩, for a randomly chosen subset of DNAi at j = 24. In evidence are the 10 points corresponding to the 10 most populous sequences. The plot shows weak or no correlation, and a relevant shift with respect to ΔGself= ΔG_DR (black line), indicating that there is no dependence of ΔG_DR on ⟨ΔGself⟩ and indicating that these two quantities reflect two independent driving forces in the selection mechanism of ADSE. Their behavior suggests that in the first stages of ADSE the selection is mostly dominated by the affinity with the resources, while in the later stages the requirement of a stronger unimolecular folding becomes more important.

A similar scatter plot analysis for ⟨ΔG_DD⟩ yields a different outcome. Figure 3d compares ΔG_DD computed for two selected DNAi (individual “k” and “l”) and the sum of ΔGself of the same two DNAi. The apparent correlation indicates that, on average, a large part of ⟨ΔG_DD⟩ simply embodies the growth of self-energy, although the difference ΔΔG ≡ ΔG_DD,kl − ⟨ ΔGself⟩_,k − ⟨ ΔGself⟩_,l (orange arrow) is non-negligible. Figure 3b shows the behavior of ΔΔG across cycles (orange squares). Despite the resulting mild growth of ΔΔG might appear not relevant, it actually indicated that ADSE, independently of ω, select strings that have higher reciprocal affinity than a random DNAi set (see Appendix 2). Indeed, as we could have obtained a significant decrease of the same quantity. It should also be noted that mutual interactions might involve the unfolding of hairpins self-structures, in which case their strength would be much larger than ΔΔG.

Self-interactions and mutual interactions can compete with the binding of DNAi to the resources either when they involve the same nucleobases or through steric hindrance. Therefore, we could expect that the increase of ΔG_DR would lead to a decrease of both ⟨ ΔGself⟩ and ΔG_DD. Their unfore-seen growth in the ADSE process is hence an indication of the selective advantage they convey. We interpret this behavior as an indirect sign that mutual DNAi interactions, more than just a impediment, are a major deadly threats for their survival. To avoid them, DNAi need screening. Indeed, the strong growth of self interactions, and the mild increment in mutual interactions, could represent the emergence of “defensive” strategies, as discussed below.

In our search for evolution quantifiers other than resource binding strength we also found that selection favors binding to resource sequences close to their resource 3’ terminal, away from the bead surface (see Figure 4 - figure supplement 2). While this is expected, being that terminal much less constrained and in a less crowded environment, it also provide useful clues to further analysis.

Figure 4.

The evolution of DNAi species

While the trends in average free energies offer some insight into the nature of the ADSE ecosystem, it has also become clear that they do not enable full understanding, especially since the dominant species are not - as one could naively expect - at the extremes of the energies distributions. To gain further insight, we examined the change in numerosity of the groups of DNAi having equal sequence, which might be regarded as the “species” of ADSE. Accordingly, the 10 most populous sequences are the 10 most successful species. The evolution of the populousness, expressed as the fraction F of the total population that belong to a given species, is shown in Figure 4a for a limited set of (mainly) successful species (for a larger selection see SI). Species that are later becoming dominant “appear” in our analysis only at generation 5 in a single or a few copies. This is reasonable since we sequence a sample of ∼ 10⁶ DNAi over the much larger population of ∼ 10¹² survived individuals. Species are named after their ω value, combined with an index (e.g. “12.0”, “12.1”, “12.2” …) expressing their ranking in populousness.

ADSE species populations do not have standard patterns of evolution. This is true even when limiting the observation to the successful species in Figure 4a whose population growths are varied and do not reflect in any simple way the energy quantifiers. Conversely, species with similar energy quantifiers may have opposite fate, either becoming dominant or rapidly extinguishing or else displaying a non-monotonic population evolution. It is worth pointing out that all species start equal, each represented in the original pool by a single molecule, which on this scale would appear as a F ≈ 10⁻¹² (black dot in Figure 4a).

We thus decided to seek further understanding on the survival of the fittest in ADSE by inspecting the “natural history” of a few species (see Tables 2, 3).

13.0 and 11.0 (orange and gray dots in Figure. 4a) are, respectively, the most and second most populous species in the last generation. It might be relevant to mention that the replica of ADSE described in the SI leads to different dominant sequences, as expected on the basis of the distinct initial pool and on the key role of randomness and sampling in the first cycles. At generation 24, 13.0 and 11.0 have grown to more than 85% and 53% of DNAi with ω = 13 and ω = 11, respectively, a further indication that ω is only one of the ingredients of the fitness. 13.0, and 11.0 likewise, have a weak self-interaction (Figure 3c), corresponding to hairpins that provide a mild defense to mutual interactions, as schematically sketched in Figure 4, panel b1. The hairpins involve most of the bases that are complementary to the resource (green shading), leaving however a few unpaired bases that might act as toehold to initiate resource binding. The resource binding takes place at its 3’ terminus (panel b2). A bit puzzled about the success of these two species we also explored their capacity of mutual interactions and found that both of them are capable of forming homodimers, as shown in panel b3 for 13.0 (binding energy marked in Figure 3d). Such dimers can bind resources at both ends (green shading regions), giving 13.0 and 11.0 self-screening and divalent binding capacity. We argue this to be a crucial drive toward success of these two species, which also explains why their growth has such an increment in the latest generations, following the probability of dimer formation.

The third most populous sequence (species 8.0) is characterized by a weak resource interaction which includes a bulge (see Figure 3c). While its survival rationale appears clear - good defensive self-binding (panel c1) and interaction targeting the 3’ end of the resource (panel c2) - it is hard to accept this could be the cause of the success. Again, inspecting the mutual interaction we found that 8.0 efficiently interacts with 13.0 (as marked Figure 3d), in a region (pink shading) that does not conflict with the capacity of 13.0 to bind to the resource (panel c3). We speculate that this parasitic capacity of 8.0 adds to the weak intrinsic binding to lead to a remarkable success of this species.

As a test to this concept, we focused on species 5.1, whose binding strength to the resource is the weakest in the top 10 species and one of the weakest among all the survivors of generation 24 (Figure 3c), a situation that cannot justify its success. However, by inspecting mutual interactions, we find its binding to 11.0 to be strong and stable (pink shading panel d1 and Figure 3d). An analogous behavior is found for species 5.0 and its interaction to 12.0 (panel d2), confirming that parasitism is an emerging successful survival solution in ADSE.

We then inspected 12.1 since it has largest ΔG_DR and ⟨ ΔGself⟩ among the top 10 species (Figure 3c). Species 12.1 has all the features to succeed in ADSE: it forms a weak hairpin with the bases involved in the binding with the resources and a strong hairpin that protects the rest of the ssDNA segment (panel e1); it strongly binds to the resource with the bond ending at its 3’ terminal (panel e2); it forms homodimers capable of double resource binding, analogous to that of 13.0 (panel e3). Remarkably, 12.1 emerges only late in ADSE (Figure 4a, blue dots), a very unusual behavior when compared to the other many species we have considered. The combination of this late appearance and of the remarkable growth in F afterward suggests that 12.1 is the outcome of one of the rare mutations that can occur even with the high quality PCR that we adopted. By assuming the initial growth as with the dominant species (black arrows in panel a), we argue that such mutation occurred around generation 9, where we have found ancestors having sequence equal to 12.1 except for 4 bases.

F has a monotonic increase for most of the species that become dominant in the latest generations. This behavior is however not at all general. Non-monotonic F, growing in the first generation and decreasing to extinction afterward, is the general behavior for the majority of species that do not become extinct in the first cycle, as also shown - as mean behavior - in Figure 2h. An examples of these is species 7.X, which we plot in panel a of Figure 4 (empty dots).

More intriguing is when a non-monotonic behavior is observed for one of the dominating species. This is the case of 11.1, whose growth can be again attributed to an effective screening and a good binding to the resources, but whose decrease is harder to justify. To understand it we again investigated mutual interactions, and we found a strong 11.0-11.1 bond (sketched in panel f), which, differently from the DNAi-DNAi complex examined above, competes with resource binding. We thus speculate that the fall of 11.1 is a consequence of the rise of 11.0, which has become frequent enough to make the formation of the complex 11.0-11.1 probable, which lowers the survival probability of both species. This could also explain why species 11.0 has such an irregular growth pattern, with a slowing down when 11.1 becomes highly populated.

The insight gained by these case-studies enables appreciating how self and mutual interactions can affect the fate of species in ADSE beyond what we could discern through distributions and average quantifiers. Indeed, this analysis shows that the success of the dominant ASDE species can be achieved through different combinations of resource and mutual DNAi interactions, a fact that makes the evolution of population so varied and generation dependent and explains why, in the latest generations of ADSE, the population distribution p(ΔG) (Figure 3a) is so irregular.

Discussion and Conclusion

We have here introduced ADSE, a synthetic molecular evolution protocol that exploits sequencing technology and DNA interaction computability to provide a test-bed for key concepts in ecology and evolution. We performed experiments by adopting the most simple environment within this protocol, i.e. by holding fixed the capture sequence and minimizing mutations, to achieve a condition dominated by competition and selection, which could enable investigating the nature of fitness in this simplified evolution process.

Our experiments and their comparison with theoretical models indicate that fitness is not a simple function of the direct competitive advantage of strong binding to the resources, as it could have been naively expected. Resource binding is indeed the dominant factor in the first part of the evolution (cycles 1-5), producing a fast growth of ⟨ω⟩ compatible with a strong dependence of the survival probability on ⟨ ω ⟩(f ∝ ω³). However, as ω reaches a value corresponding to bonds of moderate stability Woodside et al. (2006) - approximately from generation 6 onward - the selective pressure related to resource binding decreases (f ∝ ω), a condition enabling to appreciate other factors at play. Which are these factors is suggested by the different growth patterns of ΔG_DR and ⟨ ΔGself⟩, the former dominating the first cycles while the latter drives the later stages. However, the quantification of binding energies isn’t sufficient to predict the fate of individual species, that also depends on interaction details (location of binding and hairpin spots, secondary structures) and possibly on kinetics.

In fact, by design, DNAi can interact with themselves, both internally - forming hairpin-like structures - and mutually - leading to the formation of complexes. These interactions can either coexist-when they involve non-overlapping sequences - or they can become competitive. Our experiments reveal that, among the wide variety of conditions made available by the initial random seed, ADSE generally promotes those that combine a good interaction with the resources and a structure capable of screening from mutual DNAi interactions conflicting with binding to the beads. By exploring hybridisation-dependent DNAi secondary structures - the only ones accessible with simple analytical tools, we found three basic motifs: hairpins due to self-interactions (as for species 12.1), formation of DNAi homoduplets (as for species 13.0 and 11.0), and formation of dimers of DNAi belonging to distinct species. The latter can bring to distinct prototypical behavior: it can be the basis of parasitism, as in the case of the pairs 8.0-13.0 and 5.1-11.0, it can lead to extinction with no benefit for either species (as in the case of 15.0-12.0), or it can provide - at least in principle - mutualism. While we do not have a solid proof, clues suggest that cooperation might be the case of the interaction between species 5.0 and 12.0, since 12.0 has a weak self-defensive system and since the populations of 12.0 and 5.0 evolve very similarly (see Table 2). Finally, we cannot exclude a possible role of the formation of multi-DNAi multivalent complexes. This is the case of species 12.0 and 15.0 whose pattern of multivalent mutual interactions can, in principle, allow higher order interactions (see Appendix 2).

Fitness in ADSE is the outcome of a complex interplay of all these elements. The history of single species in Figure 4 indicates that the nature of fitness is beyond what can be captured by analyzing the probability distributions of the relevant parameters describing structure and interactions. Moreover, the variety of population evolution enlightens the fundamental fact that despite the constancy of resources, the progressive modification of the ecosystem brings about a change in the competition patterns and thus in fitness. A forthcoming work focused on a limited number of species will be devoted to better disentangle the nature of fitness in ADSE.

The evolution of the ADSE ecosystem, which corresponds to a marked decrement of its entropy - with 2/3 of the initial species become extinct in the first 5 generations - and potentially terminating with the indication of a single winner species leads instead, in the last cycles, to a significant number of dominant species all still growing at the expenses of subservient ones, indicating that in their direct competition none of them is strongly prevailing despite - or maybe thanks to - their distinct survival strategies and their very different population share. Hence, while the niche hypothesis could still be verified in the long run - beyond the experimental limitations due to PCR - it is clear that its drive is weak, suggesting that it could be overcome by environmental fluctuations, thus allowing for coexistence of different species even in a single niche environment.

Methods

Our experimental design takes advantage of a selective capture mechanism where magnetic beads carrying single-stranded DNA filaments of fixed length and sequence (resources) target DNA individuals (DNAi) (Figure 1 - figure supplement 1) present in a DNA library based on their level of complementarity. This process of selection is carried out through subsequent steps that are described in detail in the next paragraphs and represented in Figure 1-figure supplement 2.

Library design

The DNA library contains 100-nt-long sequences where a randomized central region of 50 nu-cleotides is flanked by 25-nt-long fixed sequences at both its 5’ and 3’ ends. The fixed regions provide an anchor point for primer annealing, required to perform PCR reactions during the amplification phase (see below). These terminal segments are made inactive by hybridization with oligomers of perfect complementarity (blockers), so that they are not involved in the selection phase. The blockers, as well as the fixed regions of the DNA library, have been designed to avoid any interactions with the resources. In addition, the blockers carry a phosphate group at their 3^′ end to prevent them from functioning as primers during the PCR amplification. Following the above described criteria, we designed two sets of sequences (Table 1), called Oligo1 (whose results are presented in the main text) and Oligo2 that was used as a replication experiment. All the oligonucleotides used in this work were purchased from Integrated DNA Technologies, Coralville, IA, USA (Table 1).

Beads preparation

The capture of DNAi within the DNA library was performed with carboxylic acid magnetic beads (M-270 Dynabeads, Invitrogen, Carlsbad, CA) coated with the resources. The resources are 20-nucleotide-long, 5’-amino-modified oligonucleotides. Their coupling to the beads surface was performed according to the manufacturer’s instructions. Following the activation and coupling procedure, the beads were washed in Tris-HCl (50mM, pH 7.4), and stored in the same buffer in single-use aliquotes.

Sample preparation

The starting samples were prepared in 1X SSC buffer (0.15M sodium chloride, 15mM sodium citrate, pH 7.0). In detail, the Oligo1 (or Oligo2) library (0.75nmol) was mixed with the blockers (2.25nmol) in 1:3 molar ratio to saturate all the available interaction sites between the blockers and the fixed regions of the DNA library. The sample was then denatured at 95°C for 5 minutes, then slowly brought to room temperature using a thermal cycler (MasterCycler Nexus Gradient, Eppendorf, Hamburg, Germany).

Selection phase

The beads, coupled with the resources, are mixed with the sample, prepared as described above. The capture of the DNAi is carried out at 40°C for 2 hours in stirring (600rpm, ThermoMixer, Eppendorf). Once the selection phase is completed, the sample is incubated for 2 minutes on a magnet and the supernatant is removed. The beads, that are now bound to the captured DNAi, are washed three times in SSC 1X buffer to eliminate the aspecific sequences. Finally, the beads are resuspended in water and incubated for 5 minutes at 60°C to recover the DNAi from the resources. The sample is quickly placed on the magnet and after 2 minutes the supernatant containing the selected DNAi is collected.

Amplification phase

The captured DNAi are then amplified by PCR with Q5® Hot Start High-Fidelity DNA Polymerase (New England Biolabs, Ipswich, MA, USA). DNA samples were diluted 10 times and 3μl were used for each PCR reaction. The PCR was performed in a final volume of 25μL using 0.25μL of polymerase, and each primer had a final primer concentration of 2.5μM. To allow the regeneration of the single strand library, the reverse primer was designed to carry a 5^′ biotin modification (see next paragraph for details). The thermal protocol was the following: i) denaturation step at 98 ? C for 2 minutes, ii) 28 cycles characterized by three thermal steps: 10 seconds at 98?C, 10 seconds at 68?C for Oligo1 (69?C for Oligo2) and 1 second at 72?C, iii) 2 minutes at 72?C. For each generation, 10 different PCR reactions were performed. The PCR products were then checked for size on a 2% TBE 1X agarose gel.

Regeneration phase

PCR products were purified by precipitation with 2.5 volumes of ethanol and 0.1 volumes of sodium acetate 2M, pH 5.2. After an overnight precipitation at -20?C, the pellet was washed with ethanol 75%, resuspended in water and quantified using the Qubit fluorometer (Invitrogen, Waltham, MA, USA). The single-strand regeneration was performed with streptavidin-coated magnetic beads (M-270 Dynabeads, Invitrogen, Carlsbad, CA) according to manufacturer’s instructions. Briefly, the PCR product was incubated for 15 minutes on a rotator with the magnetic beads. After 2 washes with binding and washing buffer (10mM Tris–HCl pH 7.5, 1mM EDTA, 0.2M NaCl) and a third one with SSC 1X buffer, an alkaline denaturation was performed by incubating the beads with 150mM NaOH for 10 minutes to induce the separation of the two DNA strands. This way it is possible to recover the unlabeled DNA strand in solution. The ssDNA is then collected and buffered with 1.25M acetic acid and TE 1X. The regenerated ssDNA is ready for the next round of selection and amplification.

Sequencing of the recovered products

To check the growth of DNAi species in our experimental model, some generations were sequenced through next-generation sequencing techniques. To prepare the sample for sequencing, the DNAi species captured as described in the “Selection phase” paragraph were first PCR-amplified using the same conditions discussed before. However, in this case, both forward and reverse primers were unlabelled. After PCR purification, performed as already described (see “regeneration phase” paragraph), the products were quantified and used to obtain libraries. About 300ng of DNA were used as starting material for the NEBNext® Ultra™ II DNA Library Prep Kit for Illumina® (New England Biolabs), and libraries were prepared following the manufacturer’s instruction. Sequencing was performed with the NextSeq™ 550 sequencer (Illumina, San Diego, CA, USA) and a paired-end strategy to obtain 75-nt long reads.

Replicas

The experimental workflow was repeated on two different sets of DNA libraries: Oligo1, whose results are described in the main text, and Oligo2, used as a replicate. The experimental conditions and procedures were the same for both sets of oligonucleotides, the only differences being the sequences of the fixed parts of the DNA libraries (and hence of the blockers and primers) and of the resources (Table 1). As a consequence, also the PCR conditions had to be adjusted, in terms of the annealing temperature that is one degree higher for Oligo2.

Eco-Evolutionary algorithm

To support experimental observation with a simple abstract model, we developed an Evolutionary Algorithm where a population of N_p sequences evolves in presence of N_r < N_p shorter sequences. In particular, N_r individuals in this population are selected at each cycle depending on their fitness, i.e. on their affinity to the resource, expressed via their ω with it. These survived sequences are then amplified by a factor of (roughly) N_p/N_r. In this work, two types of fitness functions have been explored: the first one is merely linearly proportional to the ω of each individual (i.e., where Σω is the sum of the ω values of the whole population), while the other one is a modification of the previous fitness which sets it to be beyond a threshold value ω_th. The code is written in C++, exploiting MPI Message Passing Interface Forum (2021) and Armadillo Sanderson and Curtin (2016, 2018) libraries for acceleration.

NUPACK calculations

We resorted to NUPACK for nucleotide sequences analysis, for the prediction of their free energies at equilibrium, either alone or when binding to one or two other oligomers. Schemes like those shown in Figure 4 have been obtained via the NUPACK web application Zadeh et al. (2011), while massive ΔG calculations have been performed with custom Python codes exploiting NU-PACK Python package (v4.0.0.27) Fornace et al. (2020); Dirks et al. (2007). The model is dna04 and ensemble=‘some-nupack3’; T = 40?C and [Na⁺] = 0.24 M, as in the experiments. Concentration of each species in the tube has been arbitrarily set to 10⁻⁶ M. Detection of hairpins and other secondary structures has been performed visually (thanks to the oxView software Bohlin et al. (2022); Poppleton et al. (2020)).

Acknowledgements

T.B. acknowledges support from MIUR-PRIN (Grant No. 2017Z55KCW). F.M. and S.S. acknowledge CloudVeneto (http://cloudveneto.it, Andreetto et al. (2019)) for the use of computing and storage facilities, through the SEDES and HPC-Physics projects.

Conflict of interest

The authors declare no conflict of interest.

Availability of data and materials

The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.

Code availability

The Jupyter Notebooks used for the present analysis are available at https://github.com/francescomambretti/stat_phys_synthetic_biodiversity.

Author Contributions

SS, EMP and TB conceived and designed the study; LC and EMP designed protocols and performed the experiments; FM, AT and SS analyzed the data and developed the theoretical models; LC and FM performed NUPACK analysis; LC, FM, SS and TB interpreted the data and wrote the manuscript.

Appendix 1 Data analysis protocol

Since we performed a paired-end sequencing, for each experimental cycle two FASTQ files were generated (R1 and R2).

The data shown in the paper all come from the cumulative analysis of R1 and R2 files. To analyze such files, we have developed an ad-hoc package (in Python and C++), currently maintained at https://github.com/francescomambretti/stat_phys_synthetic_biodiversity. Here follows a scheme of the main steps of the analysis:

1. Read FASTQ file, extracting the list of single-stranded DNA (ssDNA) sequences. Reverse sequences are reverted and complemented, so to have all the sequences oriented as 5^′ → 3^′. Fixed sequences at the ends are removed in this way: if the full blocker is present, then the 50 bases after it are retained. If some of the last bases of the blocker are lost (it may happen), then just remove the primer and keep the subsequent 50 bases. Slightly different choices for this cut of the sequences turned out to be substantially irrelevant for all the results.

2. Each sequence can be either saved or discarded, according to some filtering procedure. Criteria include the minimum quality of the read (≥ 10, measured according to the Phred score), the minimum and maximum allowed length, the presence of invalid nucleotides (indicated by ‘N’ in the FASTQ file) and other checks (see below in particular for PCR by-products).

3. The ‘survived’ sequences are processed, computing the related observables. Unique sequences are identified and their degeneracy is computed. This allows to monitor the time evolution of the abundance of each species (where the equivalence: DNA sequence-species holds).

4. For each of these unique ssDNAs, we compute the indicators quantifying the ‘affinity’ with the target sequence: in particular, the maximum consecutive overlap (MCO) and eventually the binding free energy with the resource and/or other oligomers.

5. Several post-processing options are available, such as: generating histograms of the ssDNA population based on their affinity with the resource; tracking the time evolution of dominant strands; determining the fraction of the total population covered by the n most abundant strands. All these quantities can be plotted by suitable Python scripts, provided in the GitHub repository.

By-Product Formation from PCR Amplification

During the amplification phase, described in our workflow, we observed the formation of PCR by-products. This phenomenon is known in SELEX systems, especially when an oligonucleotide library characterized by high heterogeneity is PCR-amplified Tolle et al. (2014). The generation of these by-products is caused by unspecific hybridization of the primers and/or of oligonucleotides, leading to PCR products abnormal in terms of length and distribution of random/fixed sequences. To limit this phenomenon we adopted blocker oligos with a phosphate group at 3’, so that they do not act as primers, and a stringent PCR protocol, characterized by: i) a high annealing temperature, ii) fewer PCR cycles, iii) a high primer concentration, iv) a very short elongation step (1s). Moreover, during the data analysis step, we developed an algorithm to identify and remove these spurious sequences.

PCR by-products are identified among the sequenced oligomers according to the following criteria:

• if either the sequence (after the cut of the blocker as previously described) contains at least n consecutive bases identical to n of the 5^′ blocker (i.e. an n-long MCO with a subsequence of the Primer1-F)

• or it has at least n consecutive bases identical to n of the 3^′ blocker then this sequence is considered as a by-product and excluded from the analysis.

Figure 2-figure supplement 1a shows how such unwanted sequences become ≈ 30% of the total sampled population at cycle 15 and then the striking majority (> 95%) in the last cycle, setting de facto an upper limit to the number of feasible experimental cycles. On the other hand, as we can see from Figure 2-figure supplement 1b, the behaviour of ⟨ω⟩ for the by-products is very similar to the one observed considering only regular oligomers, therefore their presence does not affect the main outcome of the experiment.

In Figure 2-figure supplement 2 the P (ω) for the first and last cycle of only the by-products of Oligo 1 are displayed and compared with the null model prediction. The main evolutive trend is also present in these sequences, which have however been discarded from all the analyses.

Appendix 2 Further analyses of the Eco-Evolutionary Model

We here show some further results and sensitivity analyses of the Eco-Evolutionary stochastic individual based model. In Figure 2-figure supplement 5 we show how the eco-evolutionary model can describe the emergent order and complexity, that is a loss of the initial huge diversity of the random strings towards few dominating species, during evolution. This behaviour can be quantified 1) by the ratio CR between the zipped file size Benedetto et al. (2002) and the original file size of all the sequences oligomers in the ecosystem or 2) by computing the Shannon entropy of the Relative Species Abundance Hill (1973). We can see that that the model is able to reproduce, at least qualitatively, the marked decrease of diversity observed during evolution.

In Figure 2-figure supplement 6a we show the outcome of the model if we do not include the saturation at ω_th = 10. As we can see, the effect is evident as after a few cycles the experimental ⟨ ω⟩ lies significantly under the model (which promotes too much the strongest sequences). Similarly, Figure 2-figure supplement 6b shows the evolution for the average MCO in the case we use only γ = 1 in the simulation (here, ω_sat = 11). It is very apparent how in this case the simulated initial evolution is much slower with respect to the observed from the data, while the final trend is captured. These results indicate that both saturation and the non-linear fitness time dependent fitness function are fundamental ingredients to describe the empirical evolution of ⟨ ω⟩. Note that in both cases of the above picture, there is a 10% random sampling drift.

Finally, Figure 2-figure supplement 7a and Figure 2-figure supplement 7b show the simulations results for different total number of individuals (N = 10⁶, 10⁵ and 10⁴) and resources (R = N/100, keeping this ratio fixed) and for different initial conditions (i.e., different sequences populations), respectively. We can see that the qualitative results concerning ⟨ω⟩ as a function of time are robust, beyond a given system’s size (the silver data in Figure 2-figure supplement 7a are already of a system’s size large enough, while the blue ones correspond to a too small system). Considering a different initial (random) population of sequences, keeping fixed all the other hyperparameters (in particular N = 10⁶ and R = 1064), as shown by Figure 2-figure supplement 7b preserves the qualitative agreement with the experimental data, assessing the substantial independence on the initial conditions. Overall, these further analyses show the robustness of the qualitative behavior of results of the Eco-Evolutionary model, pointing out to two relevant phases of the evolution of the oligomers, as highlighted in the main text.

Further results on the interactions between oligomers

In this section we highlight some further results on the interactions between species and resources and among oligemers, respectively.

Figure 4-figure supplement 1 shows the average pair-wise interactions between 10 random strings (as representative of the whole population) during the different evolutionary cycles (t = 0, 9, 18, 24), ⟨ΔG_pp⟩. Such averages have been computed by giving to NUPACK a pool of 1000 sequences, randomly subdivided in 100 groups of 10 sequences each. We can see that, during evolution, the selected oligomers increase their average interaction strength, a signature that the evolution favors not only strings that have higher overlap with the resource, but also that interact in some way with the other species.

Figure 4-figure supplement 2 is the 2D probability distribution p(ω, x_R), computed from experimental data, to simultaneously have a given MCO= ω and x_R, the rightmost basis of the of the resource involved in the formation of such an MCO. Noticeably, we start from a distribution peaked in ω = 4 and 10 ≤ x_R ≤ 12, essentially corresponding to purely random binding. Across cycles, the probability distribution is deeply modified, with a shift towards the bottom right corner. Along the horizontal axis this simply reflects the shift in the average ⟨ ω⟩ of the population, and also the fact that winners have large overlaps. Interestingly, we observe a significant shift of x_R towards high values, meaning the MCO are preferably formed far away from the bead surface, near the free end of target strands.

Results of the experimental replica

In Figure 2-figure supplement 3 and Figure 2-figure supplement 4 show the main outcomes of the replica experiment Oligo2. The evolution of the average MCO ⟨ ω ⟩over cycles is qualitatively the same of that one found in Oligo1: we have a initial steep grow ⟨ ω ⟩, compatible with a non linear fitness, and then we observe a slower increase for the last cycles, saturating for ⟨ ω ⟩ ≈ 7. The histogram shows the whole p(ω) at cycles t = 0 and t = 18. As for what we found in Oligo1, at t = 0 the behaviour of p(ω) is compatible with our null random model, while at t = 18 the oligomers with higher ω are selected, and the tail of p(ω) is fatter than the one expected by chance.

These results show that indeed the main outcome of Oligo1 are replicated by Oligo2.

Appendix 3 A combinatorial null model with no evolution

As a first step to understand what happens in our experiments, we have built a null model, assuming purely random attachment sites, i.e. there is not any search for the relative positions of binding filaments so to maximize the MCO. In this way we have a null expectation for the typical MCO distribution in case of oligomers insensitive to the presence of target sequences.

Let us call a the number of consecutive bonded pairs between oligomers 1 (e.g. the species) and 2 (e.g. the resource), with and R_1,2 being the relative position of the left ends of the two filaments, as also explained in Mambretti et al. (2022b) and in the main text. Consider the initial population being uniformly random selected (in the high 4^L dimensional space of all possible DNA strands). The first quantity we are interested to calculate is π(a), the probability density of a at time t = 0, i.e. before the evolution starts. This distribution is also the distribution of the subsample (i.e, the ≈ 10⁶ strands which we are able to sequence at each cycle), if we assume purely random attachment (no competition among the consumers, no optimization of the attachment site).

The relative position of the two strands is important since the target can be either “inside” (Figure 3-figure supplement 1, panel a)) or surpassing one of the borders (Figure 3-figure supplement 1, panel b) and C)) of the longer oligomer after the attachment. In the picture, R_1,2 = 0 when the first nucleotide of the target overlaps with the first one of the consumer. From now on, we will refer to R_1,2 simply as R.

Our goal consists in finding the MCO of L-long strand (the species) with a l-long target ssDNA (the resource), where l < L. In the following calculations, we therefore assume that R is a random uniform variable (−l + 1 ≤ r ≤ L − 1) and we will calculate a| _R (for this R).

Such a problem is sketched in Figure 3-figure supplement 2 (where R = 2, but the following holds for any R in the allowed range), where we indicate with d the number of mismatching bases between the target (red) and the consumer (blue), with u the number of matching ones, being d + u := l. Labeling as x_i the number of consecutive matching nucleotides comprised between two consecutive mismatching bases, as shown in the picture, the problem can be reformulated in terms of counting how many integer solutions can be found for the equation:

that corresponds to determine the number of ways in which a set of x_i solving Eq. 1 can be obtained. Note that x_i vanishes between two neighboring pairs of mismatching sites, as marked by the arrows in the picture. In the case reported in the picture, e.g., x_i = (1, 0, 1, 4, 0, 4, 0, 0, 2) and and d = 8.

Ref. Charalambides (2008) provides us with a formula to answer this question. We stress that finding the number of solutions A_d,u for the above equation is equivalent to obtaining the probability distribution π(a) for the relative occurrence of each a in a sample of DNA strands, with uniformly random extracted relative positions of the two sequences.

In our case the target string may not go beyond the consumer species boundaries (as the latter is longer than the former). But we can modify the name of some variables to adapt the theorem to our problem and make some constraints explicit:

• from i = 1, …, n it follows that n = d + 1 because we have d + 1 terms (from 0 to d).

• s_i = 0 ∀i, because each x_i ≥ 0 by definition. Therefore, s = 0.

• at least one of the x_i must be equal to a

• in principle, m_i is equal to a for each i (which means that there are no consecutive overlaps larger than a).

• it must be that a ≤ u ≤ min((d + 1)a, l) (i.e. the total overlap u - which, at most is equal to d + 1 times a - must not be larger than l)

In fact, by using that s = 0 and d + u = l, the theorem Eq. 2 can be rewritten as:

where the inner sum can be replaced by(how many ways are there to dispose the d +1 integer indexes to form r groups):

First of all, can be placed inside the summation because it corresponds to the = 0 case. Now, we have to apply the constraint that at least one of the w_j = a, which means having the first term where the w_j are always equal to a, minus the second term which includes all the cases where all the w_j < a:

Note that the previous expressions concern only the cases where 0 ≤ R ≤ L − l. However, this should not matter for practical purposes because R does not appear in the above expressions.

Similar formulas can be found when the target surpasses the left/right border:

if − l + 1 ≤ R ≤ −1

while if L − l ≤ R ≤ L – 1

Using all the previous results, the general equation for the number of ways which lead to a MCO a is given by:

where the 3 options are originated by the 3 options for the relative position of the target and of the consumer sketched in Figure 3 - figure supplement 1.

Let us now derive n_k, ∀k. In particular, if the target “resource” is completely within the “consumer” species boundaries, then we have

with (L− l + 1) is the total number of possible positions R where the consumer can attach to the resource, 4^L is the number of ways to build a L-long string (i.e. all possible species), 3^d represents the number of arrangements for the non-matching bases (3 is because we exclude the only matching base) and

Analogously,

where the summations run over the admitted integer R values as well as over the allowed d values. In this case, 4^L−l−R is the number of ways in which the nucleotides not involved in the attachment can be obtained. Here.A similar equation holds for n₃(a), where

and .

To get the distribution of the overlap measure π(a), it suffices to divide n(a, R) by the total number of ways in which we can build a string, i.e. 4^L times the number of possible positions R where the attachment can happen, i.e. L + l − 1.

The underlying probability distribution π(a) is naturally discretized and can be represented as an histogram, with the occupancy probability for each bin a (the affinity to the resource). π(a) is represented in Figure 3-figure supplement 3a as orange bars; grey bars in the same panel represent the average of 50 independent simulations with K = 10⁶ 50-mers uniformly generated, where their overlap with the target strand - once uniformly random extracted their relative position - is measured.

The analytical result for the null model and simulations yield almost exactly the same results, with the difference that for large MCO we have strong sampling effect in the simulations. In order to understand how possible biases in the initial frequency of each nucleotide in the starting population, we run N = 50 independent simulations with a non-uniform frequency for A,C,G,T in the sequences. For instance, the panel b) of Figure 3-figure supplement 3 compares the π(a) obtained with A,C,G,T having a probability of 25% (grey data) with the corresponding distribution obtained having: A=23%, C=19%, G=29%, T=29%. The average bin occupancy does not change significantly, but only the large a tails are slightly affected by the bias and the error-bars are much larger for the majority of the a values.

We furthermore set a threshold T to model the physical constraint that if a < T, the two ssDNAs cannot bind. We thus obtain a null model physically constrained distribution π^′(a), where π^′(a) = 0 if a < T, while is the same as before (after a proper normalization) for a ≥ T. The new analytical distribution π^′(a) reads as

where sgn(0) := 0. The corresponding pseudo-algorithm to generate random number from such distribution is

• sample a random number z from π(a)

• if z < T, sample another number from π(a)

• else, save it in the new distribution π^′(a)

The comparison between analytical and numerical results are reported in Figure 3-figure supplement 4a and Figure 3-figure supplement 4b (linear and logarithmic scale on y axis, respectively). The cyan bars correspond to the analytical distribution for the null model with physical constraint (π^′(a)), whereas the purple bars represented the average over N = 200 independent simulations of π^′(a) performed via Metropolis Monte Carlo accept-reject technique Metropolis et al. (1953). A custom Python script which simulates this process is available at the GitHub repository associated to this work. The null model that we present in the main text is the physical constraint one, i.e. π^′(a).