Nature, in Code

Disruptive 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)}$$
disruptive selection occurs when the heterozygous effect, h, is larger than 1. In other words, the fitness of the heterozygote A1A2 is lower 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 = 1.5):

Because Δp is negative when p is small, and positive when p is large, p will either go do 0 or 1, depending on whether its initial value is small or greater than

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

p* is an unstable equilibrium. Even though Δp is 0 at p*, the slightest deviation from p* will push the system to the stable equilibria p = 0 or p = 1:

Code

First plot

var s = 0.1;
var h = 1.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 = 1.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);
			

Third plot

var s = 0.1;
var h = 1.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 + (Math.random() * 0.000001);
}

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.