Be an Epidemiologist! Challenge

You must create a simple model of how disease spreads around a group of people.

Rules and Requirements

The model must be a 1000 by 1000 2D array with each element being a different person.

The user must input three variables using argv: the probability of transmission (how likely someone is going to infect someone else), chance of mutation and how many periods the simulation should run for.

In the first period (t=0), four people should be chosen randomly and infected with the disease.

The way the disease behaves is governed by the following rules:


The disease can only move vertically and horizontally, moving to the person next door.
The infection last for 3 periods in every person. You may not factor in immunodeficiencies.
After a person has been infected three times, they are immune and cannot be infected again. 
The disease is subject to mutation making the previously immune people vulnerable to this new mutated disease. The mutated disease has exactly the same traits and follows the same rules as the original disease.
If a mutation occurs, the whole disease does not change, just that particular 'packet' upon transmission.
Once a person has been been infected by one virus, they cannot be infected again until the current infection passes.
If a person is infected, they are infectious from the start of their infection period to its end.
There are no levels of immunity - a person is either immune or not.
To stop memory overload, there is a maximum limit of 800 mutations.


At the end of the specified number of periods, you should output the results.

The results must be a 1000 x 1000 grid showing which people are infected and which people are not. This can be output as a text file, as an image file or graphical output (where #FFFFFF is a healthy person and #40FF00 is an infected person).

Please can you include the language name and a command to run it in your answer.

Winning

The fastest code to run on my computer wins. Its time will be measured with the following piece of Python code:

import time, os
start = time.time()
os.system(command)
end = time.time()
print(end-start)


Note that when running this script, I will use the following defaults:

Probability of transmission = 1
Chance of mutation = 0.01
Number of periods = 1000
I was curious what this would look like so I made this quick and dirty hack in JavaScript: http://jsfiddle.net/andrewmaxwell/r8m54t9c/

// The probability that a healthy cell will be infected by an adjacent infected cell EACH FRAME.
var infectionProbability = 0.2

// The probability that the infection will mutate on transmission EACH FRAME.
var mutationProbability = 0.00001

// The maximum number of times a cell can be infected by the same infection.
var maxInfections = 3

// The number of frames a cell in infected before it becomes healthy again.
var infectionDuration = 3

// The width and heigh of the board
var size = 400

// The number of cells infected at the beginning.
var startingNum = 4

var imageData, // the visual representation of the board
    cells, // array of cells
    infectionCount // counter that is incremented whenever a mutation occurs

// Just some colors. The colors are re-used as the number of mutations increases.
var colors = [[255,0,0],[255,255,0],[0,255,0],[0,255,255],[0,0,255],[255,0,255],[128,0,0],[128,128,0],[0,128,0],[0,128,128],[0,0,128],[128,0,128],[255,128,128],[255,255,128],[128,255,128],[128,255,255],[128,128,255],[255,128,255]
]

// when a cell is infected, it isn't contagious until the next frame
function infect(person, infection){
    person.infect = true
    person.infectionCounts[infection] = (person.infectionCounts[infection] || 0) + 1
    person.currentInfection = infection
}

// when a mutation occurs, it is given a number and the counter is incremented
function mutation(){
    return infectionCount++
}

function reset(){

    cells = []
    infectionCount = 0
    imageData = T.createImageData(size, size)

    // initialize the cells, store them in a grid temporarily and an array for use in each frame
    var grid = []
    for (var i = 0; i &lt; size; i++){
        grid[i] = []
        for (var j = 0; j &lt; size; j++){
            cells.push(grid[i][j] = {
                infectionTime: 0, // how many frames until they are no longer infected, so 0 is healthy
                infectionCounts: [], // this stores how many times the cell has been infected by each mutation
                neighbors: [] // the neighboring cells
            })
        }
    }

    // store the neighbors of each cell, I just want to minimize the work done each frame
    var neighborCoords = [[0,-1],[1,0],[0,1],[-1,0]]
    for (var i = 0; i &lt; size; i++){
        for (var j = 0; j &lt; size; j++){
            for (var n = 0; n &lt; neighborCoords.length; n++){
                var row = i + neighborCoords[n][0]
                var col = j + neighborCoords[n][1]
                if (grid[row] &amp;&amp; grid[row][col]){
                    grid[i][j].neighbors.push(grid[row][col])
                }
            }
        }
    }

    // infect the initial cells
    for (var i = 0; i &lt; startingNum; i++){
        infect(cells[Math.floor(cells.length * Math.random())], 0)
    }
}

function loop(){
    requestAnimationFrame(loop)

    // for each cell marked as infected, set its infectionTime
    for (var i = 0; i &lt; cells.length; i++){
        var p = cells[i]
        if (p.infect){
            p.infect = false
            p.infectionTime = infectionDuration
        }
    }

    for (var i = 0; i &lt; cells.length; i++){
        var p = cells[i]

        // for each infected cell, decrement its timer
        if (p.infectionTime){
            p.infectionTime--

            // for each neighbor that isn't infected, if the probability is right and the neighbor isn't immune to that infection, infect it
            for (var n = 0; n &lt; p.neighbors.length; n++){
                var neighbor = p.neighbors[n]
                if (!neighbor.infectionTime &amp;&amp; Math.random() &lt; infectionProbability){
                    var infection = Math.random() &lt; mutationProbability ? mutation() : p.currentInfection
                    if (!neighbor.infectionCounts[infection] || neighbor.infectionCounts[infection] &lt; maxInfections){
                        infect(neighbor, infection)
                    }
                }
            }

            // colors! yay!
            var color = colors[p.currentInfection % colors.length]
            imageData.data[4 * i + 0] = color[0]
            imageData.data[4 * i + 1] = color[1]
            imageData.data[4 * i + 2] = color[2]
        } else {
            imageData.data[4 * i + 0] = imageData.data[4 * i + 1] = imageData.data[4 * i + 2] = 0
        }

        imageData.data[4 * i + 3] = 255
    }

    T.putImageData(imageData, 0, 0)
}

// init canvas and go
C.width = C.height = size
T = C.getContext('2d')
reset()
loop()