Nature, in Code

Balancing Selection

Given the equation of the change in allele frequency per generation:

$$\Delta p = \frac{pqs ( ph + q (1-h))}{(1 - 2pqhs - q^{2}s)}$$
balancing selection occurs when the heterozygous effect, h, is negative. In other words, the fitness of the heterozygote A1A2 is higher than the fitness of both homozygotes. When plotting the change of the allele A1 frequency, Δp, as a function of the allele's current frequency, we'll get the following graph (with s = 0.1, h = -0.5):

Because Δp is positive when p is small, and negative when p is large, p will go to a stable equilibrium value at

$$p* = \frac{h-1}{2h - 1}$$
We can see this when we run multiple simulations with variable initial values of p:

Code

First plot

var s = 0.1;
var h = -0.5;

var data=[];
var x_max = 1;

for (var i = 0; i <= x_max + 0.005; i = i + 0.01) {
    var p = i;
    var q = 1-p;
    var delta_p = (p*q*s * ( p*h + q*(1-h))) / (1 - 2*p*q*h*s -q*q*s);
    data.push(delta_p);
}

draw_line_chart(data,"p","\u0394 p",[],x_max,true);
			
Second plot

var s = 0.1;
var h = -0.5;
var initial_p = 0.01;
var p_values = [];

var data = [];
var generations = 300;

while (initial_p < 1) {
	p_values.push(initial_p);
	initial_p = initial_p + 0.01;
}

for (var i = 0; i < p_values.length; i = i + 1) {
	run_simulation(p_values[i]);
}

function run_simulation(p) {
	var results = [];
	for (var i = 0; i < generations; i = i + 1) {
		var q = 1-p;
		var delta_p = (p*q*s * ( p*h + q*(1-h))) / (1 - 2*p*q*h*s -q*q*s);
		p = p + delta_p;
		results.push(p);
	}
	data.push(results);
}

draw_line_chart(data,"Generation","p",[],generations);
			
Note: the draw_line_chart function is built with D3.js and can be found here.