Return to Question

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)

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.
So, 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?
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? So 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.)
Of course any advice is highly appreciated!

The code

# -*- 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.

Updates

1. 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.

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)

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

# -*- 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.

In the program the system is a list, the cells are dictionaries with 2 keys: "enzyme" (formerly "e") for the enzymes and "parasite" (formerly "p") for the parasites. The values of the keys are the numbers of the 2 variants.
Our parameters are:
population_size (formerly pop_size): the number of the cells
cell_size: the maximal number of molecules (enzymes+parasites) of cells at which cell division takes place
replication_rate_p (formerly a_p): fitness (replication rate) of the parasites relative to the fitness (replication rate) of the enzymes (for example if replication_rate_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 population_size cells with cell_size // 2 enzymes 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 "enzyme" < 2 are discarded as dead. After that if the number of viable cells is bigger than population_size, we choose population_size of them according to cell fitness ("enzyme"/("enzyme"+"parasite")), and they move on to the next generation. On the other hand, if the number of viable cells is population_size or less, they all move on to the next generation.

I've never studied programming in school. The original program is the result of heavy googling. I'm at a point where I need advice from experienced people. At certain parameter values the program gets quite slow.
So far combinations of these parameter values were examined:
population_size: 100; 200; 500; 1000
cell_size: 20; 50; 100; 200; 500; 1000
replication_rate_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 population_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 replication_rate_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 population_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.

So, 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?
In which directions should I improve my programming skills in Python to efficiently implement these kind of models?
Or am I near the full potential of Python in this regard? So 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.)
Of course any advice is highly appreciated!

The code

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

Updated code with stylistic changes based on HoboProber's answer. The original code is below the updated one.

# -*- coding: utf-8 -*-
"""
Simulation of an enzyme-parasite system,
explicitly assuming compartmentalization.
"""
from random import choices, random
import csv
import time
import numpy as np


def simulation(population_size, cell_size, replication_rate_p, mutation_rate, gen_max):
    """Function doing the simulation."""
    def fitness(population):
        """Function for calculating fitnesses of cells."""
        return [i["enzyme"] / (i["enzyme"] + i["parasite"]) for i in population]

    def population_stats(population):
        """Function for calculatig statistics of the system."""
        gyaklist_e = [i["enzyme"] for i in population]
        gyaklist_p = [i["parasite"] for i in population]
        fitness_list = fitness(population)
        return (
            sum(gyaklist_e), sum(gyaklist_p),
            sum([1 for i in population if i["enzyme"] > 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(fitness_list), np.var(fitness_list),
            np.percentile(fitness_list, 25),
            np.percentile(fitness_list, 50),
            np.percentile(fitness_list, 75)
            )

    # Creating the system with the starting state being
    # half full cells containing only enzymes.
    population = [{"enzyme": cell_size // 2,
                   "parasite": 0} for _ in range(population_size)]
    gen = 0
    yield (gen, *population_stats(population), population_size,
           cell_size, mutation_rate, replication_rate_p, "aft")
    print(f"N = {population_size}, rMax = {cell_size}, "
          f"aP = {replication_rate_p}, U = {mutation_rate}")

    while population and gen < gen_max:
        gen += 1

        # Replicator proliferation until cell_size in each cell.
        for i in population:
            while not i["enzyme"] + i["parasite"] == cell_size:
                lucky_replicator = choices(
                    ["enzyme", "parasite"],
                    [i["enzyme"], replication_rate_p*i["parasite"]]
                    )
                if lucky_replicator[0] == "enzyme" and random() < mutation_rate:
                    lucky_replicator[0] = "parasite"
                i[lucky_replicator[0]] += 1

        if gen % 100 == 0:
            yield (gen, *population_stats(population), population_size,
                   cell_size, mutation_rate, replication_rate_p, "bef")

        # Each cell divides.
        new_population = [{"enzyme": np.random.binomial(i["enzyme"], 0.5),
                           "parasite": np.random.binomial(i["parasite"], 0.5)}
                          for i in population]
        for i in zip(population, new_population):
            i[0]["enzyme"] -= i[1]["enzyme"]
            i[0]["parasite"] -= i[1]["parasite"]
        population += new_population

        # Discarding dead cells.
        new_population = [i for i in population if i["enzyme"] > 1]

        # Choosing surviving cells according to their fitnesses
        # if there are viable cells.
        # If there are no viable cells, the simulation is over.
        if new_population:
            fitness_list = fitness(new_population)
            fitness_sum = np.sum(fitness_list)
            fitness_list = fitness_list / fitness_sum
            population = np.random.choice(
                new_population, min(population_size, len(new_population)),
                replace=False, p=fitness_list
                ).tolist()
        else:
            population = new_population
            for i in range(2):
                yield (gen+i, *(0, 0)*9, population_size,
                       cell_size, mutation_rate, replication_rate_p, "aft")
            print(f"{gen} generations are done. Cells are extinct.")

        if gen % 100 == 0 and population:
            yield (gen, *population_stats(population), population_size,
                   cell_size, mutation_rate, replication_rate_p, "aft")

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


def write_out_file(result, n_run):
    """Function writing csv output file."""
    local_time = time.strftime("%m_%d_%H_%M_%S_%Y", time.localtime(time.time()))
    with open("output_data_" + local_time + ".csv", "w", newline="") as out_file:
        out_file.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"
            )
        out_file = csv.writer(out_file, delimiter=";")
        counter = 0
        print(counter, "/", n_run)
        for i in result:
            out_file.writerows(i)
            counter += 1
            print(counter, "/", n_run)


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

The original code:

# -*- 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.

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.

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.
So, 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?
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? So 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.)
Of course any advice is highly appreciated!

The code

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.

Updates

1. 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.

# -*- coding: utf-8 -*-
"""
Simulation of an enzyme-parasite system,
explicitly assuming compartmentalization.
"""
from random import choices, random
import csv
import time
import numpy as np


def simulation(population_size, cell_size, replication_rate_p, mutation_rate, gen_max):
    """Function doing the simulation."""
    def fitness(population):
        """Function for calculating fitnesses of cells."""
        return [i["enzyme"] / (i["enzyme"] + i["parasite"]) for i in population]

    def population_stats(population):
        """Function for calculatig statistics of the system."""
        gyaklist_e = [i["enzyme"] for i in population]
        gyaklist_p = [i["parasite"] for i in population]
        fitness_list = fitness(population)
        return (
            sum(gyaklist_e), sum(gyaklist_p),
            sum([1 for i in population if i["enzyme"] > 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(fitness_list), np.var(fitness_list),
            np.percentile(fitness_list, 25),
            np.percentile(fitness_list, 50),
            np.percentile(fitness_list, 75)
            )

    # Creating the system with the starting state being
    # half full cells containing only enzymes.
    population = [{"enzyme": cell_size // 2,
                   "parasite": 0} for _ in range(population_size)]
    gen = 0
    yield (gen, *population_stats(population), population_size,
           cell_size, mutation_rate, replication_rate_p, "aft")
    print(f"N = {population_size}, rMax = {cell_size}, "
          f"aP = {replication_rate_p}, U = {mutation_rate}")

    while population and gen < gen_max:
        gen += 1

        # Replicator proliferation until cell_size in each cell.
        for i in population:
            while not i["enzyme"] + i["parasite"] == cell_size:
                lucky_replicator = choices(
                    ["enzyme", "parasite"],
                    [i["enzyme"], replication_rate_p*i["parasite"]]
                    )
                if lucky_replicator[0] == "enzyme" and random() < mutation_rate:
                    lucky_replicator[0] = "parasite"
                i[lucky_replicator[0]] += 1

        if gen % 100 == 0:
            yield (gen, *population_stats(population), population_size,
                   cell_size, mutation_rate, replication_rate_p, "bef")

        # Each cell divides.
        new_population = [{"enzyme": np.random.binomial(i["enzyme"], 0.5),
                           "parasite": np.random.binomial(i["parasite"], 0.5)}
                          for i in population]
        for i in zip(population, new_population):
            i[0]["enzyme"] -= i[1]["enzyme"]
            i[0]["parasite"] -= i[1]["parasite"]
        population += new_population

        # Discarding dead cells.
        new_population = [i for i in population if i["enzyme"] > 1]

        # Choosing surviving cells according to their fitnesses
        # if there are viable cells.
        # If there are no viable cells, the simulation is over.
        if new_population:
            fitness_list = fitness(new_population)
            fitness_sum = np.sum(fitness_list)
            fitness_list = fitness_list / fitness_sum
            population = np.random.choice(
                new_population, min(population_size, len(new_population)),
                replace=False, p=fitness_list
                ).tolist()
        else:
            population = new_population
            for i in range(2):
                yield (gen+i, *(0, 0)*9, population_size,
                       cell_size, mutation_rate, replication_rate_p, "aft")
            print(f"{gen} generations are done. Cells are extinct.")

        if gen % 100 == 0 and population:
            yield (gen, *population_stats(population), population_size,
                   cell_size, mutation_rate, replication_rate_p, "aft")

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


def write_out_file(result, n_run):
    """Function writing csv output file."""
    local_time = time.strftime("%m_%d_%H_%M_%S_%Y", time.localtime(time.time()))
    with open("output_data_" + local_time + ".csv", "w", newline="") as out_file:
        out_file.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"
            )
        out_file = csv.writer(out_file, delimiter=";")
        counter = 0
        print(counter, "/", n_run)
        for i in result:
            out_file.writerows(i)
            counter += 1
            print(counter, "/", n_run)


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

# -*- coding: utf-8 -*-
"""
Simulation of an enzyme-parasite system,
explicitly assuming compartmentalization.
"""
from random import choices, random
import csv
import time
import numpy as np


def simulation(population_size, cell_size, replication_rate_p, mutation_rate, gen_max):
    """Function doing the simulation."""
    def fitness(population):
        """Function for calculating fitnesses of cells."""
        return [i["enzyme"] / (i["enzyme"] + i["parasite"]) for i in population]

    def population_stats(population):
        """Function for calculatig statistics of the system."""
        gyaklist_e = [i["enzyme"] for i in population]
        gyaklist_p = [i["parasite"] for i in population]
        fitness_list = fitness(population)
        return (
            sum(gyaklist_e), sum(gyaklist_p),
            sum([1 for i in population if i["enzyme"] > 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(fitness_list), np.var(fitness_list),
            np.percentile(fitness_list, 25),
            np.percentile(fitness_list, 50),
            np.percentile(fitness_list, 75)
            )

    # Creating the system with the starting state being
    # half full cells containing only enzymes.
    population = [{"enzyme": cell_size // 2,
                   "parasite": 0} for _ in range(population_size)]
    gen = 0
    yield (gen, *population_stats(population), population_size,
           cell_size, mutation_rate, replication_rate_p, "aft")
    print(f"N = {population_size}, rMax = {cell_size}, "
          f"aP = {replication_rate_p}, U = {mutation_rate}")

    while population and gen < gen_max:
        gen += 1

        # Replicator proliferation until cell_size in each cell.
        for i in population:
            while not i["enzyme"] + i["parasite"] == cell_size:
                lucky_replicator = choices(
                    ["enzyme", "parasite"],
                    [i["enzyme"], replication_rate_p*i["parasite"]]
                    )
                if lucky_replicator[0] == "enzyme" and random() < mutation_rate:
                    lucky_replicator[0] = "parasite"
                i[lucky_replicator[0]] += 1

        if gen % 100 == 0:
            yield (gen, *population_stats(population), population_size,
                   cell_size, mutation_rate, replication_rate_p, "bef")

        # Each cell divides.
        new_population = [{"enzyme": np.random.binomial(i["enzyme"], 0.5),
                           "parasite": np.random.binomial(i["parasite"], 0.5)}
                          for i in population]
        for i in zip(population, new_population):
            i[0]["enzyme"] -= i[1]["enzyme"]
            i[0]["parasite"] -= i[1]["parasite"]
        population += new_population

        # Discarding dead cells.
        new_population = [i for i in population if i["enzyme"] > 1]

        # Choosing surviving cells according to their fitnesses
        # if there are viable cells.
        # If there are no viable cells, the simulation is over.
        if new_population:
            fitness_list = fitness(new_population)
            fitness_sum = np.sum(fitness_list)
            fitness_list = fitness_list / fitness_sum
            population = np.random.choice(
                new_population, min(population_size, len(new_population)),
                replace=False, p=fitness_list
                ).tolist()
        else:
            population = new_population
            for i in range(2):
                yield (gen+i, *(0, 0)*9, population_size,
                       cell_size, mutation_rate, replication_rate_p, "aft")
            print(f"{gen} generations are done. Cells are extinct.")

        if gen % 100 == 0 and population:
            yield (gen, *population_stats(population), population_size,
                   cell_size, mutation_rate, replication_rate_p, "aft")

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


def write_out_file(result, n_run):
    """Function writing csv output file."""
    local_time = time.strftime("%m_%d_%H_%M_%S_%Y", time.localtime(time.time()))
    with open("output_data_" + local_time + ".csv", "w", newline="") as out_file:
        out_file.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"
            )
        out_file = csv.writer(out_file, delimiter=";")
        counter = 0
        print(counter, "/", n_run)
        for i in result:
            out_file.writerows(i)
            counter += 1
            print(counter, "/", n_run)


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