Genetic Diffusion Across a Geographic Barrier

To contact us Click HERE
Imagine some red flowers growing in a meadow at the base of a mountain range. Suppose one of the plants has a mutation that makes its flowers blue. If the new color makes it easier for insects to find and pollinate it, then as the years go by, the mutation will likely spread and replace the original red color in nearby plants. But will the mutation be able to spread to the far side of the mountain range? How long will it take before you find blue flowers there?

A lot of readers have enjoyed An Evolving Ecosystem Game, so today we are going to modify that agent-based simulation of an evolving ecosystem to model the diffusion of mutations across a geographic landscape with a barrier.

Let's start by describing the simulation conceptually.

We will represent land by a large grid of squares, 100 by 100.We can think of each of the 10,000 cells in our grid as having asingle plant in it. Or, for a more realism, maybe think of each cellas a square mile (or square kilometer), and assume that each cellholds a relatively homogeneous population of many plants.

In our first experiments, all the cells will hold plants. Inlater experiments, we will block off some of the cells to form abarrier, representing the inhospitable mountain range.

Each cell (plant, or "agent") has a color gene that we willrepresent as a number between zero and one, corresponding to aspectrum of colors from red to blue.

Initially, all the plants have red flowers, so we will color the cellsred. As time goes on, the plant in a given cell dies and is replacedwith a new plant grown from seed.

The color of the new plant will be inherited from its parents,which are the plant it replaced and a randomly chosenneighbor. Initially, we will assume color has no impact on fitness, sothat all neighbors are equally likely to be the parent, but later, wewill assume that colors closer to the blue end are favored bypollinating insects over colors at the red end.

Initially, I intended to let the new color be the average ofthe two parent colors, plus a small random component representingmutations that can shift it either toward red or toward blue.

However, averaging quickly wipes out the mutation, even whenit is favorable. The gene with value 1 (the initial blue plant) mixeswith one of the nearby 0's to make 0.5, which then mixes with anothernearby 0 to make 0.25, and so on. Within a few generations, everythingis back to red. This is not the situation I wanted to model.

Instead, we let the new color be inherited from one of the twoparents, chosen at random (50-50). And, with a small probability, wealso allow a mutation that shifts the gene up or down by 0.1 (thoughwe also bound it to stay between zero and one).

Implementation in Python

So - ready to try it out? As usual, today we will code up theexperiment using free, open-source software. We will once again usethe free, open-source, cross-platform programming language calledPython, which is one of the mostwidely used and easy to learn programming languages in the world.

Start by downloading and installing Python on your machine via the Python 3.2.2 Windows installer. You will also need thePyGame add-on package,usingthe PyGameWindows installer for Python 3.2, which is also free and opensource. Use the version of PyGame that most closely matches theversion of Python you are using.

My article on An Evolving Ecosystem Game has some additional pointers on getting started with Python and PyGame. I also highly recommend reading a tutorial such as the free book-length PDF file Introduction to Computer Science Using Python and PyGame.

Start Python. From the File menu, choose "New Window". This brings upa window you can use for editing Python scripts. Now copy and pastethe Python code listing below into your new editor window. Save theresulting file with a name like "flowers.py", and put it somewhere youcan find it again, such as on your desktop.

import pygameimport randomimport mathfavorable = Falsebarrier = FalsemutationRate = 0.01name = "neutral"makeMovie = Falsesize = 100 # boxes per side of gridpixels = 4 # pixels per side of boxdef breed(A):    if(A.color < 0):        return A    a = Agent(A.i, A.j)    a.color = A.color    if(random.random() < 0.5):        a.color = A.chooseNeighbor().color    if(random.random() < mutationRate):        diff = math.copysign(0.1, random.random() - 0.5)        a.color = max(0, min(1, a.color + diff))    return aclass Agent:    def chooseNeighbor(self):        a = self.randomNeighbor()        if(favorable):            b = self.randomNeighbor()            if(a.color > b.color):                return a            else:                return b        return a    def randomNeighbor(self):        i = bound(self.i + random.choice([-1,0,1]))        j = bound(self.j + random.choice([-1,0,1]))        if (i==self.i and j==self.j) or pop[i][j].color < 0:            return self.randomNeighbor() # try again        else:            return pop[i][j]    def __init__(self, i, j):        self.i = i        self.j = j        self.color = 0def bound(x):    return max(0, min(size-1, x))def rainbow(x):    z = pygame.Color(0,0,0,0)    if(x < 0):        return z    z.hsva = [x*250, 90, 90, 0]    return zslide = 1def draw(pop):    global slide    for row in pop:        for cell in row:            R = pygame.Rect((cell.i*pixels, cell.j*pixels), (pixels, pixels))            pygame.draw.rect(screen, rainbow(cell.color), R, 0)    if makeMovie:        pygame.image.save(screen, "img/slide%(n)03d.jpg" % {"n":slide})    slide += 1    pygame.display.flip()def update(oldpop):    pop = [] # new population matrix    for oldrow in oldpop:        row = []        pop.append(row)        for oldcell in oldrow:            row.append(breed(oldcell))    return pop# initialize graphics window:pygame.init()nPix = pixels*sizescreen = pygame.display.set_mode([nPix, nPix])pygame.display.set_caption("Genetic Diffusion")white = [255, 255, 255]black = [0, 0, 0]screen.fill(white)pygame.display.flip()# initialize agents:pop = [] # population matrixfor i in range(0,size):    row = []    pop.append(row)    for j in range(0,size):        row.append(Agent(i, j))if(barrier):    b = int(size/2)    for i in range(1,size-1):        pop[i][b].color = -1    for i in range(0,4):        for j in range(0,4):            pop[i][j].color = 1# main animation loop:gen = 1done = Falseclock = pygame.time.Clock()while done==False:    draw(pop)    gen += 1    clock.tick(30)    for event in pygame.event.get():        if event.type == pygame.QUIT:            done = True    pop = update(pop)pygame.image.save(screen, name+".jpg")pygame.quit()

Now, press F5 in the editor window to run the code (or go to the Runmenu if you prefer, and choose "Run Module"). The graphics windowopens and you can watch the simulation evolve. If you are using Linux,just type

python flowers.py

at the command line instead.

You should see a graphics window. Initially it will be all red, butsoon, orange and then yellow dots will appear. As time goes on, areaswith a lot of yellow will start to turn green, and if you keepwaiting for a minute or so, you will even start to see some cyan(light blue) and dark blue emerging. Ultimately, you will wind up witha picture something like this. Because of different random dice rolls,your pattern will differ in specific details, but should look prettysimilar in a qualitative sense: localized patches of all the differentcolors, all co-existing, because in this first run, weset favorable = False, so no one color is preferred.

If the code does not work for you, the first thing tocheck is whether the indentation matches what you see above:in Python, unlike in most other programming languages, indentationmatters. Python uses indentation to group statements thatlanguages like C or Java would surround with curly braces.

With this as a baseline, let's see what happens if we changethe first line of code (after the import statements) from favorable= False to favorable = True. Also, change name ="neutral" to name = "favorable". First, let's look at thecode so we know what to expect. After the "import" lines, which justload the Python packages we are using (PyGame, random numbers, and themath library), we set various global variables like the mutation rateand the size of the grid.

The breed function is central to the simulation, so we willcome back to it in a moment.

Since each agent only needs to keep track of a single number(the color gene), we don't really need the programminginfrastructure of agent based modeling here - we could just store amatrix of colors and be done with it. However, I have written the codein a more agent-based style so that it will be easier for you to use it asa starting point for making more complicated simulations.

By defining an Agent class, we create a new data type, likenumbers and strings, called "Agent". These are pretty simple agents:they have coordinates (i and j) showing where they live on the grid,and a color that starts at zero. They also have a method (function)for selecting a random neighbor, and a more interesting methodcalled chooseNeighbor. This is the function that breedwill use. It picks a random neighbor, and then if favorable istrue, it picks a second random neighbor and of the two, returns theneighbor with the largest color value.

The breed function starts with an agent called A, and checkswhether its color is negative, which means it is a black barrier cellrather than a colored plant cell. Assuming it is a real plant, webreed it by choosing a neighbor to be the second parent, theninheriting the color of one of the parents. Occasionally (based on themutation rate) we also shift the color up or down by 0.1. Finally, wereturn the new agent (plant) we have just grown from seed.

In this game, the agents are synchronized: every winter they all dieand are replaced by a whole new generation. This means we need twomatrices: one to hold the old generation and one to hold the new one,so that when we compute the parent colors for a new plant, we readfrom the old generation, not from a partial mixture of old and new.That's what the update function does, a little further down inthe code.

We also have a rainbow function for mapping the zero toone color gene into an actual Color on screen. And a drawfunction for re-drawing the screen once per new generation.

The rest of the code just puts an initial plant (agent) ineach cell of the grid, then possibly adds a barrier (mountain range),then finally runs the main loop, in which we repeatedly compute a newgeneration of plants and draw them on the screen. When you finallyclick the Close box on the graphics window, we save the final image toa JPG graphics file for future reference. For instance, we made onecalled "neutral.jpg" from the first run. When you re-run the code, youwill create "favorable.jpg", as shown below.

Actually, I stopped the simulation early - after a few moreseconds the screen would be solid dark blue. Why the dramaticdifference from the first experiment? Clearly, now that colors closeto blue are favored over colors close to red, they tend to spread. Wehave broken the symmetry, so the various colors no longer coexist inequilibrium with each other.

For out third and final experiment, set up the first few lines likethis:

favorable = Truebarrier = TruemutationRate = 0name = "barrier"

I also set makeMovie to True, and created a folder called "img"in the same folder where I have the code. This makes the code write aJPG image file for each generation, which I then combined togetherinto an animated movie. For anyone not actually running the code, youcan play the following movie to see what happens.

As you can see, we are now adding the mountain range barrier (thehorizontal black line in the middle of the picture). Up in the topleft, we create an initial mutation (in blue). In order to focus onthe spread of the initial mutation, we set the mutation rate to zeroso there will not be any others: the only way to get blue squares onthe far side of the mountain will be if the original mutation spreads.

While the screen does ultimately fill with blue, the path to gettingthere is quite different from our second experiment. This time, thegeographical barrier plays a key role in breaking the symmetry.

In case you cannot see the movie (modern browsers like Firefox shouldwork fine, but older ones from Microsoft may fail), here are threepictures from early, middle and late in the simulation, to give youthe idea.

The first time I created an animated movie from the results of asimulation was back in the early 1990's. At that time, I had to fly toCalTech inPasadena, California to use a super-computer lab with highlyspecialized (and highly expensive) equipment. The process took acouple of days (not to mention thousands of dollars). Two decadeslater, we can all do this sort of thing in seconds on our ordinaryhome computers!

I have fond memories of that trip to CalTech. One verymemorable highlight was wandering the wonderful10-acre DesertGarden at the nearby Huntington Library. If you ever have theopportunity, I highly recommend it. Coming originally from moistnorthern forests, I found it an extremely moving experience to be inthe dry desert, out of sight and sound of civilization, walking pasttowering, majestic cacti. Like seeing the Milky Way galaxy on a darknight in the deep woods, it is an experience that helps you understandyour place in the universe.

I hope you enjoyed this discussion. I encourage you to try your ownexperiments. If you change something and it fails to run, read theerror message; it often gives a good clue about what line the problemis in (usually the line you just changed), and indentation is oftenthe issue. It is best to use an editor that understands Python, suchas the one built into IDLE, since it can color code the text and helpmanage indentation for you. For instance, text following a # isignored, so we can put comments in the code to remind us whatit does.

You can click the "M"button below to email this post to a friend, or the "t" button toTweet it, or the "f" button to share it on Facebook, and so on.

As usual, please post questions, comments and other suggestions using the box below, or G-mail me directly at the address mentioned in the Welcome post. Remember that you can sign up for email alerts about new posts by entering your address in the widget on the sidebar. If you prefer, you can follow @ingThruMath on Twitter to get a 'tweet' for each new post. The Contents page has a complete list of previous articles in historical order. You may also want to use the 'Topic' and 'Search' widgets in the side-bar to find other articles of related interest. See you next time!