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<r 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<R<1. Generate a random number x and if x<R 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<B<1. Generate a random number x and if x<D 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.

Screen Shot 2015-08-04 at 11.52.27

 

 


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()

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s