Abstract
Aptamers are single-stranded oligonucleotides selected by evolutionary approaches from massive libraries with significant potential for specific molecular recognition in diagnostics and therapeutics. A complete empirical characterisation of an aptamer selection experiment is not feasible due to the vast complexity of aptamer selection. Simulation of aptamer selection has been used to characterise and optimise the selection process; however, the absence of a good model for aptamer-target binding limits this field of study. Here, we generate theoretical fitness landscapes which appear to more accurately represent aptamer-target binding. The method used to generate these landscapes, selective phenome growth, is a new approach in which phenotypic contributors are added to a genotype/phenotype interaction map sequentially in such a way so as to increase the fitness of a selected fit sequence. In this way, a landscape is built around the selected fittest sequences. Comparison to empirical aptamer microarray data shows that our theoretical fitness landscapes more accurately represent aptamer ligand binding than other theoretical models. These improved fitness landscapes have potential for the computational analysis and optimisation of other complex systems.
1. Introduction
1.1. Background
Aptamers are single-stranded nucleic acid sequences capable of specific high-affinity binding [1–4]. This makes them attractive candidates as recognition molecules in diagnostics and therapeutics. Aptamers are isolated by systematic evolution of ligands by conventional exponential enrichment (SELEX), which involves several iterative steps of incubation with target, washing away of weak binders, and PCR amplification of strong binders.
Aptamer selection is complex. Many variables such as library size, quantity of target, temperature, selection buffer, pH, degree of PCR amplification, and use of mutation or recombination diversification need to be considered. Due in part to these factors, less than 30% of classical SELEX experiments are successful in isolating aptamers with dissociation constants less than 30 nM [5]. Understanding the dynamics of the selection process is extremely important; but what does this entail? The DNA required to fully represent the number of permutations in a 75-base aptamer library would be roughly equal to the mass of the moon [6]. In order to represent this immense sequence space, an initial SELEX library may contain up to 1015 molecules. A rigorous empirical analysis of anything close to this number of library members is simply not feasible.
Nevertheless, empirical analyses of smaller fractions of a DNA aptamer libraries have been undertaken. The two main empirical selection analysis techniques used are high density DNA microarrays and high-throughput sequencing (HTS). Briefly, high density microarrays can contain up to approximately 1 million features, each representing an aptamer in a library. After array incubation with fluorescent target and a washing step, the binding affinity score of all aptamers on the array can be measured by fluorescence scanning. Platt et al. [7], Knight et al. [8], and Rowe et al. [9] used microarrays to both evolve aptamers and gain insight into an aptamer binding landscape. Additionally, DNA microarray data has been applied to aptamer specificity landscapes [10], fitness landscape morphology [11], and aptamer affinity maturation [12]. In comparison, the possible sequence space coverage using HTS is much greater, with Illumina’s HiSeq HTS capable of yielding sequence data for more than 70 million sequences from a single lane [13]. Using this approach the copy number of a sequence is used as a proxy for its target binding strength so that the fitness of individuals in the library pools can be estimated. Cho et al. used this HTS approach to monitor microfluidic aptamer selection rounds and gauge enrichment [14]. PCR bias may distort this copy number/binding correlation but, by using a motif based statistical framework such as MPBind [15], the binding potential of aptamers can be predicted, eliminating error from PCR bias. Although both DNA microarrays and HTS led to major breakthroughs in understanding library sequence space fitness and selection, these techniques are only capable of analysing a small fraction of a given library’s sequence space. Another approach to analysing aptamer selection is via computational simulation. The challenge in simulating aptamer selection is the design of a suitable model for aptamer binding fitness.
Computational approaches to model aptamer fitness by virtue of folding include secondary structure prediction by minimum free energy [16] and inverse folding [17]. These approaches can be computationally expensive and while being excellent models for folding they may not capture the higher complexity of molecular binding. Hoinka et al. coded a program to simulate the aptamer selection process called “AptaSim” [18]. The binding model used chose aptamer affinities at random without relevance of sequence. While AptaSim was an important step forward in simulating selection enrichment and mutation copy number, AptaSim cannot appropriately represent heritability or represent binding correlation between related sequences required for the study of genetic systems. Oh et al. used a string matching function as a binding fitness model to simulate aptamer selection [19]. This model does include heritability and binding correlation between related sequences, but as only close range epistasis is possible by using one “optimal solution” aptamer the landscape is cone shaped and would often be unrepresentative of a true aptamer binding landscape.
Kauffman’s model is a robust mathematical model, related to study of autocatalytic sets, which serves as an objective function relating genotypic sequence to phenotypic fitness [20]. Derivations of the original model have been used to describe complex interacting systems in areas as diverse as immunology [21], evolutionary biology [22], and economics [23]. The model describes a fitness landscape whose size is determined by the number of components in system () and the ruggedness of the landscape can be tuned using the degree of interaction of these components (). This system is perhaps best described when used to represent problems in evolutionary biology as originally intended by Kauffman. A population of genomes where each contains N genes are given fitness values based on the sum of fitness contributions from each of their genes. The fitness contribution of each gene is determined by its interaction with other genes within its own genome. The interacting genes can be positioned sequentially, randomly, or by some other gene interaction pattern predetermined by an interaction map. The allelic sequence of these interacting genes is designated a fitness contribution, usually from a generated random distribution. In this way, the allelic substitution of one interacting gene means there will be a completely different fitness contribution score for the entire collection of interacting genes. In the model by increasing , the number of interactions between genes, the complexity of the system, and the ruggedness of the landscape are increased. In addition to the value of , the position of these interactions on the interaction map is of great importance to the fitness landscape.
Typically the model is used as a scoring system for a population of genomes which can evolve via diversifications such as mutation or recombination. In this way the model is an objective function for a complex system. As mentioned earlier the model can be adapted to many other areas of study. In this paper we use the model to represent binding of an aptamer to a target analyte. In this representation would be equal to length of the aptamer in the library and would be equal to the interactions of bases within each aptamer. Many modifications to the original model have been made, some to optimise the model for a given research area. Herein we will describe some modifications to the model which are aimed at optimising the model to represent binding of an aptamer to a target analyte. The model was believed to resemble molecular fitness landscapes similar to the binding of an aptamer to an analyte [24]. In the model mutational additivity usually holds for noninteracting positions in sequences. This mutational additivity is biologically accurate as has been demonstrated for several proteins [25–33].
Wedge et al. used an model for the simulation of protein directed evolution (DE) [36], a similar field to aptamer selection. Binary strings of length 40 and 100 were used with random epistatic interactions varying from = 0 to 10. Genetic algorithms utilising mutation, crossover, different library sizes, and selection pressures were simulated and compared to deduce general rules for protein directed evolution, which are of great use to DE experiments. As noted in this study, the “No Free Lunch Theorem” (NFL) [37] establishes that all search algorithms perform exactly the same when averaged over all possible problems. This infers that, for an optimisation algorithm, any elevated performance in one class of problem is exactly paid for in performance in another class. If there is discrepancy between a real life system and a model used to describe it, any elevated performance in optimisation using simulation of the model is exactly paid for in performance for the real life system. This illustrates the need for an accurate model when using simulation results to improve empirical ligand selection experiments.
Despite this biological accuracy in regard to mutational additivity, the classical model may have limitations in representing some biological systems. The model’s greatest utility is that ruggedness can be tuned using the epistasis variable . However, this epistasis is quite uniform throughout the sequence. For some biological applications, a higher amount of epistasis is desirable. As increases the landscape tends to become more multipeaked and rugged, to the point where it is too chaotic to allow adaptation. Kauffman refers to this phenomena as the “complexity catastrophe” [38]. Kauffman goes further to say “the complexity catastrophe is averted in the model for those landscapes which are sufficiently smooth to retain high optima as increases” [38]. Thinking along these lines provided a solution to the complexity catastrophe, creating complex landscapes which retained a smoother surface.
1.2. Constructional Selection of Landscapes
Altenberg developed an evolutionary approach to selecting epistatic interaction, thereby creating landscapes which were smoother than classic landscapes with the same degree of epistatic interaction [34]. Altenberg achieved this using selective genome growth, a type of constructional selection, to create modular interaction matrices. These selected matrices have reduced epistasis which give rise to smoother fitness landscapes [34, 39]. Selective genome growth is a process by which the genome of the fittest individual is expanded one gene at a time (Figure S1a in Supplementary Material available online at https://doi.org/10.1155/2017/6760852). The new gene is only kept if the fitness of a selected optimum genome is increased. In this way the probable global optima of the landscape is constructed and all other points on the landscape are relative to this optimisation. A similar method for creating landscapes was devised by Hebbron et al., which uses preferential attachment growth process to add genes to a genome [40]. A problem with these two approaches, when applied to specific applications, is that due to the increasing returns of the selection system these methods attribute extremely high pleiotropy to a handful of genes (vertical lines in Figure 2(c)). This phenomenon of increasing returns of gene control is biologically appropriate and accurate for a system describing a group of genomes, but when describing the binding of an aptamer to an analyte this high aggregated pleiotropy is not biologically appropriate. Each base in an aptamer has a relatively low number of interactions due to its spatial capacity, meaning that high aggregated pleiotropy is not biologically representative for an aptamer.
Herein we have created a new model that we have termed “selective phenome growth.” Selective phenome growth is a constructional selection technique in which phenotypic contributing factors are added to a genotype-phenotype interaction map incrementally (Figure S1b) in such a way that each new phenotypic contributing factor increases the fitness of global or local optima. Additionally, comparison is made between selective phenome growth landscapes and aptamer binding landscapes.
2. Model and Methods
2.1. Selective Phenome Growth to Create a Genotype/Phenome Interaction Map
Selective phenome growth is a new method of constructing an interaction matrix one phenotypic contributor at a time. The method or representing the interaction map is the same as Altenberg’s [34], with slight modification to represent aptamers, and is as follows:(1)The aptamer consists of binary valued bases that have influence over phenotypic functions, each of which contributes a component to the total fitness.(2)Each base controls a subset of the fitness components, and, in turn, each fitness component is controlled by a subset of the bases. This genotype-phenotype map can be represented by a matrix, of indices , where indicates that base affects fitness component .(3)The columns of , or epistasis vectors, , give the bases controlling each fitness component .(4)The rows of , or pleiotropy vectors, , give the fitness components controlled by each base .(5)If any of the bases controlling a given fitness component is altered, the new value of the fitness component will be uncorrelated with the old one. Each fitness component is a uniform pseudo-random function of the genotype, .(6)If a fitness component is affected by no genes, it is assumed to be zero.(7)The total fitness is the normalized sum of the fitness components:Similarly to Altenberg’s selective genome growth [34], a test sequence of length is randomly generated. In the phenotype selection loop (Figure 1), positions from a total of are selected at random and added as a phenotypic contributor to the interaction matrix. If the new addition decreases the overall fitness of the test sequence it is removed; if the new addition increases or does not change the overall fitness of the test sequence it is kept in the interaction matrix. The selection loop is usually repeated until phenotypic contributors are selected.
(a) Local = 6
(b) Random = 6
(c) Genome selected
(d) Novel phenome selected
Selective phenome growth can be seen as randomly selecting one sequence as the fittest member of a prospective landscape and evolving the interaction map and therefore the fitness landscape around this fittest sequence. In this way phenotypic contributors are sequentially added in such a way that each increases the fitness of the selected fittest sequence.
2.2. Characterising Landscapes
Adaptive walks were used to find the two measures, accessibility of local optima and hamming distance from fittest optima similarly to Kauffman 1993 [38]. 100,000 adaptive walks from randomly selected points on each landscape were performed and the frequency each local optimum reached was recorded. The more adaptive walks terminating at an optimum, the larger the basin of attraction. In this way the accessibility of a local optimum is a measure of its basin of attraction. For the same 100,000 adaptive walks the hamming distance from the fittest optimum found for all other optima found was calculated. Hamming distance from fittest optima is a measure of peak clustering.
Mean path divergence (MPD) is a measure developed by Lobkovsky et al., 2013, which quantifies the degree of similarity among monotonic evolutionary trajectories with the same endpoints [11]. MPD was used to gauge the roughness of a given landscape.
2.3. Breaking the One-Gene One-Phenotypic Contributor Paradigm
Previously in model literature the number of phenotypic contributors has been equal to the number of genes. This rigid framework is not representative of biological systems as the number of phenotypes may be lower or may exceed the number of genes in a genome. This one phenotype per gene paradigm made sense with earlier interpretations of the model with local interactions, as each numbered row in an interaction matrix corresponded to a phenotypic contributor determined by its gene interacting with other genes. This one phenotype per gene paradigm where each locus always contributes to its own fitness contribution can be seen as the red diagonal lines in Figures 2(a) and 2(b).
Altenberg’s genome growth constructional selected interaction maps do not have this pattern of each locus always contributing to its own fitness contribution as all interacting genes for a given phenotypic contributor are selected at random. In Altenberg’s genome growth interaction maps the number of phenotypic contributors is usually predetermined and equal to the number of genes [34]; however this is not necessarily the case for our phenome selected interaction maps.
In the case of our phenome selected interaction maps the phenotypic contributors are added sequentially and so their number can be varied. This in effect varies the complexity of the landscape in such a way that the landscape can be tuned by adding or subtracting phenotypic contributors. The effect of this landscape tuning via adding or subtracting phenotypic contributors was analysed in Results and Discussion (Figure 6).
3. Results and Discussion
3.1. Comparison of Epistatic Maps
From Figure 2 the different types of interaction map can be compared. For locally distributed interactions (Figure 2(a)) the rows, or epistasis vectors, are equal to the columns, or pleiotropy vectors, which are equal to the value plus one. For randomly distributed interactions (Figure 2(b)) the epistasis vectors are equal to the value plus one and the pleiotropy vectors are a Poisson distribution with a mean of .
For genome selected interactions (Figure 2(c)) the epistasis vectors are equal to the selected values plus one which yields a mildly skewed distribution (Fisher-Pearson coefficient of 0.13). The pleiotropy vectors for genome selected interactions are heavily skewed (Fisher-Pearson coefficient of 1.90) towards a handful of loci which can be seen as the vertical lines in Figure 2(c). For phenome selected interactions (Figure 2(d)) the epistasis vectors are equal to the selected value plus one which yields a skewed distribution (Fisher-Pearson coefficient of 0.69). The pleiotropy vectors for phenome selected are a normal distribution (Fisher-Pearson coefficient of 0.07). From Figure 2 it can be seen that phenome selected growth produces an interaction map with less aggregated pleiotropy when compared to genome selected maps.
In order to make a relevant comparison between genome selected and phenome selected landscapes the complexity of phenome selected growth interaction map was limited to be equal to that of the genome selected landscape. This was achieved by limiting the number of phenotypic contributors such that the total number of base interactions was limited to that of a previously generated genome selected interaction map. As the number of interactions a phenotype has is randomly selected, the actual number of interacting bases selected for phenome selected growth is greater than the enforced limit relating to the genome selected interaction map. Consequently, as the genome selection map was involved in 577 interactions, the phenome selected map was involved in 588 interactions.
3.2. Comparison of Theoretical Landscapes
The accessibility of local optima was used as a comparison for different landscapes. Explained briefly, 100,000 adaptive walks from randomly selected points on each landscape were performed and the number of times each local optimum was reached was recorded. The larger the basin of attraction of a local optimum is, the more adaptive walks will terminate at that optimum. In this way the accessibility of local optimum is a measure of an optimum’s basin of attraction.
From Figure 3 the comparison of the accessibility of local optima demonstrates the difference in basin size between different types of binding landscapes. For low value landscapes such as (Figure 3(a)), the model exhibits a range of basin sizes, with an optimum’s basin size loosely correlating to its fitness. At higher value (Figure 3(b)) the correlation between an optimum’s basin size and its fitness is abolished. The basin size was reduced sufficiently that during the adaptive walk simulation only a single adaptive walker can reach each local optimum.
(a) Random
(b) Random
(c) Genome selected growth
(d) Phenome selected growth
The landscape derived from selective genome growth (Figure 3(c)) exhibits a wider range of basin sizes than a random landscape (Figure 3(b)) with a maximum walker tally of 9 and a mean walker tally of 1.006 with a standard deviation of 0.113. In this landscape an optimum’s basin size correlates to its fitness.
The landscape derived from our selective phenome growth (Figure 3(d)) exhibits a much wider range of basin sizes than the genome selected landscape (Figure 3(c)) with a maximum walker tally of 407 and a mean walker tally of 1.010 with a standard deviation of 1.735. These statistical results indicate that the phenome selected landscape has a larger range of basin sizes when compared to genome selected landscapes. Furthermore, the distribution of these basin sizes has more variation in phenome selected landscapes when compared to genome selected landscapes.
In the phenome selected landscape there is a stronger correlation between an optimum’s basin size and its fitness when compared to genome selected landscapes. Additionally there are a lower number of unique optima for phenome selected (98977) compared to genome selected (99381) landscapes. This lower number of unique optima found for phenome selected landscape indicates a lower number of basins corresponding to a smoother, less chaotic landscape.
One major criticism of landscapes when used as fitness landscapes is that with larger values they become excessively chaotic and an unrealistic representation of reality. This observed smoothness of the genome and phenome selected landscapes despite the landscape complexity suggests that the model is more similar to that of a real binding landscape.
From Figure 4 the correlation between an optimum’s fitness and its similarity to the fittest optima found can be seen for different landscapes. Typically hamming distance from fittest optima plots only includes other points from optima reached during adaptive walks, shown as red in Figure 4 [35, 40], but here we have included the two additional groups in order to better represent the landscapes. The first is a random sequence group, shown as yellow in Figure 4, which gives some indication of the average sequence within the landscape. The second is a minima group, shown as blue in Figure 4, which gives some indication of the troughs or valleys in the landscape.
(a) Local
(b) Local
(c) Genome selected growth
(d) Novel phenome selected growth
For the landscape (Figure 4(a)) the optima group, displayed as red points, show an inverse correlation between fitness of an optimum and its distance from the fittest optima. This correlation is indicative of the “Massif Central” global structure described by Kauffman [35]. The random group, displayed as yellow points, show a range of fitness values and no discernable relationship between fitness of a random sequence and its distance from the fittest optima. The minima group, displayed as blue points, are all of a similar fitness which indicates a convergence on a small number of related low fitness sequences. This may be an artefact of the low complexity of the landscape.
For the higher complexity landscape (Figure 4(b)) the optima group, displayed as red points, show no correlation between fitness of an optimum and its distance from the fittest optima. The absence of this correlation shows a disruption of the “Massif Central” type landscape observed with higher values.
In such a landscape, optima’s position in sequence space occurs more randomly in a chaotic fashion. The random group displayed as yellow and the minima group displayed as blue show sequences with the same similarity of the fittest optima as the optima group. This is indicative of a chaotic landscape.
For the genome selected landscape (Figure 4(c)) the optima group, displayed as red points, show an inverse relation between fitness of an optimum and its distance from the fittest optima. The line of best fit for this relation shows a decreasing exponential-like trend with a decreasing rate of fitness increase with closer hamming distance to the fittest optima found. This seems to represent a landscape in which many similar fitness optima exist at a relatively distant hamming proximity from the fittest optima found, indicative of a plateaued “Massif Central” fitness landscape.
The genome selected random group displayed as yellow points show quite an even distribution with a slight skew towards the inverse relation between fitness of a random point and its distance from the fittest optima. The genome selected minima group displayed as blue points show a relation between fitness of an optimum and its distance from the fittest optima. This relation is the exact opposite to that the optima group displayed. Intriguingly, this indicates that lowest minima are found closer to the fittest optima than the average minima. This pattern was observed on multiple genome selected interaction maps, so it is not an artefact.
One explanation for this effect is that for genome selected interaction maps the high pleiotropy attributed to a few loci means that there exist sequences with great sequence homology to the fittest optima but with relatively low fitness as these high pleiotropy and high fitness contributing loci are not optimised. To test this the optimum sequence was mutated at each of the high pleiotropy loci according to the interaction map and it was found that just 9 mutations yielded a sequence with a fitness score of 0.35. This sequence was not a minimum as it was not represented by one of the blue minima points in Figure 4(c). This example shows that it is useful to investigate minima and not just optima when characterising a landscape.
For the phenome selected landscape in Figure 4(d) the optima group, displayed as red points, show an inverse relation between fitness of an optimum and its distance from the fittest optima. The line of best fit for this correlation shows an increasing exponential-like trend with a higher rate of fitness increase with closer hamming distance to the fittest optima found. This seems to represent a landscape in which fitness level drops steeply from the fittest optima and gradually levels out. This phenome selected steeple type “Massif Central” landscape shows the opposite effect of the plateaued “Massif Central” described for the genome selected landscape. As mutational additivity is observed in empirical landscapes [25–33], the most realistic theoretical aptamer binding landscapes are Mount Fuji-like in nature [42] with fewer distinct fittest sequences. The steeple nature of the phenome selected landscape fits this Mount Fuji-like approximation well. Furthermore mutational additivity, or gradient of the red lines of best fit in Figure 4, holds for phenome selected landscape sequences close to the fittest sequence (Figure 4(d)); however, for genome selected landscapes mutational additivity for sequences close to the fittest sequence does not hold (Figure 4(c)). This difference indicates that phenome selected landscapes are a more accurate representation of aptamer binding fitness landscapes.
The phenome selected random group displayed as yellow points show a slight skew towards the inverse relation between fitness of an optimum and its distance from the fittest optima, similar to that of the optima group. The phenome selected minima group displayed as blue points show an even distribution with no discernable correlations.
3.3. Optimising Phenome Selected Landscapes
As described earlier, phenome selected landscapes have the unique capability to alter the number of phenotypic contributors thereby dismissing the one-locus one-phenotype contributor paradigm. By changing the number of phenotypic contributors in an interaction map another mechanism, aside from altering value, can be used to tune the landscape. Similarly to using value to tune a landscape, the number of phenotypic contributors can be used to tune the landscape from states of relative order to states of relative chaos. In Figure 5(a) four interaction maps with varying numbers of phenotypic contributors were used to create fitness landscapes. These landscapes were assessed in two ways, firstly by their accessibility of optima (Figure 5(b)) and secondly by their optima’s hamming distance from the fittest optima found (Figure 5(c)).
(a)
(b)
(c)
For a low number of phenotypic contributors (Figure 5(a)(i)) local optima are inaccessible (Figure 5(b)(i)). At such low numbers of phenotypic contributors there exist loci in the interaction map which are not represented in any phenotypic contributors. This means that they have no bearing on sequence fitness. Mutational relatives in which these loci are altered are, from a phenotypic scoring point of view, equal but, from a sequence point of view, different. Therefore convergence on a single sequence representing a scoring optima is more difficult. Additionally, for a low number of phenotypic contributors, the optima’s hamming distance from the fittest optima (Figure 5(c)(i)) shows an inverse relation between fitness of an optimum and its distance from the fittest optima. However, despite this relation there is a relatively large fitness and hamming distance gap between the fittest optima and the next fittest optima. This gap reflects a large evolutionary jump from the optima to the fittest optima found.
For a low to intermediate number of phenotypic contributors (Figure 5(a)(ii)) the local optima (Figure 5(b)(ii)) are accessible with a reasonably large basin size for the fittest individuals and a range of basin sizes which correlate to fitness. This observed reasonably large basin size and basin size correlation to fitness are indicative of a well-tuned landscape. Additionally, the hamming distance from fittest optima (Figure 5(c)(ii)) shows a strong inverse relation between fitness of an optimum and its distance from the fittest optima. The distance between the fittest optima and the other optima is minimal indicating good mutational additivity and a well-tuned landscape.
For the intermediate to high range of phenotypic contributors (Figure 5(b)(iii)) the local optima (Figure 5(b)(iii)) are accessible albeit with a skewed distribution of basin sizes. The hamming distance from fittest optima (Figure 5(c)(iii)) shows an inverse relation between fitness of an optimum and its distance from the fittest optima. There is a relatively large fitness and hamming distance gap between the fittest optima and the next fittest optima. This gap indicates overoptimisation towards the lead aptamer sequence, with the landscape becoming more rugged and optima being of a lower score.
For the high number of phenotypic contributors (Figure 5(a)(iv)) only one local optimum (Figure 5(b)(iv)) is accessible and its basin size is extremely large. This is indicative of landscape that has been overoptimised for the lead aptamer sequence to the point that is has converged on a singularity. The hamming distance from fittest optima (Figure 5(c)(iv)) shows an inverse relation between fitness of an optimum and its distance from the fittest optima. However there is a larger fitness and hamming distance gap between the fittest optima and the next fittest optima. This gap indicates severe overoptimisation towards the lead aptamer sequence, with the landscape becoming increasingly rugged and the optima being of a lower score.
This progression through the number of phenotypic contributors depicted in Figure 5 shows that there is an optimal number of phenotypic contributors when designing an interaction map for a fitness landscape.
3.4. Comparison to an Empirical Binding Landscape
To compare theoretical binding models to empirical aptamer-protein binding, DNA microarray data was used. Previously, Rowe et al. analysed a complete DNA-protein affinity landscape using DNA microarrays [41]. In this study all possible DNA oligomer variants of 10 bases were synthesised onto a microarray. The microarray was incubated with fluorescent allophycocyanin protein and fluorescent scanning performed to reveal the entire protein-binding landscape [41]. We used this complete empirical binding dataset for comparison to our theoretical landscapes.
To compare theoretical and empirical aptamer binding data we used mean path divergence of monotonic trajectories, a method of measuring smoothness of a landscape [11]. Briefly described, a monotonic function is entirely nonincreasing or nondecreasing; ergo a monotonic trajectory starts at point A and finishes at point B without changing the polarity of its slope. Path divergence is calculated by averaging the deviation within a set of monotonic trajectories between two set points, one point being a local maximum. Mean path divergence (MPD) is the average of path divergence values between all possible points on a landscape [11]. MPD is defined by Lobkovsky et al. as where is the probability the occurrence of and is the distance between the trajectories and [11]. The current selected local maximum is . In this way MPD is a measure of the smoothness of landscape, one important feature for which there may be disparity between a theoretical model and empirical data. Each individual path divergence measurement is binned according to its hamming distance between the two given points before averaging so a mean path divergence profile is produced. The mean path divergence for selected and empirical landscapes is similar (Figure 6). Closer to the local maximum lower, MPD is observed. MPD increases with distance from the local maximum, levelling out at a local maximum distance of around 3.
The phenome selected landscape is more similar in path divergence to the empirical data than the previously described genome selected landscape (Figure 6). This similarity to empirical data in terms of mean path divergence shows that the phenome selected landscape is more similar in smoothness to the empirical landscape and therefore a more realistic representation of aptamer-protein binding.
4. Conclusions
Herein we have described a method for generating genotype-phenotype interaction maps with lower aggregated pleiotropy vectors which yield smooth fitness landscapes. These fitness landscapes have a steeple type “Massif Central” structure which appears to be more biologically accurate when representing the binding of an aptamer to an analyte.
Furthermore we have removed the one phenotype per gene paradigm that is the norm for model literature and used the varying number of phenotypic contributors to tune our new fitness landscapes. Comparisons between phenome selected landscapes and the similarly constructed genome selected landscapes have been made. The most striking difference is the complementary correlation line of best fit for hamming distance to fittest sequence (Figure 4). Perhaps future work would combine both approaches, constructionally selecting both genotypic and phenotypic contributors, to yield a landscape with a linear correlation line of best fit for hamming distance to fittest sequence. This difference seen in the correlation line of best fit for hamming distance to fittest sequence between genome and phenome relates to both mutational additivity and the global landscape structure. In situations where a higher number of divergent fitness optima should occur, such as genome evolution, genome selected landscapes seem more representative. In situations where a lower number of convergent fitness optima should occur, such as the binding of an aptamer to a specific analyte, phenome selected landscapes seem more representative. Our phenome selected landscape seems to model aptamer-target binding better than any other tested model. A future application of the phenome selected fitness landscape model would be its use in simulation of aptamer selection. As the model seems to be more accurate, the simulation would be a more accurate representation of real life aptamer selection so any selection condition optimisations should be more applicable and transferable to real life aptamer selections. Although phenome selected landscapes seem to represent this biological system more accurately than their standard landscape counterparts, caution should be used when applying them to other biological systems as they may not be well described.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was funded by Hong Kong RGC GRF Grants 17127515 and 17119814 awarded to Julian Alexander Tanner, by HKU Outstanding Young Researcher Award awarded to Julian Alexander Tanner, and by Seed Fund Grant 2014090160039 awarded to Julian Alexander Tanner. The authors thank Dr. William Rowe for providing the allophycocyanin-DNA microarray binding data.
Supplementary Materials
Constructional selection of interaction maps. A) Selective genome growth of an interaction map. Genes (columns) with many phenotypic contributors are added sequentially. B) Selective phenome growth of an interaction map. Phenotypic contributors (rows) with any gene interactions are added sequentially.