Population dynamic simulation on biological information maintenance

Question

Background

Using this simulation I investigate a system in which enzymes proliferate in cells. During the replications of enzymes, parasites can come to be due to mutation. They can drive the system into extinction. I'm interested in where in the parameter space coexistence is possible.

In the program the system is a list, the cells are dictionaries with 2 keys: "e" for the enzymes and "p" for the parasites. The values of the keys are the numbers of the 2 variants.

Our parameters are:

• pop_size: the number of the cells
• cell_size: the maximal number of molecules (enzymes+parasites) of cells at which cell division takes place
• a_p: fitness of the parasites relative to the fitness of the enzymes (for example if a_p = 2, the parasites' fitness is twice as that of the enzymes)
• mutation_rate: the probability of mutation during a replication event
• gen_max: the maximal number of generations (a generation corresponds to one
• while cycle; if the system extincts, the program doesn't run until gen_max)

We start with pop_size cells with cell_size // 2 enzimes and 0 parasites. In each cell the molecules proliferate until their number reaches cell_size. Each cell divides, the assortment of the molecules happens according to binomial distributions (\$p=0.5\$). Cells with "e" < 2 are discarded as dead. After that if the number of viable cells is bigger than pop_size, we choose pop_size of them according to cell fitness ("e"/("e"+"p")), and they move on to the next generation. On the other hand, if the number of viable cells is pop_size or less, they all move on to the next generation.

My request

I've never studied programming in school. This program is the result of heavy googling. Now I've reached a point where I need advice from experienced people. At certain parameter values the program gets quite slow.

1. What better solutions exist performance-wise than my solutions for the manipulations of the list's items throughout the program and for writing data to file? And algorithm design-wise?

2. In which directions should I improve my programming skills in Python to efficiently implement these kind of models? Or am I near the limit of Python's capabilities in this regard?

3. Should I change to a more appropriate programming language in order to achieve significantly better performance at these kind of tasks? If yes, which languages should I consider? (My guess is C.)

The program consists of two functions. simulation() does the simulation, writeoutfile() writes the data to file.

# -*- coding: utf-8 -*-
from random import choices, random
import csv
import time
import numpy as np


def simulation(pop_size, cell_size, a_p, mutation_rate, gen_max):
    def fitness(pop):
        return [i["e"] / (i["e"] + i["p"]) for i in pop]

    def output(pop, gen, pop_size, cell_size, mutation_rate, a_p, boa_split):
        if pop:
            gyaklist_e = [i["e"] for i in pop]
            gyaklist_p = [i["p"] for i in pop]
            fitnesslist = fitness(pop)
            return (
                gen,
                sum(gyaklist_e), sum(gyaklist_p),
                sum([1 for i in pop if i["e"] > 1]),
                np.mean(gyaklist_e), np.var(gyaklist_e),
                np.percentile(gyaklist_e, 25),
                np.percentile(gyaklist_e, 50),
                np.percentile(gyaklist_e, 75),
                np.mean(gyaklist_p), np.var(gyaklist_p),
                np.percentile(gyaklist_p, 25),
                np.percentile(gyaklist_p, 50),
                np.percentile(gyaklist_p, 75),
                np.mean(fitnesslist), np.var(fitnesslist),
                np.percentile(fitnesslist, 25),
                np.percentile(fitnesslist, 50),
                np.percentile(fitnesslist, 75),
                pop_size, cell_size, mutation_rate, a_p, boa_split
                )
        return (
            gen,
            0, 0,
            0,
            0, 0,
            0, 0, 0,
            0, 0,
            0, 0, 0,
            0, 0,
            0, 0, 0,
            pop_size, cell_size, mutation_rate, a_p, boa_split
            )

    pop = [{"e": cell_size // 2, "p": 0} for _ in range(pop_size)]
    gen = 0
    yield output(
        pop,
        gen, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
        )
    print(
        "N = {}, rMax = {}, aP = {}, U = {}".format(
            pop_size, cell_size, a_p, mutation_rate
            )
        )

    while pop and gen < gen_max:
        gen += 1

        for i in pop:
            while not i["e"] + i["p"] == cell_size:
                luckyreplicator = choices(
                    ["e", "p"], [i["e"], a_p*i["p"]]
                    )
                if luckyreplicator[0] == "e" and random() < mutation_rate:
                    luckyreplicator[0] = "p"
                i[luckyreplicator[0]] += 1

        if gen % 100 == 0:
            yield output(
                pop,
                gen, pop_size, cell_size, mutation_rate, a_p, boa_split="bef"
                )

        newpop = [
            {"e": np.random.binomial(i["e"], 0.5),
             "p": np.random.binomial(i["p"], 0.5)}
            for i in pop
            ]
        for i in zip(pop, newpop):
            i[0]["e"] -= i[1]["e"]
            i[0]["p"] -= i[1]["p"]

        pop += newpop
        newpop = [i for i in pop if i["e"] > 1]

        if newpop:
            fitnesslist = fitness(newpop)
            fitness_sum = np.sum(fitnesslist)
            fitnesslist = fitnesslist / fitness_sum
            pop = np.random.choice(
                newpop, min(pop_size, len(newpop)),
                replace=False, p=fitnesslist
                ).tolist()
        else:
            pop = newpop
            for i in range(2):
                yield output(
                    pop,
                    gen+i, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
                    )
            print("{} generations are done. Cells are extinct.".format(gen))

        if gen % 100 == 0 and pop:
            yield output(
                pop,
                gen, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
                )

        if gen % 1000 == 0 and pop:
            print("{} generations are done.".format(gen))


def writeoutfile(simulationresult, runnumber):
    localtime = time.strftime(
        "%m_%d_%H_%M_%S_%Y", time.localtime(time.time())
        )
    with open("output_data_" + localtime + ".csv", "w", newline="") as outfile:
        outfile.write(
            "gen"+";" +
            "eSzamSum"+";"+"pSzamSum"+";" +
            "alive"+";" +
            "eSzamAtl"+";"+"eSzamVar"+";" +
            "eSzamAKv"+";" +
            "eSzamMed"+";" +
            "eSzamFKv"+";" +
            "pSzamAtl"+";" + "pSzamVar" + ";" +
            "pSzamAKv"+";" +
            "pSzamMed"+";" +
            "pSzamFKv"+";" +
            "fitAtl"+";"+"fitVar"+";" +
            "fitAKv"+";" +
            "fitMed"+";" +
            "fitFKv"+";" +
            "N"+";"+"rMax"+";"+"U"+";"+"aP"+";"+"boaSplit"+"\n"
            )
        outfile = csv.writer(outfile, delimiter=";")
        counter = 0
        print(counter, "/", runnumber)
        for i in simulationresult:
            outfile.writerows(i)
            counter += 1
            print(counter, "/", runnumber)


RESULT = [simulation(100, 20, 1, 0, 10000)]
RESULT.append(simulation(100, 20, 1, 1, 10000))
N_RUN = 2
writeoutfile(RESULT, N_RUN)
# Normally I call the functions from another script,
# these last 4 lines are meant to be an example.

On parameter values

So far combinations of these values were examined:

• pop_size: 100; 200; 500; 1000
• cell_size: 20; 50; 100; 200; 500; 1000
• a_p: 0.75; 1; 1.25; 1.5; 1.75; 2; 3
• mutation_rate: 0-1
• gen_max: 10000

Primarily I would like to increase pop_size and above 1000 cells the program is slower than I would prefer. Of course that's somewhat subjective, but for example a million cells would be a perfectly reasonable assumption and at that order of magnitude I think it's objectively impossibly slow.

The program also gets slower with the increase in cell_size and slightly slower with a_p, but for the time being I'm happy with the values of the former and the effect of the latter is tolerable.

The effect of the mutation rate on speed is also tolerable.

In addition to pop_size, gen_max should be increased and has significant effect on run time. I know I don't catch every extinction events with 10000 generations. 20000 would be better, 50000 would be quite enough and 100000 would be like cracking a nut with a sledgehammer.

Answer 1

Numpy can be extremely fast, nigh on as fast as C or other low level languages (because it uses C!). But this is on the condition that the slow stuff is actually done in Numpy. By which I mean, you can't keep looping through lists and dictionaries then do select actions in Numpy, you have to stick to Numpy arrays and element-wise operations.

I will give some comments on style then return to that.

• First, there are zero comments throughout your entire code. I recommend both """docstrings""" at the start of your functions and short # Comments between lines where code is a little confusing.

• f-strings are a python 3.6+ feature which greatly improve readability. They are used in place of .format() and string concatenation. For example:

print(f'{gen} generations are done. Cells are extinct.')

• You spread a lot of code over several lines when really, longer lines would be cleaner. You don't have very-highly nested code so the lines won't even be that long.

• Good uses of yield. This is something new programmers often skip over and it's good to see it being used to effect here.

• Your imports are clean, minimal and well separated from the rest of the code.

• Some of the naming could use some work to help clarity. Just name your keys enzyme and parasite, rather than e and p. What is a_p? Try not to use built-in function names as argument names (pop) as it can cause issues and be confusing. Here, it is clearly short for population but be careful with it. Use snake_case for naming lower-cased objects ratherthanthis.

• You are frequently returning a huge number of values. If you're always printing 0s to the file you don't need them to be returned, just write them to the file every time, then write the rest of the return values. Some things like gen should be kept track of externally, rather than it being returned every time. If something is static, you probably don't need to feed it into a function then spit it back out unchewed.

• Multi-line strings can be achieved with triple quotes:

example = """
          Like
          This
          """

Back to Numpy

• As I say, to be fast, you need to use Numpy start-to finish in your slow sections. If you generate a list with pure python, then cast it to an array, then put it back to pure python, you often save no time. It can even be slower than just pure python.

• You fitness function for example should instead use element-wise operations.

• If you replace the slowest sections of pure python with pure Numpy, you should see some good improvements. You could try a Code Profiler to find exactly where the hang-ups are.