A toy model for evolution, of a “unicellular mathematical organism”

Here is an idea, with which I wanted to play with for a little while.

Let us look at the evolution of a “unicellular mathematical organism”. We’ll start to say that each organism is defined by its “mathematical DNA” – Here we are going to consider 2 pairs of integers (or chromosomes) only. The first pair is defined by 2 integers $a$ and $a'$. We will call it the $A$-pair. The second pair, call it the $B$-pair, is defined by $b$ and $b'$. For simplicity we will define $A=(a,a')$ and $B=(b,b')$.

We will say that $A$ is the pair of genes which govern the reproduction, while $B$ is the pair of genes which defines their fitness, and ultimately their death.

It follow that a cell is define by the couple of pairs $(A,B)=((a,a'),(b,b'))$ .

Now here are the rules which govern  the reproduction, birth and death of these cells.

1 – Reproduction and Birth:

When two cells meet [cells $1$ and $2$, with genomes $(A_1, B_1 )$ and $(A_2, B_2)$], they can reproduce (generate a new cell) with probability $R(A_1,A_2)$$R$ is a function of the pairs $A_1$ and $A_2$.

The genome of the new cell is formed by two halfs of the genome of its parents. It can be any of the following outcome:

$((a_1,a_2),(b_1,b_2))$$((a_1,a'_2),(b_1,b_2))$$((a_1,a'_2),(b_1,b'_2))$$((a'_1,a_2),(b_1,b_2))$, … and all other possibilities.

2 – Death

We will simply say that at any time a cell can die with probability by time unit $D(B)$.

Defining the rates $R$ and $D$.

The reproduction rate $R$ is function of the $A$-genes of the two organisms. One option is to choose

$R(A_1,A_2) = (a_1+a'_1+a_2+a'_2)/K$ (for some constant $K$), so that the larger are the integers of each cell, the most likely they are to reproduce. One possible alternative solution would be $R(A_1,A_2)=\min(a_1+a'_1)+\min(a_2+a'_2)$ or the same using $\max$ instead of $\min$. This would be equivalent as introducing a dominant or recessive gene.

Identically, $D$ could be define by $D(B)=b+b'$. In this way cell with small integers on the $B$-gene would be most likely to survive and therefore to transmit their genes to the next generations.

A first algorithm:

Start with $N_0$ initial cells. Each cell has genome $((a,a'),(b,b'))$ where $a,a',b,b'$ are integers randomly distributed between $1$ and $P$.

Let $r,d$ be two real numbers such that $r+d=1$. With probability $r$ we are going to select two cells for the reproduction process. With probability $d=1-r$ we are going to select a cell and test its survival skills. So we generate a random number $x$ between $0$ and $1$. If $x try the reproduction process, else try the death process.

1 – Reproduction Process:

Select two cells and calculate the probability $R(A_1,A_2)$ one could choose

$R=\frac{a_1+a'_1+a_2+a'_2}{4P}$, so that $0. Generate a random number $x$ and if $x the reproduction process is successful and we generate a new cell.

2 – Death Process:

Select one cell and calculate the probability $D(B)$ one could choose

$B=\frac{b+b'}{2P}$, so that $0. Generate a random number $x$ and if $x the cell dies.

Below you see the evolution of the Population size as a function of the time, when choosing parameters: $r=0.3, d=0.7$ and $P=10$, and initial condition $N_0=100$. Note that the evolution is non-monotonic.

THE PYTHON CODE TO PLAY WITH

import sys , math , cmath , random , shelve

import numpy

proba_r = 0.5

proba_d = 0.5

N_step = 140

N_specie_0 = 100       # Number of initial species

class organism() :        # we here define the class which will generate organisms as objects

def __init__(self):

a = random.randint( 1 , 10 )

ap = random.randint( 1 , 10 )

b = random.randint( 1 , 10 )

bp = random.randint( 1 , 10 )

self.A = [ a , ap ]

self.B = [ b , bp ]

Average_A = []

Average_B = []

Pop_size = []

Tab_of_object = []

for ii in range( N_specie_0 ):

Tab_of_object += [ organism() ]

aaa = 0

bbb = 0

for obj in Tab_of_object :

aaa += obj.A[ 0 ] + obj.A[ 1 ]

bbb += obj.B[ 0 ] + obj.B[ 1 ]

aaa = aaa * 1.0 / (2 * len( Tab_of_object ))

bbb = bbb * 1.0 / (2 * len( Tab_of_object ))

Average_A += [ aaa ]

Average_B += [ bbb ]

Pop_size += [ len( Tab_of_object ) ]

for tt in range( N_step ) :

LL = len(Tab_of_object)

for ppp in range( LL ) :      # a time step is the realisation of LL updates where LL is the number of total organisms

if ( len( Tab_of_object ) == 1 ) : #the code quit is we have only one object left as reproduction is then impossible

print “quiting”

quit()

x = random.random()

if ( x < proba_r)  :

index_1 = random.choice( range( len( Tab_of_object ) ) )

index_2 = index_1

while ( index_2 == index_1 ) :

index_2 = random.choice( range( len( Tab_of_object ) ) )

obj_1 = Tab_of_object[ index_1 ]

obj_2 = Tab_of_object[ index_2 ]

R = ( obj_1.A[0] + obj_1.A[1] + obj_2.A[0] + obj_2.A[1] ) * 1.0 / ( 4.0 * 10.0 )

y = random.random()

if ( y < R ) :

obj_baby = organism()

obj_baby.A[ 0 ] = obj_1.A[ random.choice( [ 0 , 1 ] ) ]

obj_baby.A[ 1 ] = obj_2.A[ random.choice( [ 0 , 1 ] ) ]

obj_baby.B[ 0 ] = obj_1.B[ random.choice( [ 0 , 1 ] ) ]

obj_baby.B[ 1 ] = obj_2.B[ random.choice( [ 0 , 1 ] ) ]

Tab_of_object += [ obj_baby ]

if ( proba_r < x < ( proba_r + proba_d ) ) :

index_of_obj_to_kill = random.choice( range( len( Tab_of_object ) ) )

obj = Tab_of_object[ index_of_obj_to_kill ]

D = ( obj.B[0] + obj.B[1] ) * 1.0 / ( 2.0 * 10.0 )

y = random.random()

if ( y < D ) :

Tab_of_object.pop( index_of_obj_to_kill )

aaa = 0

bbb = 0

for obj in Tab_of_object :

aaa += obj.A[ 0 ] + obj.A[ 1 ]

bbb += obj.B[ 0 ] + obj.B[ 1 ]

aaa = aaa * 1.0 / (2.0 * len( Tab_of_object ) )

bbb = bbb * 1.0 / (2.0 * len( Tab_of_object ) )

Average_A += [ aaa ]

Average_B += [ bbb ]

Pop_size += [ len( Tab_of_object ) ]

import matplotlib.pyplot as plt

plt.plot( Pop_size )

plt.xlabel(‘time step’)

plt.ylabel(‘Pop_size’)

plt.show()

plt.plot( Average_A )

plt.plot( Average_B )

plt.xlabel(‘time step’)

plt.ylabel(‘Average_A Average_B’)

plt.show()