2 Genotypes and phenotypes
2.1 Overview: From fitness numbers to folded RNA
In practical 1, you studied evolution where fitness was just a simple number (0-1). But real organisms don’t have abstract fitness values; their fitness comes from their physical and chemical properties. This practical builds a bridge from abstract models to real biology by exploring how genotypes (DNA sequences) translate into phenotypes (proteins and RNA structures) that then determine fitness.
We’ll progress through increasing biological realism:
- Fitness as a number, where we left of last time.
- DNA sequences, individuals are DNA strings evolving toward a target sequence.
- Proteins, individuals are DNA strings that translate into target protein sequences. Now the genotype is no longer what is directly selected on!
- RNA with structure — Use the Vienna package to fold RNA and see how these structures evolve. Fitness depends on whether the structure matches a target-fold.
This progression shows why the genotype-phenotype mapping is so important in evolution. A small DNA change might have no effect on the protein (silent mutation), a big effect, or even a beneficial effect. You’ll find that evolution can “explore” this large genotype space in quite surprising ways.
Let’s start:
2.2 Part 1: Fitness as a number
Copy the code below. Notice that calculate_diversity() is empty — it just returns 0. Your job is to fill it in! The diversity line is already plotted, waiting for you to implement the function so it will track the diversity within the population.
import random
import math
import matplotlib.pyplot as plt
random.seed(42)
# --- PARAMETERS ---
population_size = 500
generations = 200000
mutation_rate = 0.001
no_reproduction_chance = 1
sample_interval = 100
# --- SETUP ---
def mutate(fitness):
if random.random() < mutation_rate:
return min(1.0, max(0.0, fitness + random.uniform(-0.1, 0.1)))
return fitness
def calculate_diversity(population):
"""Calculate diversity of the population - NOT YET IMPLEMENTED!"""
return 0
# Initialize
population = [0.1] * population_size
avg_fitness = []
diversity_over_time = []
# Plotting
fig, ax1 = plt.subplots(figsize=(12, 6))
# Left axis: fitness
ax1.set_xlabel("Generation")
ax1.set_ylabel("Average Fitness", color='tab:blue')
ax1.set_ylim(0, 1)
line1, = ax1.plot([], [], color='tab:blue', linewidth=2, label='Fitness')
ax1.tick_params(axis='y', labelcolor='tab:blue')
# Right axis: diversity
ax2 = ax1.twinx()
ax2.set_ylabel("Diversity", color='tab:green')
line2, = ax2.plot([], [], color='tab:green', linestyle=':', linewidth=2, label='Diversity')
ax2.tick_params(axis='y', labelcolor='tab:green')
fig.suptitle("Evolution of Fitness and Diversity")
fig.tight_layout()
fig.legend(loc='upper right')
plt.grid(True)
plt.ion()
# --- EVOLUTION LOOP ---
for gen in range(generations):
best = max(population)
# Moran process: roulette wheel with no-reproduction option
competitors_with_dummy = list(range(population_size)) + [None]
probs = population + [no_reproduction_chance]
winner = random.choices(competitors_with_dummy, weights=probs, k=1)[0]
if winner is not None:
dead_idx = random.randint(0, population_size - 1)
population[dead_idx] = mutate(population[winner])
# Sample and plot
if gen % sample_interval == 0:
avg_fitness.append(sum(population) / population_size)
diversity_over_time.append(calculate_diversity(population))
x_vals = [i * sample_interval for i in range(len(avg_fitness))]
line1.set_data(x_vals, avg_fitness)
line2.set_data(x_vals, diversity_over_time)
ax1.relim()
ax1.autoscale_view()
ax2.relim()
ax2.autoscale_view()
fig.suptitle(f"Generation {gen}: Best Fitness = {best:.3f}")
plt.pause(0.001)
plt.show()2.3 Diversity patterns
Diversity is an important component of evolution; it’s the fuel on which selection acts. But diversity can be measured at meany levels. Let’s think about this for a bit, by modifying the code above.
Exercise 2.1 (The dynamics of diversity - Biology)
Modify the function calculate_diversity to calculate the diversity of the population by taking the standard deviation of the fitness values.
- Use the standard (low) mutation rate and study the dynamics of diversity. Describe the pattern verbally.
2.4 Evolving a DNA sequence
A big problem with the previous model is there is no true distinction between genotype (that which mutates) and phenotype (that which is selected). Let us try and adapt the model to be more biologically meaningful, by making each individual represented by a DNA sequence. For this, copy the code below1:
import random
import math
import matplotlib.pyplot as plt
from collections import Counter
# set the random number seed
random.seed(0)
plt.ion() # Enable interactive mode
# Parameters
alphabet = "ATCG"
target_sequence = "GATGCGCGCTGGATTAAC" # Example target sequence
dna_length = len(target_sequence)
target_length = len(target_sequence)
# Simulation settings
population_size = 500
generations = 20000
mutation_rate = 0.00001
sample_interval = 10
sample_size = 10
no_reproduction_chance = 0.1
# Core functions
def fitness(dna):
"""Calculate fitness as % of correct bases"""
matches = 0
for i in range(len(dna)):
if dna[i] == target_sequence[i]:
matches += 1
return matches / target_length
def mutate(dna, rate=mutation_rate):
return ''.join(
random.choice([b for b in alphabet if b != base]) if random.random() < rate else base
for base in dna
)
def count_beneficial_mutations(dna, num_tests=100):
"""Test random mutations (faster than testing all possibilities)"""
f0 = fitness(dna)
count = 0
for _ in range(num_tests):
i = random.randint(0, len(dna) - 1)
b = random.choice([x for x in alphabet if x != dna[i]])
mutant = dna[:i] + b + dna[i+1:]
if fitness(mutant) > f0:
count += 1
# Extrapolate: if 10 out of 100 random tests were beneficial,
# estimate ~(10/100) * (len(dna) * 3 bases) = total beneficial mutations
total_possible = len(dna) * 3
return count * (total_possible / num_tests)
def diversity(pop):
counts = {}
for ind in pop:
counts[ind] = counts.get(ind, 0) + 1
total = len(pop)
return -sum((c/total) * math.log(c/total + 1e-9) for c in counts.values()) if total > 0 else 0
# Initialize population
initial_sequence = "GATAGCGAAGTTTAGCCG" # far from target (only first 3 are correct)
population = [initial_sequence for _ in range(population_size)]
avg_fitness = []
avg_beneficial = []
diversity_over_time = []
best_individuals = []
# Initialize interactive plot
fig, ax1 = plt.subplots(figsize=(12, 8))
ax1.set_xlabel("Generation")
ax1.set_ylabel("Average Fitness", color='tab:blue')
ax1.set_ylim(0, 1)
line1, = ax1.plot([], [], color='tab:blue', linewidth=2, label='Fitness')
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax2 = ax1.twinx()
ax2.set_ylabel("Beneficial Mutations / Diversity", color='tab:purple')
line2, = ax2.plot([], [], color='tab:purple', linestyle='--', linewidth=2, label='Beneficial Mutations')
line3, = ax2.plot([], [], color='tab:green', linestyle=':', linewidth=2, label='Diversity')
ax2.tick_params(axis='y', labelcolor='tab:purple')
fig.suptitle("Evolution Toward Target Sequence")
fig.tight_layout()
ax2.set_ylim(0, 20)
fig.legend(loc='upper right')
plt.grid(True)
plt.draw()
# Keyboard controls
paused = False
quit_sim = False
def on_key(event):
global paused, quit_sim
if event.key == ' ': # Spacebar
paused = not paused
status = "PAUSED" if paused else "RUNNING"
fig.suptitle(f"Evolution Toward Target Sequence - {status}\n(Press SPACE to pause, Q to quit)", fontsize=10)
plt.draw()
elif event.key.lower() == 'q': # Q key
quit_sim = True
plt.close()
fig.suptitle("Evolution Toward Target Sequence\n(Press SPACE to pause, Q to quit)", fontsize=10)
fig.canvas.mpl_connect('key_press_event', on_key)
best_seq = ""
best_score = -1
found = False
# Evolution loop
for gen in range(generations):
if quit_sim:
break
# Pause handling
while paused and not quit_sim:
plt.pause(0.1)
fitnesses = [fitness(ind) for ind in population]
total_fit = sum(fitnesses)
best = max(fitnesses)
if(best == 1 and not found):
found = True
print("Found perfect solution at generation", gen)
avg_fitness.append(sum(fitness(ind) for ind in population) / population_size)
if gen % sample_interval == 0:
sample = random.sample(population, sample_size)
avg_beneficial.append(sum(count_beneficial_mutations(ind) for ind in sample) / sample_size)
diversity_over_time.append(diversity(population))
# Update plot data (only plot the sampled data points)
x_sample = [i * sample_interval for i in range(len(avg_beneficial))]
x_all = list(range(len(avg_fitness)))
line1.set_data(x_all, avg_fitness)
line2.set_data(x_sample, avg_beneficial)
line3.set_data(x_sample, diversity_over_time)
ax1.relim(); ax1.autoscale_view()
ax2.relim(); ax2.autoscale_view()
best = max(population, key=fitness)
fig.suptitle(f"Best: {best} (target: {target_sequence})\n(Press SPACE to pause, Q to quit)", fontsize=10)
plt.pause(0.001)
# Tournament selection (as in evolving_fitness_final.py)
new_pop = []
num_competitors = 10 # can be adjusted
for _ in range(population_size):
# Select num_competitors individuals randomly
competitors = random.sample(population, num_competitors)
# Pick the one with highest fitness
fitness_values = [fitness(ind) for ind in competitors]
total = sum(fitness_values)
probs = [f / total for f in fitness_values]
winner = random.choices(competitors, weights=probs, k=1)[0]
# Mutate winner to produce offspring
offspring = mutate(winner)
new_pop.append(offspring)
population = new_pop
if gen % 250 == 0:
best = max(population, key=fitness)
best_individuals.append((gen, best))
plt.show()Answer the following questions using the options available in the model:
Exercise 2.2 (Evolving DNA - Biology / abstract thinking)
- Run the code. What does the new (dashed blue) line represent? Do you understand how it changes over time?
The program reports after how many generations it manages to find the target sequence. With default settings this can take a long time… so you may never see the print
- Try increasing the mutation rate to see how it affects the time to find the target sequence. Try mutation rates between 0.0001 and 0.1. Keep track of both how long (number of generations) it takes to find the target, and how fit the population is once the target is found. What do you observe?
- Diversity is no longer calculated as the standard deviation in fitness, but as the Shannon diversity of all present sequences (although it is not super complex, you do not need to fully understand this function). Because of this, the exact number (quantities) cannot be compared to our earlier model. However, do you see a qualitative differences?
- Study how fitness is calculated in this model. Is there a distinction between genotype en phenotype? Why/why not?
2.5 Evolving a protein sequence
Next we will extend the simulation a little more. The individual genotypes will still be represented as a DNA sequence, but before evaluating fitness this will be translated into a protein sequence. To do so, the code first defines the codon table (which we of course all know by heart =)), and then translates the DNA sequence into a protein sequence. The protein sequence is then used to calculate the fitness of the individual, which is based on how well the protein sequence matches a target protein sequence. The code is as follows:
import random
import math
import matplotlib.pyplot as plt
from collections import Counter
# set the random number seed
random.seed(0)
plt.ion() # Enable interactive mode
# Codon table
codon_table = {
'TTT': 'F', 'TTC': 'F', 'TTA': 'L', 'TTG': 'L', 'CTT': 'L', 'CTC': 'L', 'CTA': 'L', 'CTG': 'L',
'ATT': 'I', 'ATC': 'I', 'ATA': 'I', 'ATG': 'M',
'GTT': 'V', 'GTC': 'V', 'GTA': 'V', 'GTG': 'V',
'TCT': 'S', 'TCC': 'S', 'TCA': 'S', 'TCG': 'S', 'AGT': 'S', 'AGC': 'S',
'CCT': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P',
'ACT': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T',
'GCT': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A',
'TAT': 'Y', 'TAC': 'Y', 'CAT': 'H', 'CAC': 'H',
'CAA': 'Q', 'CAG': 'Q', 'AAT': 'N', 'AAC': 'N',
'AAA': 'K', 'AAG': 'K', 'GAT': 'D', 'GAC': 'D',
'GAA': 'E', 'GAG': 'E', 'TGT': 'C', 'TGC': 'C',
'TGG': 'W', 'CGT': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'AGA': 'R', 'AGG': 'R',
'GGT': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G',
'TAA': '*', 'TAG': '*', 'TGA': '*'
}
# Parameters
alphabet = "ATCG"
target_protein = "DARWIN"
dna_length = len(target_protein)*3
target_length = len(target_protein)
# Simulation settings
population_size = 625 # must be a square number for grid mode
generations = 20000
mutation_rate = 0.0005
sample_interval = 5
sample_size = population_size
no_reproduction_chance = 0.01
# Core functions
def translate(dna):
return ''.join(codon_table.get(dna[i:i+3], '?') for i in range(0, len(dna), 3))
def fitness(dna):
protein = translate(dna)
matches = 0
for i in range(len(protein)):
if protein[i] == target_protein[i]:
matches += 1
return matches / target_length
def mutate(dna, rate=mutation_rate):
return ''.join(
random.choice([b for b in alphabet if b != base]) if random.random() < rate else base
for base in dna
)
def count_beneficial_mutations(dna, num_tests=100):
"""Test random mutations (faster than testing all possibilities)"""
f0 = fitness(dna)
count = 0
for _ in range(num_tests):
i = random.randint(0, len(dna) - 1)
b = random.choice([x for x in alphabet if x != dna[i]])
mutant = dna[:i] + b + dna[i+1:]
if fitness(mutant) > f0:
count += 1
# Extrapolate: if 10 out of 100 random tests were beneficial,
# estimate ~(10/100) * (len(dna) * 3 bases) = total beneficial mutations
total_possible = len(dna) * 3
return count * (total_possible / num_tests)
def diversity(pop):
counts = Counter(pop)
total = len(pop)
return -sum((c/total) * math.log(c/total + 1e-9) for c in counts.values()) if total > 0 else 0
# Initialize population
initial_sequence = ''.join(random.choice(alphabet) for _ in range(dna_length))
population = [initial_sequence for _ in range(population_size)]
avg_fitness = []
avg_beneficial = []
diversity_over_time = []
best_individuals = []
# Grid setup
side = int(math.sqrt(population_size))
assert side * side == population_size, "Population size must be a square number for grid mode"
def get_neighbors(i, j):
return [(x % side, y % side)
for x in range(i-1, i+2)
for y in range(j-1, j+2)]
# Initialize interactive plot
fig, ax1 = plt.subplots(figsize=(12, 8))
ax1.set_xlabel("Generation")
ax1.set_ylabel("Average Fitness", color='tab:blue')
ax1.set_ylim(0, 1)
line0, = ax1.plot([], [], color='black', linewidth=1, label='Max fitness')
line1, = ax1.plot([], [], color='tab:blue', linewidth=2, label='Fitness')
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax2 = ax1.twinx()
ax2.set_ylabel("Beneficial Mutations / Diversity", color='tab:purple')
line2, = ax2.plot([], [], color='tab:purple', linestyle='--', linewidth=2, label='Beneficial Mutations')
line3, = ax2.plot([], [], color='tab:green', linestyle=':', linewidth=2, label='Diversity')
ax2.tick_params(axis='y', labelcolor='tab:purple')
fig.suptitle("Evolution Toward " + str(target_protein))
fig.tight_layout()
ax2.set_ylim(0, 5)
fig.legend(loc='upper right')
plt.grid(True)
plt.draw()
# Keyboard controls
paused = False
quit_sim = False
def on_key(event):
global paused, quit_sim
if event.key == ' ': # Spacebar
paused = not paused
status = "PAUSED" if paused else "RUNNING"
fig.suptitle(f"Evolution Toward {target_protein} - {status}\n(Press SPACE to pause, Q to quit)", fontsize=10)
plt.draw()
elif event.key.lower() == 'q': # Q key
quit_sim = True
plt.close()
fig.suptitle(f"Evolution Toward {target_protein}\n(Press SPACE to pause, Q to quit)", fontsize=10)
fig.canvas.mpl_connect('key_press_event', on_key)
best_seq = ""
best_score = -1
best_fitnesses = []
found = False
# Evolution loop
for gen in range(generations):
if quit_sim:
break
# Pause handling
while paused and not quit_sim:
plt.pause(0.1)
fitnesses = [fitness(ind) for ind in population]
total_fit = sum(fitnesses)
best = max(fitnesses)
best_fitnesses.append(best)
if(best == 1 and not found):
found = True
print("Found perfect solution at generation", gen)
if gen % sample_interval == 0:
sample = random.sample(population, sample_size)
avg_beneficial.append(sum(count_beneficial_mutations(ind) for ind in sample) / sample_size)
diversity_over_time.append(diversity(population))
# Update plot data
line0.set_data(range(len(avg_fitness)+1), best_fitnesses)
line1.set_data(range(len(avg_fitness)+1), avg_fitness + [sum(fitnesses)/population_size])
line2.set_data(range(len(avg_beneficial)), avg_beneficial)
line3.set_data(range(len(diversity_over_time)), diversity_over_time)
ax1.relim(); ax1.autoscale_view()
ax2.relim(); ax2.autoscale_view()
best = max(population, key=fitness)
fig.suptitle(f"Best: {translate(best)} (target: {target_protein})\n(Press SPACE to pause, Q to quit)", fontsize=10)
plt.pause(0.01)
else:
avg_beneficial.append(avg_beneficial[-1])
diversity_over_time.append(diversity_over_time[-1])
new_pop = []
num_competitors = 10 # can be adjusted
for _ in range(population_size):
# Select num_competitors individuals randomly
competitors = random.sample(population, num_competitors)
# Pick the one with highest fitness
fitness_values = [fitness(ind) for ind in competitors]
total = sum(fitness_values)
probs = [f / total for f in fitness_values]
winner = random.choices(competitors, weights=probs, k=1)[0]
# Mutate winner to produce offspring
offspring = mutate(winner)
new_pop.append(offspring)
population = new_pop
avg_fitness.append(sum(fitness(ind) for ind in population) / population_size)
if gen % 250 == 0:
best = max(population, key=fitness)
best_individuals.append((gen, translate(best)))
plt.show()Answer the following question about the model:
Exercise 2.3 (Evolving protein sequences - Biology / abstract thinking)
- Another line was added to the model. What new information can you obtain from analysing this line?
- Study carefully how the other lines (also present in previous models) change over time. What do you observe? Try and capture what you see into words.
- In biology, multiple genotypes can translate to the same phenotype (this is called a many-to-one genotype-phenotype map), or alternatively, one genotype can produce multiple phenotype (this is called phenotypic plasticity, or a one-to-many genotype-phenotype map). Which genotype-phenotype (GP) mapping applies to this model? Why?
2.6 Part 4: RNA structure evolution
Now we’ll take the genotype-phenotype map even further. RNA sequences fold into 3D structures that determine function. The genotype is a sequence, the phenotype is a structure, and fitness depends on how well that structure does its job.
First, you need to install the ViennaRNA package. This is a real bioinformatics tool used by researchers to determine how RNA folds.
Installing ViennaRNA
What is pip? Pip is the Python package manager, it downloads and installs Python libraries. You already have it if you have Python.
Install ViennaRNA: Open a terminal and run:
pip install ViennaRNAThis will not take super long. If it fails, try pip3 install ViennaRNA instead.
Test if it works
Copy this code into a Python file and run it:
from RNA import fold
rna_sequence = "GCGCUUCGAGAGAGGATAGCATCGCCG"
structure, energy = fold(rna_sequence)
print("Sequence:", rna_sequence)
print("Structure:", structure)
print("Energy:", energy)You should see something like:
Sequence: GCGCUUCGAGAGAGGATAGCATCGCCG
Structure: ((.((((....))))...)).......
Energy: -4.599999904632568The dots and parentheses are a simple way to show the 2D fold of an RNA molecule. Brackets show the maching base pairs, and dots are unpaired bases.
Exercise 2.4 (Evolving RNA structures - Coding practice) Use the same evolution framework as in Parts 2 and 3, but instead of comparing sequences or proteins, compare RNA structures. Fitness = how brackets/dots are matching, rather than comparing the nucleotides directly. For example, if your target is ((((....)))) and your RNA folds to ((((....)))), fitness = 1.0. If it matches 10 out of 12 positions, fitness = 10/12 ≈ 0.833, etc.
Start with a random RNA sequence of length 32, and try to evolve it to match a target secondary structure. To find a possible target structure, you can fold another random sequence, and define this as the end-target.
Try alternating between two different target structures, and describe what you observe (if anything, I personally don’t know if anything exciting happens).
Hint: Structure comparison is quite simple:
def structure_fitness(rna_seq, target_structure):
structure, _ = fold(rna_seq)
# Count how many positions match the target
matches = 0
for i in range(len(structure)):
if structure[i] == target_structure[i]:
matches += 1
return matches / len(target_structure)Note that the code is getting longer and longer. Later in the course, we will also show you how code can be organised across files/folders↩︎