Testing biological hypotheses with embodied robots: adaptations, accidents, and by-products in the evolution of vertebrates

Introduction

When properly designed, physically embodied robots may be used to test hypotheses about the workings of animals. This is the goal of biorobotics, as demonstrated and codified by Webb (1995, 2001). For example, biomechanists interested in fish propulsion have used self-propelled, embodied robotic models to test hypotheses about the interactions of mechanical stiffness and control (McHenry et al., 1995; Lauder et al., 2007; Tangorra et al., 2010; Long et al., 2011b; Esposito et al., 2012; Shelton et al., 2014). Biorobotics, with its explicit biological goals, is a natural complement to evolutionary robotics (Harvey et al., 1997; Lipson and Pollack, 2000; Nolfi and Floreano, 2000, 2002; Pollack et al., 2000; Pfeifer and Bongard, 2006; Vargas et al., 2014). Evolving embodied robots to test biological hypotheses offers new modeling tools to biologists and new insights about evolutionary mechanisms to engineers, computer scientists, and cognitive scientists (Floreano and Keller, 2010; Long, 2012; Bongard, 2013; Eiben, 2014a,b).

Evolutionary robotics can help address difficult questions in evolution. For events that occurred in extinct species, a major challenge is the nature of the physical evidence. Fossils, if they are even available, provide only limited information about morphology and environment; they are silent regarding physiology and behavior. Phylogenetic reconstruction can readily infer patterns of trait evolution but says nothing about the real-time dynamics of the specific microevolutionary mechanisms creating those changes. In general, historical analyses of adaptation are hampered by the change in a trait’s function over time and the lack of information about the individual and population-level genetics, behavior, ecology, and selection environment (Brandon, 1990). If one is interested in the relationship between morphology, behavior, performance, and evolutionary mechanisms, then the only recourse is to judiciously design and test models.

Digital simulation, where models of organisms or mobile robots operate and evolve autonomously in a virtual environment, has proven to be a powerful approach in evolutionary robotics. Modeling a quadruped and hexapod, Bongard (2014) found that evolution toward a targeted behavior accelerated when the neural network controllers occupied a succession of different body plans that changed over developmental and evolutionary timeframes. Lipson and Pollack (2000) explored vast regions of morphological space with digital simulation in multiple replicates of mobile organisms, creating physically embodied versions of a few evolved designs to validate their physics engine. Lenski et al. (2003) pioneered the use of digital organisms, populations of computer programs that evolve and can be used to directly test evolutionary theory and the origins of complexity over thousands of generations. Digital simulations of the coupled, closed-loop mechanics of swimming vertebrates tested biological hypotheses of motor control (Ijspeert and Cabelguen, 2006; Ijspeert et al., 2007; Tytell et al., 2010). Most importantly, digital simulation makes possible models with large population size, large numbers of generations, and multiple replicates and starting conditions (Long et al., 2010).

In the face of these obvious benefits of digital simulation why bother evolving physically embodied robots? First, embodied robots offer physical interactions that – by virtue of being instantiated rather than simulated – are more accurate (Pfeifer et al., 2007). For example, digital simulation of the physical dynamics of swimming is particularly challenging given the hydroelastic interactions of a bending body and non-linear fluid forces (Tytell et al., 2010). Second, embodied robots cannot violate the laws of physics. Digital simulations may – for the sake of simplicity, invention, abstraction, ignorance, or elegance – assume the physically impossible, like frictionless joints, infinitely stiff bones, or a smooth, featureless environment. The failure of embodied robots to behave or function as predicted is immediately informative about the workings of the system. Third, the physical is logically and empirically prior to the digital if one’s goal is to model an actual physical system. Usual engineering practice is to validate a digital simulation with a working physical model; differences between the two point to problems with assumptions or mechanisms of the digital model. Fourth, the evolution of physical things, notably robots, not only allows for the testing of scientific hypotheses, but also for the automatic design of purpose-built and intelligent machines, overcoming the limitations of systems designed directly by humans (Eiben, 2014b).

While they offer important benefits over digital simulations, physically embodied robots have costs that are magnified in the context of evolutionary robotics. When we fail to create embodied robots that autonomously reproduce and develop, we must engage in the labor-intensive process of constructing by hand each individual. This problem can be overcome by using a hybrid approach (Eiben, 2014b), evolving and developing just the software controller (Nolfi and Floreano, 2000), or by conducting a portion of the experiments in simulation (Elfwing and Doya, 2014). But if one seeks to evolve morphology in complex physical circumstances, where no valid physics engine exists, then the new bodies must be fabricated each generation. Thus manufacturing time currently imposes a high labor cost on embodied evolution experiments involving morphology (Long et al., 2006; Doorly et al., 2009; Long, 2012). The experimental testing of physically embodied robots is also time-consuming: robots must be built, maintained, and fixed; a strict protocol must be developed and followed for all procedures; and data from each individual in each experiment must be checked for quality, concatenated with other data, and used appropriately in algorithms for calculating fitnesses and the genotypes of the next generation.

In this study, we use physically embodied robots to model trait evolution via three co-operating mechanisms: (1) direct selection, (2) random effects, including mutation, mating, and genetic drift, and (3) indirect selection via genetic or functional correlation with other traits. Changes in a trait that result from the direct selection on that trait are called adaptations. Changes in a trait that result from neither direct nor indirect selection are accidents, and they come about by random effects. Mutation is the most commonly modeled random effect, as it is required to generate the genetic variation upon which selection operates. Another random effect, often avoided by modelers because of its tendency to swamp out selection effects, is genetic drift. Drift is essentially a sampling error that comes about because of the mathematics of small numbers. A small population size allows for genetic drift to operate alongside mutation. Finally, changes in a trait that result from the indirect effects of selection on other traits correlated with the focal trait are often called by-products. The presence of other traits allows for selection to alter its target from the focal trait and to work indirectly, through correlated effects, or to avoid the focal trait altogether (Connor, 2012).

Our goal is to test the following hypothesis about trait evolution: selection for enhanced feeding and also fleeing from predators in early vertebrates was sufficient to increase the number of vertebrae. We test this hypothesis using an evolutionary robotic approach: we employ an embodied biorobotic model of an early vertebrate fish species acting as a prey in a predator-prey ecology. The robots, based on the Tadro-class of robots (Long et al., 2006; Doorly et al., 2009; Long et al., 2011a,b; Long, 2012), are behaviorally autonomous, seeking light as a proxy for feeding, while avoiding a pursuit predator.

Fossil records show that between 550 and 400 million years ago, early fish-like vertebrates evolved at least three morphological traits that have been retained in modern fish species: vertebrae, caudal fins, and lateral lines (Janvier, 2008). Each of these three traits appears in the first 100 million years Paleozoic Era, with proto-vertebrae first observed in Haikouichthys during the early Cambrian (Shu et al., 2003), a caudal-fin-like structure making its appearance in Haikouella during the Cambrian (Chen et al., 1999; Holland and Chen, 2001), and lateral lines evident in jawless fishes from the Ordovician (Janvier, 2008). Because of this shared early history, we selected caudal fins and lateral lines as the companion traits to evolve along with vertebrae.

Vertebrae, the bones that run in series from head to tail to form the vertebral column, show a particularly interesting pattern of evolution. Centra, the structures that form the compression-resisting bodies of vertebrae (Porter and Long, 2010), have evolved and been lost multiple times [for review see Koob and Long (2000)]. Might each case of evolutionary origin have been driven by the same selection pressures? Might vertebrae evolve under multiple selection pressures? Might their repeated origin and loss suggest that they also evolve as by-products of selection on other traits?

To understand the mechanisms that might have evolved vertebrae in vertebrates, we tested three related hypotheses:

Hypothesis 1

Vertebrae evolve as a direct target of selection. We predict that: (a) the number of vertebrae will increase and (b) the selection gradient for vertebrae will be positive under constant selection. If this prediction is upheld, then the evolution of more vertebrae is an adaptation for enhanced feeding and fleeing.

Hypothesis 2

Vertebrae evolve via random processes. We predict: (a) no particular pattern for the evolution of vertebrae, (b) the number of vertebrae will not be correlated with another trait under selection, and (c) the selection gradient for number of vertebrae will be zero or will be small and vary in sign from generation to generation. If these predictions are upheld, then the evolution of more vertebrae is an accident (random) with respect to selection and the correlated effects of other traits.

Hypothesis 3

Vertebrae evolve as an indirect target of selection through correlation among traits. We predict that: (a) the number of vertebrae will increase in positive or negative correlation with a change in another trait, (b) the selection gradient for vertebrae will be zero or negative, and (c) the selection gradient will be non-zero for another trait. If this prediction is upheld, then the evolution of more vertebrae is a by-product of selection on some other trait.

The evolutionary mechanisms described in these three hypotheses may alternate in relative magnitude over time. For example, direct selection may predominate in one generation and random processes may do so in the next. This is possible because even though the fitness function – and hence the selection pressure – remains constant over generational time, the individuals that make up the population do not. By allowing for three different mechanisms to operate in parallel, we are increasing the probability that selection will be refuted as the primary evolutionary driver of this system. Under this condition, failure to refute selection as a primary driver in this small population would constitute stronger evidence for selection than in a large population where the magnitude and importance of drift is reduced. We remain neutral as to whether or not vertebrates evolved primarily in small or large populations; however, it is worth noting that small isolates of larger populations often create conditions amenable to rapid evolution (Ridley, 1996). This study models evolutionary dynamics in small populations, and hence any claims as to the generality of the results are limited to that context.

Materials and Methods

Autonomous Robots

We developed two surface-swimming, physically embodied, autonomous robots to simulate a vertebrate prey (“Preyro”), predator (“Tadiator”), and their behavioral interactions. Only the population of Preyros evolved; the single Tadiator was held constant in morphology and coding, providing a consistent force of selection over the generations. Preyro and Tadiator were first introduced in evolutionary experiments by Doorly et al. (2009) and their design was elaborated in Long (2012). We briefly summarize and update their design here.

Preyro was modeled after a jawless fish from the Paleozoic Era, Drepanaspis gemuendenensis (Figure 1). We examined a fossil (Figure 1A) and a 3D physical model by Anton Fürst (Figure 1B) of Drepanaspis at the Natural History Museum, Vienna, Austria (Figure 1A). Drepanaspis is reconstructed as having a rigid, pan-shaped body propelled by a flexible tail with a backbone that lacks bony vertebral centra; it lacks paired fins but possesses eyes and a lateral line (Tarlo, 1964). Preyro is similar in size and shape to Drepanaspis, with a rigid, pan-shaped body propelled by a flexible tail (Figure 1C).

The body of Preyro was a circular 7.7 cm diameter plastic container (Tupperware, Orlando, FL, USA). To facilitate straight swimming and to damp lateral rocking, a Plexiglas™ keel was added to the bottom of the hull; the keel was trapezoidal (9.3 cm at the base, 5 cm at the top, and 5 cm deep). For tracking via video, Preyro had two LEDs of different colors mounted on the hull at the bow and the stern. Lead weights were used as ballast and for trimming.

For detection of the light source, Preyro had two photoresistors, mounted above the waterline ±45° from the midline, serving as eyespots (Figure 2). Two IR emitter–receivers mounted on the sides were used for predator detection and served as proxies for the lateral line of Drepanaspis. Just as the lateral line informs fish of moving objects in the near vicinity, the IR emitter–receivers informed Preyro of objects within a pre-defined distance determined by ζ, the predator detection threshold, in centimeters.

Preyro also carried a three-axis accelerometer (Wireless Dynamic Sensor System, Vernier Software & Technology, Beaverton, OR, USA), a digital microcontroller (MIT HandyBoard, Newton Labs, Renton, WA, USA), an omni-directional IR transmitter, a servo motor, driveshaft, and a submerged biomimetic tail (see below for details).

The software controller for Preyro was written in Interactive C (v. 4; Newton Labs) as a simple subsumption hierarchy consisting of two behaviors: feeding and fleeing. When feeding, Preyro swam forward with a constant tailbeat amplitude and frequency. To turn toward the light, the center of the tailbeat was shifted in proportion to the light differential detected by the two photoresistors. This resulted in the robot approaching the light source, which served as a proxy for food. Feeding behavior continued until interrupted by the fleeing behavior. The fleeing behavior was modeled after a fish’s fast-start escape response and consisted of three steps: (1) high-amplitude tailbeat to rotate Preyro toward the side receiving the stimulus; (2) high-amplitude swing of the tail back in the other direction to rotate Preyro away from Tadiator and propel it forward; and (3) return of the tail to center and resumption of the regular feeding behavior. A flee attempt was initiated when an object was detected by the lateral IR emitter–receivers, which served as a proxy for the lateral line.

Tadiator differed from Preyro in the following ways. The only sensory input for Tadiator came from an array of four IR receivers mounted on its top (Pololu, Las Vegas, NV, USA). This array monitored the IR signal from the omni-directional beacon placed on top of Preyro. The direction of Preyro from Tadiator was determined by the relative strength of the beacon’s signal in the four quadrants of the array. The bearing of Preyro was used as a turning signal, with Tadiator continuously adjusting its heading to keep the bearing of Preyro near zero degrees, dead ahead. Instead of a biomimetic tail, Tadiator swam with a flapping tail of square plastic. In preliminary trials, Tadiator’s performance was adjusted so that it swam at a speed and with a maneuverability similar to that of Preyro, thus ensuring that the competition between the two was initially even.

Biomimetic Tails

Preyro was propelled by a biomimetic tail (Figure 2). This tail, submerged 6 cm below the waterline, was attached to a 8.25 cm long drive shaft rotated by the servo motor. The tail varied from individual to individual by the morphology of its vertebral column. The morphology and biomechanical properties of these biomimetic vertebral columns are based on those of living shark species (Long et al., 2011a,b). Construction details are reviewed briefly here.

The biomimetic vertebral column consisted of a hydrogel notochord with ring vertebrae. Hydrogels were made of a 0.1 g one-to-one concentration of powdered gelatin (porcine skin, Type A; Sigma, St Louis, MO, USA) dissolved in distilled boiling water (Long et al., 2006). The liquid gelatin was poured into cylindrical Delrin molds of 9.3 mm inner diameter and hardened at 4°C for 1 h. After solidifying, the hydrogels were carefully extracted from the molds and inspected for cracks or tears. Intact hydrogels were cross-linked in a solution of 2.5% glutaraldehyde (stock, 25% EM grade; Polysciences, Warrington, PA, USA) in phosphate buffer solution (0.1 mol 1-1 NaH₂PO₄, 0.15 mol 1-1 NaCl, pH 7.0) on a shaker bed at room temperature for 1 h. Immediately following cross-linking, hydrogels were rinsed thoroughly in distilled water and stored in an aqueous solution of 19.98% ethanol for 12–24 h before use in tail assembly. Once fixed and rinsed, the hydrogel rods were trimmed to a uniform length of 84 mm. Each hydrogel was used within 24 h.

Ring vertebrae were constructed of Delrin tubing (9.6 mm inner diameter; 1.5 mm thickness) cut to a length of 5 mm. The number of ring vertebrae, N, depended on an individual’s genotype (see Section “Evolvable Traits and Their Genetics”). Ring vertebrae were evenly spaced on the 84 mm long notochords and affixed with cyanoacrylate adhesive; as N increased, the number of joints increased and their length decreased.

The biomimetic vertebral column was then used as the axial skeleton of the biomimetic tail. One end of the biomimetic vertebral column was first glued to a rectangular piece of Plexiglas that served as a grip for the drive shaft (Figure 2). To its other end, the biomimetic vertebral column was affixed using thermoplastic glue to a trapezoidal Plexiglas caudal fin. Caudal fins varied in the magnitude of their span, b, as indicated by an individual’s genotype (see Section “Evolvable Traits and Their Genetics”). Finally, this construct was enclosed in a vertical septum, a bilayer of Press and Seal™ (Glad, Oakland, CA, USA). This bilayer septum prevented the cantilevered tail from bending vertically. To account for unavoidable variability in production, we constructed three replicate tails of each genotype. Thus each genotype was tested in the selection trials (see below) by three replicates, each of which was, in turn, tested three times.

Evolvable Traits and Their Genetics

In addition to the focal trait, the number of vertebrae, N, two additional traits were selected to be evolvable in Preyro: the span of the caudal fin, b, and the predator detection threshold, ζ. Each trait had a fixed range in which it could evolve (Figure 3). Preyro took on individual identities in a given generation as defined by the population’s six diploid genotypes, G_i. Each G_i was a 3-tuple, with each element an integer specifying N, b, and, ζ. These are quantitative traits, each modeled as a polygenic system with the number of genes determined by the resolution of the phenotypic increment (see Figure 3). To create independent genetic assortment, and thus avoid genetic correlation among traits, we placed the polygenes for one trait on its own separate chromosome pair.

In each generation, the three of the six G_i with the best relative fitnesses, ω – as determined by the selection experiments and the fitness function (see sections below) – were selected to reproduce. The number of gametes contributed by each G_i to the gene pool was proportional to ω rank: six, four, and two gametes, respectively, for the first, second, and third-place robots. The other three individuals of the population contributed 0 gametes to the gene pool.

To create gametes, the diploid G_i was split into two haploid genotypes, H_i, by dividing by two the value of each entry in the individual’s defining 3-tuple (Figure 4). Each H_i was then mutated by an amount drawn from a Poisson distribution, making the probability biased toward mutations of small magnitude. For mating, the 12 mutated haploid gametes, $H_i^{\ast },$ making up the gene pool, were then randomly combined, creating the six new diploid G_i for the offspring in the next generation. Gametes from the same or different parents had the same probability of mating. Gametogenesis and mating were conducted in simulation using a custom genetic algorithm program written in Java (Version SE 6).

Population size was held constant at six individual G_i. The three replicate biomimetic tails created and tested for each G_i can be thought of as clones.

Selection Experiments

During each experiment a single Preyro and a single Tadiator swam for 3 min in a circular tank with a 3.2 m diameter. The tank was matte black and sat in a black, light-proofed room to prevent any light pollution from influencing the swimming direction of any individuals. A single 100-watt light hung from the ceiling 1 m over the center of the tank, providing the only illumination in the room and serving as the concentrated food source for the light-seeking Preyro (Figure 5).

In each generation, j, three clones of each of six individual, i, G_i were instantiated in Preyro; these 18 different instantiations were pitted separately against Tadiator in three, 3-min trials. The order of the trials with respect to G_i and clone was randomized. To prepare the Preyro for each G_i it was outfitted with the biomimetic vertebral column carrying the appropriate trait morphology for b and N as specified by the G_i. In addition, Preyro’s software was altered to use the appropriate ζ as specified by the G_i.

At the start of each trial, Tadiator and Preyro were released from the same locations on opposite sides of the tank. Trials ran for 3 min without interference from the experimenters. Biomimetic tails were inspected after each trial to ensure that they had not broken. If the tail had broken, the trial was discarded, and a new tail was made, and the trial was repeated at the end of the random trial queue. Each trial was videotaped by an overhead camera (JVC digital video; 30 Hz temporal resolution; 1.2 cm spatial resolution).

Performance Measurements

Performance of each trial of a G_i clone was measured from video tape and the on-board accelerometer on Preyro. Five specific types of performance were measured: (1) average speed, U, (2) average distance to the light source, R, (3) average distance to the predator, D, (4) peak acceleration, a, during fast-start trials (see below), and (5) the number of successful escape responses, θ. These types of performance were chosen to characterize feeding (U, R) and fleeing from a predator (D, a, θ).

At a time resolution of 30 images per second, we manually digitized the locations of the predator and prey robots throughout the trials by tracking the two LEDs on each robot using LoggerPro (v. 3.6.0, Vernier Software & Technology, Beaverton, OR, USA). The absolute position of each robot was determined by taking the average between the positions of the front and back LEDs.

The U was calculated by finite differences of the position of Preyro from frame to frame; a single average U was used to represent this performance in each trial. The R was the average distance of Preyro from the light source over the 3-min trial. The D was the average distance between Preyro and Tadiator over the 3-min trial. The θ was the difference between the total number of escape responses and those not initiated by the predator; this corrected for escapes initiated by interaction with the wall of the tank.

The a was the only performance variable not determined during the selection trials. Because of the noise introduced to the accelerometer by collisions with walls and with Tadiator, we could not determine which parts of the accelerometer record during the experiments corresponded to true predator-mediated escapes. Instead, we measured a immediately after each selection trial by letting Preyro swim without Tadiator; when Preyro was swimming steadily in the center of the tank, we initiated an escape by placing an object near one of the IR sensors. We did this three times for the left IR sensor and three times for the right IR sensor, in randomized order. From the three-axis accelerometer trace for each escape, a was calculated as the resultant of the smoothed peak acceleration of the two vectors in the horizontal plane.

Fitness Function

The fitness, ω, of an individual, i, in generation, j, is defined as the chance of survival of the genotype, G_i (Ridley, 1996). The individual relative fitness, ω_ij, is relative to that of other individuals in a given generation. To calculate ω_ij, we used the average performance values from the three trials of the three clones (n = 9 trials for each G_ij): U_ij, R_ij, D_ij, a_ij, θ_ij. To standardize the different scales of the performance measures, the ω_ij was calculated as a sum of z-scores:

where s denotes the standard deviation of the performance in the population in that generation. Note that the z-score for R is negative; this sign rewards smaller distances from the light source. All other performance measures reward larger magnitudes. The process of selection and mating was repeated for 10 generations.

We created this fitness function for two reasons. First, we sought a fitness function that would reward enhanced feeding and predator avoidance while allowing for trade-offs between the two. From previous work evolving tail stiffness in light-seeking embodied and digital Tadros (Long et al., 2006, 2010), we determined that U and R, among a host of variables, were the performance variables that best predicted the ability to get to and stay near a light source. The other variables used in this current fitness function, D, a, and θ, were chosen to characterize predator avoidance. In fish, D has clear implications, since one cannot be ingested by a predator without being in immediate proximity; a has been shown to be a central performance variable in the fitness of fishes [for review see Ghalambor et al. (2003)]. Finally, θ measures the ability to initiate an escape response; we have shown in guppies that the likelihood of capture by a predator increases if the prey fails to fully initiate a fast-start escape response (Jones et al., 2008).

Second, by using a compound fitness function composed of five performance variables, we are able to study the relationship between morphology, performance, and fitness, and this network is extremely important in evolutionary theory (Kingsolver and Huey, 2003). By defining the causal relationship between performance and fitness, we create a simpler relationship between morphology and fitness that allows us to more easily understand, as a first approximation, the complicated dynamics of the fish predator-prey evolutionary system (Ghalambor et al., 2003). Alternatively, we could have measured those five performance variables (and others) and used a single, independent measure of fitness, like the actual amount of light gathered.

Measuring Selection, Chance, and Correlation among Traits

The selection differential, S, is the difference between the mean value for a trait in the population as a whole and the mean value for that trait in the three G_ij selected to breed. The S estimates the total strength of selection acting on the trait, including the selection acting directly on that trait and indirectly through selection on correlated traits (Lande and Arnold, 1983). The linear selection gradient, β, measures the strength of directional selection acting directly on a character; it is the coefficient for a given trait from a multiple regression of fitness, ω_ij, onto all of the traits. The β is usually given as a standardized coefficient to permit comparison among traits and different studies. We report both raw and standardized values.

To measure the magnitude of random effects on the change in the mean value of a trait, we calculated the strength of chance, C, as the difference between the predicted change in the value of the mean of the trait if S alone were operating and the actual realized change in mean. S and C are component vectors of the resultant evolutionary change.

Correlation among the traits was measured using the Pearson correlation coefficient, r, and the partial correlation coefficient, ρ. The ρ removes the effect of the correlation with other variables as it measures the correlation of a particular pair of traits. To find the relative magnitude of the indirect correlated selection, χ, for a focal trait, x, in a given generation, we calculated:

where i indicates the other traits. All of these measures are reported from the generation in which they were calculated; note that they pertain to the changes occurring in the mean value from that generation to the next. Although we conduct standard statistical tests, please note that we are working with the entire population, and not a sample. With the whole population in hand, S, β, C, χ, and differences in means from generation to generation are therefore the actual values and we need not rely on significance tests to estimate the presence of a difference. However, in order to provide a conservative interpretation of our results we provide statistical tests with α = 0.05 for analysis of all genotypes (n = 66, 6 for each of 11 generations), α = 0.10 for analysis of correlation between traits within a generation (n = 6), and α = 0.05 for analysis of correlations among S and β across generations (n = 10). JMP (version 10.0.0) was used for all statistical analyses.

Results

Evolution of Traits by Chance Alone

To ascertain the magnitude and possible patterns of evolution by chance alone, we ran our genetic algorithm without inputs from selection experiments. With the same mean and variation for the starting condition as the trials for the embodied Tadros, we ran three different simulation trials using the genetic algorithm without selection (Figure 6). Note that in some cases purely random evolution will produce patterns that are directional over multiple generations. Thus selection is not the only way to achieve directional evolution over short time frames. Moreover, the genetic mechanisms that we include in our model – mutation, mating, and genetic drift – can combine to create intergenerational changes in the traits of the same magnitude that we see in the experiments with the robots.

Evolution of Traits by Selection, Chance, and Correlation

Over 10 generations of selection, the population’s mean number of vertebrae, N, increased significantly (Figure 7A; n = 66, r² = 0.125, F = 9.159, p = 0.004). Likewise, the caudal fin span, b, increased significantly as determined by linear regression (Figure 8A; n = 66, r² = 0.109, F = 7.795, p = 0.007). However, predator detection threshold, ζ, did not increase significantly in a linear regression (Figure 9A; n = 66, r² = 0.005, F = 0.325, p = 0.5706). Instead, ζ was significantly fit by a model that included a linear and a quadratic term for the independent variable generation (n = 66, F = 5.664, p = 0.0055, r² = 0.152; both linear and quadratic terms were significant, p < 0.05, with the linear term of positive sign and the quadratic of negative). When a quadratic term was included for N and b, in neither model was the quadratic term significant; hence, we used the original linear models. In summary, for all three traits, univariate regression analysis detects significant positive and linear terms for each, which is evidence of directional evolution on average over the 11 generations analyzed.

Please note that the statistical analysis represents a conservative and simplified approach to the complexities of the evolutionary system. Since we have all individuals in the population, one could argue that statistics that assume a sample as representative of the larger, unsampled population are unnecessary. From that perspective, if we look at the data from generation to generation (Figure 7), we see that the population’s mean value of N increases in generations 2 and 3, oscillates, and then achieves a steady state in generation 7. Likewise, the population’s mean value of b initially decreases before increasing from generation 4–11 (Figure 8). The two approaches are compatible: the statistical analysis detects simple, longer-term patterns while we can discuss changes from generation to generation without statistics.

Selection, Chance, and Correlation

The selection differential for N, S_N, is positive in generations 1 and 2, and that corresponds to intergenerational increases in the mean N from 1 to 2 and 2 to 3, respectively (Figure 7). In generation 1, direct selection alone is operating on N. While S_N is positive, the strength of chance, C_N, is zero. The standardized selection gradient, β_N, of 1.39 indicates the relative strength of the direct selection on N. The lack of any correlated indirect selection, χ_N, rules out the third class of possible mechanism.

By generation 2, the picture is more complicated. While S_N is positive, and the mean of N increases from generation 2 to 3, β_N is nearly zero, ruling out a role for direct selection on N. Instead, indirect selection is a driver, as indicated by a total χ_N of 1.09, with all of that positive magnitude contributed by the correlation of N and ζ. Chance, too, plays a role, as shown by C_N, which has the same positive value as S_N.

In generation 3, we see yet another type of interaction between selection, chance, and correlation. While the S_N value of 0 would seem to indicate that no selection, direct or indirect, is acting, the β_N of 1.03 shows positive direct selection. This direct selection is counterbalanced by a negative C_N and a very strong negative χ_N, once again due mostly to the correlation of N and ζ.

Generation 4 is notable because of the small role of β_N and C_N; a small χ_N of 0.18 drives up the mean N of the population and accounts for the positive S_N. Generation 5 is the first in which we see a strong negative β_N, which wins the day in terms of decreasing the mean of N, in the face of positive C_N and positive χ_N. Generation 6 shows an immediate reversal of the sign of selection, with a strong positive β_N and an increase in the mean of N, in the face of a very strong negative χ_N, with little effect of C_N. Generation 7 shows another reversal in the sign of β_N, now strongly negative; however, in this case the mean N of the population is stable inter-generationally; hence the large-magnitude of C_N, coupled with a positive χ_N, balances β_N. The stable value of the mean N for the remaining generations is explained by equal and opposite values of S_N and C_N coupled with small, oscillating values of β_N and small positive values of χ_N.

While neither b nor ζ is the focal trait, each has the potential to be under direct selection. The overall statistical pattern of the mean of b increasing over the 11 generations comes about by variable contributions from selection, chance, and correlation (Figure 8). The largest values of β_b occur in generation 5 and 6 and are negative (Figure 8C). Generation 6 is interesting since the large negative β_b is accompanied by a large negative C_b and the largest χ_b, which is positive. Even though the total selection, S_b, is positive with a value of 1.0, the C_b is larger still, with a negative value of −1.17; the result is that the mean b in generation 7 decreases slightly. Generation 8 is also dominated by chance but in a different way than generation 6. Here, the C_b is large and positive, but S_b is small and negative. The β_b is small and positive, and χ_b is close to zero. Thus C_b is the only strong effect and it accounts for the abrupt increase in mean b from 24.67 to 26.33 in generations 8 and 9, respectively. As with N, the evolution of b is a mix of selection, chance, and correlation that changes from generation to generation.

Compared to N and b, the quadratic pattern of the change in the mean of ζ differs in its balance of underlying mechanisms in important ways (Figure 9). First, the variance in the value of ζ drops dramatically from generation 2, is close to 0 in generation 5, and then increases again. Without variance in a trait, selection cannot operate; that is what we see in generation 5, with near-zero values of S_z. Second, in generation 5 we see little or no impact of selection, chance, or correlation. The value of C_z is 0 and the value of standardized β_ζ is 0.97, the smallest absolute measured for β_ζ with the exception of that for generation 10 (0.09); χ_ζ is also of small magnitude. It is no surprise then that the mean of ζ in generation 6 is virtually identical to that in generation 5. Third, we find the largest negative value of β, −5.79, in any of the traits, in generation 3. Fourth, with the exception of generations 6 and 7, the effect of indirect, correlated selection on ζ is small, as measure by χ_ζ. Thus unlike N or b, the pattern of evolution of ζ is largely explained by an interplay of just direct selection and chance.

The variation in the dominant evolutionary mechanism from generation to generation is explained, in part, by the changing pattern of correlation among the traits (Figure 10). For example, while the partial correlations, ρ_{N- b} and ρ_N-ζ are negative in generation 5, the ρ_N-b switches sign to become strongly positive in generation 6 while ρ_N-ζ becomes more strongly negative. Because of this, and in conjunction with changes in β_β and β_ζ, we can explain the rapid oscillation of indirect correlation selection on N (Figure 7D).

Evolution of Performance

The five different types of performance that we measured for the fitness function underwent different patterns of evolutionary change (Figure 11). Generations 5, 6, and 7, which show the strongest direct selection on N, as measured by β_N (Figure 7C), show a corresponding connection between negative β_N in generations 5 and 7 and the resulting increases in the population mean of swimming speed, U, in generations 6 and 8. The positive β_N in generation 6 results in a decrease in the mean U in generation 7. Overall, mean U decreases initially, oscillates, and then returns from the negative range to where it started in generation 10. The population mean for the distance to the light, R, increases and stays in the positive range, relative to its starting point, until generation 10. The population mean for the distance to the predator, D, drops initially and only by generation 8 is positive relative to where it started. The population mean for peak acceleration, a, increases in generation 3, drops back to 0 in generation 5, before staying in the positive range. The population mean for the number of escapes, θ, decreases into negative range at generation 4.

Although we have focused on the mean performance values, keep in mind that only the genotypes from breeders in each generation (indicated by gray markers) are those selected to pass on traits to the next generation. Thus the position of the breeders relative to the mean value is a proxy for the contribution of that type of performance to fitness. If the performance types uniformly mapped onto fitness we would expect all of the means of breeders to be positive relative to the mean of the population, with one exception: R, by virtue of smaller values being rewarded, should be negative. We do not always see that simple pattern. For example, sometimes we see no difference between the mean of the breeders and the mean of the population; in those instances that performance type is not contributing to the differences in fitness among individual genotypes.

Generation 2 serves as an example. The performance values for R and θ differ in the expected directions from the population means; U, D, and a do not. Thus in generation 2, when N increases without direct selection (Figure 7C), we would postulate that selection on N requires a connection to U, D, and/or a. The positive χ_b, mostly from the N-z correlation (Figure 7D), suggests that the changes in θ are directly linked to the positive direct selection on ζ (Figure 9C) and, through the correlation, to N. By examining connections like these among morphology, performance, and fitness we begin to gain insight into the causal feedback between selection on performance and the pathway through which selection, both directly and indirectly, act on a trait. We will undertake a full analysis of this network in a future study.

Discussion

In a population of autonomous, embodied, fish-like robots, the mean number of vertebrae, N, in the vertebral column of the propulsive tail evolves under selection for enhanced performance in feeding and fleeing. The population’s mean value of N increases from 4.7 to 5.7 over 10 generations (Figure 7). From biomechanical studies, we know that increasing N increases the apparent storage modulus, E′ (MPa), and apparent loss modulus, E″ (MPa), of a biomimetic vertebral column (Long et al., 2011b). From that same study, we know that as N increases, so too does the steady swimming speed and peak acceleration of a swimming robot using the biomimetic vertebral column in its propulsive tail. Thus, at one level our explanation of adaptation seems complete using our robotic system as a model for early vertebrates. We were unable to refute our biological hypothesis: selection for enhanced feeding and fleeing in early vertebrates was sufficient to increase the number of vertebrae, N.

But we must interpret cautiously. First and foremost, please keep in mind that we have not run multiple replicates of this population of embodied robots in this study. We have, however, published a preliminary study using the identical set-up; this population ran for only six generations (Doorly et al., 2009) and yielded nearly identical patterns of evolution for N as seen in the first six generations in this study. In addition, we have used digital simulations to examine the repeatability of evolutionary results on the Tadro-class of robots (Long et al., 2010). We ran 32 replicates of a population of 70 light-seeking, one-eyed digital Tadros for more than 500 generations. Selecting for enhanced light-seeking, without a predator, we found that structural stiffness of the tail, proportional to the number of vertebrae (Long et al., 2011a), repeatedly evolved to a single, stable global optimum.

Second, these results pertain to populations of very small size, where mutation and drift, random effects of chance, may operate at a magnitude equal to that of selection. In vertebrates in the wild, small populations under new predation pressure are noteworthy for their rapid and large-magnitude evolutionary change (Reznick et al., 1997; Grant and Grant, 2002; Losos et al., 2006). This empirical evidence supports the hypothesis that small, reproductively isolated populations of a species may offer important opportunities for speciation (for review, Losos and Glor, 2003). Thus populations of small sizes are biologically relevant.

Third, because of the time it takes to build, run, and analyze a population of embodied robots with evolving morphology (see “Introduction”), we have allowed the population to evolve for only 10 generations. Thus claims of testing evolutionary phenomena must be qualified with the important caveat that we are allowing very little time. However, we note the rapid evolutionary changes in living populations of vertebrates in the wild. With as little as one generation of selection by introduced predators, Anolis lizards show rapid evolution of behavior and morphology (Losos et al., 2006). Small, isolated populations of Geospiza ground finches show dramatic changes in morphology annually (Grant and Grant, 2002). Male Poecilia guppies in small populations exposed to an introduced predator show significant changes in size and life history traits in just 4 years or about 12–20 generations (Reznick et al., 1997). Thus in the wild significant evolutionary events in vertebrates may occur in the smallest time scale possible, between single generations.

Selection, Chance, and Correlation

Our initial interpretation of the results of this evolutionary robotic model is incomplete. First, we need to consider in more detail other aspects of the model, in particular the alternative evolutionary mechanisms to direct selection. Next, we need to consider what those more complicated results tell us about the target. In short, the evolution of vertebrae in robots and, likely, vertebrates, is not simply about direct selection acting to drive the adaptation of a single trait.

The mechanisms that constrain adaptation are as important as those that create it (Kingsolver and Huey, 2003). Correlations among performance elements indicate the degree of functional constraint or facilitation of the evolution of morphology (Lande and Arnold, 1983; Ghalambor et al., 2003; Walker, 2007). Genetic correlations may also constrain or facilitate the evolution of morphology; however, in this study we have removed the possibility of genetic correlations from our model by allowing independent assortment of alleles. To allow for the possibility of functional constraint by correlation, we modeled the evolution of N with two other evolving traits, tail span, b, and predator detection threshold, z.

What appears to be a simple, linear increase in the mean value of N is not. First, the linear regression, while significant, misses the fact the mean N for the last five generations is constant (Figure 7A). Second, direct selection on N, as measured by the selection gradient, β, varies from generation to generation in its magnitude (Figure 7C) and in its relative importance as an evolutionary driver. In the first generation, direct, positive selection is the primary cause of the increase in mean N in generation 2. But in generation 2, direct selection is negligible, and chance and indirect correlated selection drive the mean of N to the value seen in generation 3. In generation 3, all three mechanisms operate. Depending on the generation in which one looks, you could say that the change in N is an adaptation, an accident, or a by-product.

These results show the value of evolving a single trait in conjunction with other traits in a small population. Local optima are less likely to be terminal for the focal trait if indirect selection, acting on other traits, is present. Small populations guarantee the presence of drift and hence chance. One might worry that direct selection would be swamped under these conditions. But we show that is not the case necessarily. For example, strong direct selection on N in generations 5, 6, and 7, switching signs from negative, to positive, to negative, may dominate in the first two cases and be counterbalanced in the third by chance and correlation.

In the face of all three of these evolutionary mechanisms in operation, a trait can only be said to be an adaptation for a particular type of performance without qualifications if the change in its mean value in a population, over generational time, can be shown to be caused solely or primarily by direct selection on that trait. If we had designed this study to use a very large population, then we could have reduced the relative impact of chance via drift on the evolution of N. If we had allowed only one single trait, N, to evolve, then we would have eliminated the opportunity for indirect correlation selection. While these caveats are well-known to evolutionary biologists, they are worth keeping in mind as we use robots to model biological evolution.

Testing the Mechanistic Hypotheses

We tested the predictions that followed from three mechanistic hypotheses. While applied here to the evolution of N, they provide a starting point for understanding the mechanisms driving the evolution of any single trait. We review each in turn.

Hypothesis 1

Vertebrae evolve as a direct target of selection. We predicted that (a) the N would increase and (b) the selection gradient for N would be positive under constant selection. This was true in some generations but not others, so it would seem unwise to make the blanket claim that a greater N are an adaptation for enhanced feeding and fleeing.

One argument in favor of the adaptation hypothesis is that a mean N of 5.7 vertebrae may be a strong local optimum. Our biomechanical studies of the impact of N on swimming performance show an abrupt drop-off in swimming speed of 20% and peak acceleration of 30% when N increases from 6 to 7 (Long et al., 2011b). Looking just at the selection gradients (Figure 7C), they alternate sign from generation 4 through 10, with four reversals. Thus it is tempting to characterize this pattern as a case of stabilizing selection. From that perspective, the large amount of chance present would act to maintain genetic variance, which is usually lost under stabilizing selection. Compared to a population with stabilizing selection but diminished variance, this population would respond more quickly to a change in the environment, showing greater evolvability.

Another way to test the adaptation hypothesis is to repeat the experiment. If the mean of N evolves in the same way, with chance and correlation present, then we can assert that the total selection, acting directly and indirectly as measured by the selection differential, S, (Figure 7B) is the primary driver. In an earlier preliminary experiment (Doorly et al., 2009), we saw a nearly identical pattern, with a positive increase in N to just under 6.0; however, the experiment lasted only for five generations, so we cannot compare later generations. Also, in the previous study we did not calculate S or β, so we do not know the impact of direct selection.

Hypothesis 2

Vertebrae evolve via random processes. We predicted (a) no particular pattern for the evolution of N, (b) N would not be correlated with another trait under selection, and (c) the β_N would be zero or would be small and vary in sign from generation to generation. Over 11 generations these predictions were not upheld. Thus, in spite of the clear presence of strong effects of chance, as measured by C_N, the impact on evolutionary response was intermittent. In addition, N was correlated with the other traits and β_N was present. The evolution of N appears then to be accidental only in a few generations.

We were concerned on the one hand that the random effects of our model were of sufficient magnitude to generate and maintain variance, while on the other hand not swamping selection. The impact of accidental evolution may be increased by increasing the magnitude of mutation and decreasing the size of the population. Our data on the evolutionary patterns created by chance alone (Figure 6) indicate that these random effects were titrated properly, such that they could create directional patterns equal in magnitude to those we saw in our experiments. Despite the opportunity for chance to drive the evolution of N to the same degree that we see with all three mechanisms present (Figure 7), chance did not dominate.

Hypothesis 3

Vertebrae evolve as an indirect target of selection through correlation among traits. We predicted that (a) N would increase in positive or negative correlation with a change in another trait, (b) the selection gradient for N would be zero or negative, and (c) the selection gradient will be non-zero for another trait. As with our chance hypothesis, indirect correlated selection occurs in some generations and not in others (Figure 7D). Thus we cannot say that increased N is always a by-product of selection on b or ζ. Indirect correlated selection on N is present, important in some generations, but not dominant overall.

One reason that indirect selection varies is that the correlation among the traits varies from generation to generation (Figure 10). We did not anticipate this changing pattern, in spite of the complexities of functional constraint (Walker, 2007) and co-variation (Walker, 2010). Because both N and b are morphologies of the tail that impact propulsion, we expected that they would be positively and constantly correlated; a stiffer backbone (higher N) is needed to withstand the greater hydrodynamic loads imposed by a larger caudal fin (higher b). This positive correlation between N and b occurs in only 6 of 11 generations. Meanwhile, we expected ζ, a sensory trait, to be uncorrelated with either N or b. The fact that we do not yet understand the reasons behind these functional correlations shows the importance of linking morphology, and correlations among morphologies, to different types of performance (Figure 11). This work is forthcoming using the analytical framework codified by Walker (2010).

One way to increase the complexity of indirect selection is to allow genetic correlation. By design, we did not. Genetic correlation – in its many forms – is likely to provide more and different patterns of indirect selection than with functional correlation alone. In this study, we have a simple one-to-one mapping of the genotype to the phenotype. Interactions among genes may be modeled as epistasis, with large possible impacts on the evolvability of the population.

Robots as biological models

How do the results from evolving robots relate to ancient vertebrates? One could argue that there is no relation between the two and thus we learn nothing about vertebrates by studying robots. We beg to differ. The initial point of discussion is this: what constitutes a useful model and does ours qualify as such? With biorobotics in mind, Webb (2001) proposed seven dimensions to characterize and evaluate any model simulation. Using slightly modified terminology (Long, 2007), those dimensions are (1) biological relevance in terms of testing a hypothesis, (2) match of the performance of the model compared to the target, (3) accuracy of mechanisms represented by the model, (4) abstractness of the model, (5) level in structural hierarchy of the target that the model represents, (6) specificity as to the number of targets represented by the model, and (7) substrate out of which the model is created. For Webb, useful biorobotic models must, at a minimum, have biological relevance and be of physical substrate.

The basic logical requirement for usefulness is to establish that a model represents its target. The nature of that representation is made explicit by use of Webb’s hyperspace (Webb, 2001). Once representation is established, it follows that tests of the model are tests of the target. Thereafter, fruitful interpretation of the experiments is not about whether the model is “good” or “bad” or whether we learn anything about the biological target, but, instead this: what exactly does the model tell us about the target, what level of confidence do we place in those results, and how might we improve the model?

We assert that the model of an evolving population of Preyros represents an evolving population of early vertebrates in the following ways. First, we tested a hypothesis about the selection pressures evolving vertebrae, thus fulfilling the need for the biological relevance of the model. Second, while we do not have the ancient vertebrates available to match performance of behaviors, we do see feeding, fleeing, and chase behaviors in our robots, and those behaviors are found in living fishes. Third, we use mechanisms seen in living fish: the biomechanics of the biomimetic vertebral column, sensory-motor responses, the neurally inspired subsumption architecture for the controller, and the evolutionary mechanisms themselves. Fourth, our models are simplifications of the target, concrete in the sense of modeling specific traits and processes. Fifth, the lowest level of our model is at tissues and materials of the backbone; the highest at the level of populations. Sixth, we attempt to be specific in our modeling by focusing on a single target, the early vertebrate Drepanaspis. Seventh, our robotic models are physically embodied.

What exactly does the evolution of a population of Preyros tell us about the evolution of vertebrae in early vertebrates? First, models have explanatory limits (Rosenblueth and Weiner, 1945). An evolutionary robotic model of past events, no matter how accurate the modeled agents, can never reconstruct the actual events with complete certainty. Fossils simply do not provide enough information about physiology, behavior, and ecology. Rather than claiming that we have discovered the exact pattern of causes driving the evolution of a biological trait, the best that we can do is work to refute specific hypotheses and, in so doing, eliminate the less-likely explanations. This leaves us with what Brandon (1990) calls “how-possible” explanations. These, in turn, are subject to further testing and possible falsification.

In our first work on the evolution of the embodied Tadro-class robots, we re-analyzed our original work (Long et al., 2006) and refuted the hypothesis that the stiffness of the backbone, which is proportional to N (Long et al., 2011b), is an adaptation for enhanced feeding behavior in robots or in early vertebrates (Long, 2012). Feeding alone as selection force evolving N is thus placed in the “unlikely” category of explanations. To the fundamental behavior of feeding we added predator avoidance and were unable to refute the hypothesis that N evolved in robots and vertebrates as an adaptation for enhanced feeding and fleeing performance (Doorly et al., 2009). This current follow-up study, with twice the number of generations, also was unable to refute the adaptation hypothesis for feeding and fleeing. However, as we have shown, the picture, even in this simple model system, is more complicated than answering, “Adaptation: yes or no?”. Direct selection on N is not the sole evolutionary driver. Chance and indirect selection, through character correlation, both act to augment, constrain, or supersede direct selection.

Given these important qualifications, we tentatively place the main research hypothesis of this paper into the “how-possibly” category for the evolution of early vertebrates: selection for enhanced feeding and fleeing in early vertebrates was sufficient to increase the number of vertebrae, N.

Looking ahead, we see at least eight different ways to improve our model. First, Tadros swim on the surface of the water, in two dimensions. To improve behavioral match, these models must be submerged and capable of movement and maneuvers in three-dimensions. Second, we have not explored how the gene regulatory networks might constrain or augment evolvability of our population. Allowing for sign epistasis among traits, for example, would allow for the possibility of dramatic genetic correlation among traits. Third, given the importance of indirect, correlated selection (Walker, 2007, 2010), we could add more evolvable traits to the Preyro. The tricky research question is: which ones? Changes two and three would interact to increase the chance for the evolution of by-products. Fourth, the population size could be increased, thus reducing the effects of genetic drift and testing the hypothesis that direct selection on N could become a more dominant evolutionary driver by changing this factor alone. Fifth, we could simply allow more generations of evolution, thus allowing us to test, for example, if the population has really reached an equilibrium state of N or if it can migrate from what looks like a local optimum? Sixth, we could test new types of selection pressures. Letting the predator co-evolve (Nolfi, 2012; Elfwing and Doya, 2014), for example, would add in realistic complexity to selection landscape, as would variation on other environmental parameters (Kingsolver and Gomulkiewicz, 2003). Seventh, and related to the sixth point, we could use different fitness functions. The most biological would be the number of offspring reproduced. Reproduction in evolving robots can be linked directly to the amount of energy gathered (Elfwing and Doya, 2014). Eighth, we could allow the predators and the prey to learn in response to their interactions, as this ability has been shown to have important effects on the evolution of behavior in robotic systems (Elfwing and Doya, 2014).

Conclusion

Autonomous, physically embodied robots permit physical interactions, behavior, ecological interactions, performance, and fitness to be outputs of rather than inputs to the evolutionary model. With physical models, physics need not be modeled because it is inherent. Autonomy is created (1) by the continuous feedback between sensory input and movement output that creates behavior, and (2) by removing a human operator from that control loop. In summary, evolutionary robotic models that allow adaptations, accidents, and by-products offer enhanced realism and complexity that permits more robust tests of biological hypotheses compared to models that focus on single traits and the single mechanism of adaptation by selection.

Eiben (2014a) has identified three grand challenges for evolutionary robotics: (1) to demonstrate that evolution can create a surprisingly novel robot; (2) to create physical robots that can reproduce; and (3) to evolve physical robots in the real world in an open-ended process. While this study succeeds at none of those challenges, we have made some headway toward the “evolution of things” (Eiben, 2014b). Our robotic things, the Tadro-class Preyros, are physically embodied and behaviorally autonomous. While they behave in vivo, as it were, Preyros start the reproductive cycle in silico, with mutation and mating algorithms creating new genotypes for offspring. With those genotypes, the equivalent of development is done per manum, by the hand of a human builder. Thus the evolution of a population of Preyros, while succeeding in evolving physical robots, fails to do so in a hands-off, open-ended process.

What is a grandly challenged evolutionary roboticist to do? We can frame the challenges in terms of the three hierarchically nested timescales described by Pfeifer and Bongard (2006): (1) the behavioral here-and-now, (2) the ontogenetic (reproductive and developmental), and (3) the phylogenetic (evolutionary). To solve the grand challenges of evolutionary robotics, our physical robot individuals and populations need four kinds of hierarchically nested autonomy operating simultaneously: (1) behavioral autonomy, (2) reproductive autonomy, (3) developmental autonomy, and (4) evolutionary autonomy.

Author Contributions

All authors (1) contributed to the design of the work; (2) made revisions regarding methodology, analysis, and/or interpretation; (3) approve of the final version to be published; and (4) agree to be accountable for all aspects of the work. In addition to those contributions, Sonia F. Roberts wrote the first draft; Sonia F. Roberts and Jonathan Hirokawa supervised data reduction and initial analysis and statistical testing; Marianne E. Porter supervised the manufacturing of the biomimetic vertebral columns and undertaking of the selection experiments by Sonia F. Roberts, Jonathan Hirokawa, Hannah G. Rosenblum, Hassan Sakhtah, and Andres A. Gutierrez; John Howard Long Jr. conducted the final data analysis and wrote the final draft.

Conflict of Interest Statement

John Long was paid to write a book, Darwin’s Devices (2012, Basic Books), that discusses the robots and the experiments in this manuscript. Other than that potential conflict of interest, the authors declare that this research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by the URSI Program at Vassar College and the National Science Foundation through grants INSPIRE 1344227 and DBI 0442269. Lively discussions with Chun Wai Liew, Rob Root, Ken Livingston, and Nick Livingston were instrumental in developing this project. Carl Bertsche, John Vanderlee, and Larry Doe provided valuable assistance constructing the robots. Nicole Doorly, Gigi Engel, Elise Stickles, Gianna McArthur, Carina Frias, Kira Irving, Nicole Krenitsky, and Keon Combie were part of the Fish Fellows program and supported this work.

Abbreviations

α, significance threshold; β, directional selection gradient; χ, indirect selection gradient; ρ, partial correlation coefficient; θ, number of escapes; ζ, predator detection threshold; ω, individual relative fitness; a, acceleration, peak; b, tail span; C, strength of chance; D, distance to predator, average; E′, storage modulus; E″, loss modulus; G, genotype, diploid; H, genotype, haploid; H*, mutated haploid genotype; N, number of vertebrae; r, Pearson correlation coefficient; R, distance to light source, average; S, selection differential; U, swimming speed, average.

References