This Jupyter notebook can be used to create input files for the simulation program for all parameter combinations of

1. Mutation type (i-iv)

2. fitness parameters smax (type i-iii), s (type iv) respectively

3. Plasmid copy number n from nRange

Furthermore it can be used to test the simulation programs rescue_maintext, rescue_alternative1, rescue_alternative2

For execution of the program it is recommended to use a cluster

from subprocess import Popen, PIPE
ipyparallel_execution=True # flag for directly executing programm

Set parameters (used to obtain results from Figure 2)

u=0.000001 # mutation probability
N0=100000 # number of initial wild-type individuals
s0=-.1 # fitness parameter of wild-type cells
nMax=20 # maximum plasmid copy number
smaxArray=[1.00,0.50,0.10,0.01] # array of fitness parameters type i-iii
sGDEArray=[0.20,0.10,0.04,0.02] # array of fitness parameters type iv
Nsim=1e4 # number of stochastic simulations

Define helper variables and helper functions

MutationTypes=["Dom","Rec","Int","GDE"] # mutation types i-iv
sArray={
    "Dom": smaxArray,
    "Rec": smaxArray,
    "Int": smaxArray,
    "GDE": sGDEArray
}
nRange=range(1,nMax+1)
N=[N0, 0, 0] # Define array for initital population
mu=1.0 # Death rate (constant)
la0=mu+s0 # Define birth rate of wild-type individuals
 # Define function to obtain birth rate from malthusian fitness
def la(s):
    return mu+s

# Define functions
# to obtain array of birth (la (lambda)) and death rates (mu)
# for cell-types i=0..n for given mutation type _MutType_,
# pcn *n*, fitness parameter *s* (either *s_max* for 
# type i-iii or *s* for type iv)
def laArray(MutType, n, s): # function for birth rate array
    laArray=[la0]
    def la1(i):
        if MutType=="Dom": 
            return la(s)
        elif MutType=="Rec":
            return la0*(i!=n)+la(s)*(i==n)
        elif MutType=="Int":
            return la0+(la(s)-la0)*(i/n)
        elif MutType=="GDE":
            return la0+i*s
#    [print(la1(i)) for i in range(1,n+1)]
    [laArray.append(la1(i)) for i in range(1,n+1)]
    return laArray
def muArray(n): # function for death rate array
    return [mu for i in range(0,n+1)]

# Define function to write input file for scenario with 
# PCN n, array of birth rates la, array of death rates mu
def writeInputFile(filepath, MutType, s, n):
    la=laArray(MutType,n,s)
    mu=muArray(n)
    with open(filepath, 'w') as f:
        f.write("%d #n\n" % n)
        f.write("%f #u\n" % u)
        [f.write("%d " % item ) for item in N]
        f.write("#N \n")
        [f.write("%f " % item ) for item in la]
        f.write("#l \n")
        [f.write("%f " % item ) for item in mu]
        f.write("#m \n")

Write input files for parameter combinations

def inputfilepath(MutType, s, n):
    return "inputfiles/"+MutType+"_s_"+str(s)\
                +"_n_"+str(n)+".inp.dat"
fpArray=[] # List containing all filepaths of inputfiles
for MutType in MutationTypes:
    for s in sArray[MutType]:
        for n in nRange:
            fp=inputfilepath(MutType, s, n)
            fpArray.append(fp)
            writeInputFile(fp, MutType, s, n)

Test simulation with rescue program

# Define function to execute C++ simulation. 
# Program is chosen from variable model 
# (maintext,alternative1,alternative2)
# Parameters obtained from : inputfile.
# Number of simulations : Nsim
def runsimulation(model,inputfile,Nsim):
    process = Popen(["bin/rescue_"+model,inputfile,str(Nsim)], stdout=PIPE)
    (output, err) = process.communicate()
    exit_code = process.wait()
    if exit_code!=0:
        print("simulation error")
    return float(output)
# Define function to save output of simulation to file
def runsimulationSaveToFile(model,inputfile,Nsim):
    output=runsimulation(model,inputfile,Nsim)
    with open(inputfile+'.'+str(model)+'.out', 'w') as f: 
        f.write("%f" % float(output))
    return float(output)

# Test of simulation with parameters: Mutation type "Dominance", s_max=1.0, n=10, model: maintext, 10 Simulations 
# Result: 0.8 
fp=inputfilepath("Dom", 1.0, 10)
NsimTest=10
runsimulation("maintext",fp,NsimTest) 

0.8