In [11]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("savefig", dpi=100)
In [12]:
# Simulates 1 trajectory of the model
class Rings:
#Inicialization
def __init__(self, ResInit, MutInit, FRes, FMut, BirthRate, DeathRate):
#System constants
self.FitA = FRes
self.FitB = FMut
#Initial conditions
self.Size = ResInit + MutInit
self.SizeA = ResInit
self.SizeB = MutInit
self.BRate = BirthRate
self.DRate = DeathRate
#Time variables
self.Time = 0
self.Size_t = self.Size
self.SizeA_t = self.SizeA
self.SizeB_t = self.SizeB
self.BRate_t = self.BRate
self.DRate_t = self.DRate
#Birth time variables
self.BProb_t = self.BRate_t / ( self.BRate_t + self.DRate_t )
self.BProbA_t = (self.SizeA_t * self.FitA)/( self.SizeA_t * self.FitA + self.SizeB_t * self.FitB )
self.BProbB_t = (self.SizeB_t * self.FitB)/( self.SizeA_t * self.FitA + self.SizeB_t * self.FitB )
#Death time variables
self.DProb_t = (self.Size_t * self.DRate_t) / ( self.BRate_t + self.Size_t * self.DRate_t )
self.DProbA_t = self.SizeA_t / self.Size_t
self.DProbB_t = self.SizeB_t / self.Size_t
#Graphing
self.TimeEvol = np.array([0])
self.SizeEvol = np.array([self.Size])
self.SizeAEvol = np.array([self.SizeA])
self.SizeBEvol = np.array([self.SizeB])
#fixation probability
self.fB = 0
#simulation variables
AverageSize = 0
FixationProb = 0
#Sets time dependent variables to initial values
def reset(self):
#Time variables
self.Time = 0
self.Size_t = self.Size
self.SizeA_t = self.SizeA
self.SizeB_t = self.SizeB
self.BRate_t = self.BRate
self.DRate_t = self.DRate
#Birth time variables
self.BProb_t = self.BRate_t / ( self.BRate_t + self.DRate_t )
self.BProbA_t = (self.SizeA_t * self.FitA)/( self.SizeA_t * self.FitA + self.SizeB_t * self.FitB )
self.BProbB_t = (self.SizeB_t * self.FitB)/( self.SizeA_t * self.FitA + self.SizeB_t * self.FitB )
#Death time variables
self.DProb_t = (self.Size_t * self.DRate_t) / ( self.BRate_t + self.Size_t * self.DRate_t )
self.DProbA_t = self.SizeA_t / self.Size_t
self.DProbB_t = self.SizeB_t / self.Size_t
#graphing
self.TimeEvol = np.array([0])
self.SizeEvol = np.array([self.Size])
self.SizeAEvol = np.array([self.SizeA])
self.SizeBEvol = np.array([self.SizeB])
#fixation probabilty
self.fB = 0
def update_prob(self):
#Birth time variables
self.BProb_t = self.BRate_t / ( self.BRate_t + self.DRate_t )
if self.SizeA_t > 0:
self.BProbA_t = ((self.SizeA_t - 1) * self.FitA + self.FitB)/( self.SizeA_t * self.FitA + self.SizeB_t * self.FitB )
else:
self.BProbA_t = 0
self.BProbB_t = (self.SizeB_t * self.FitB)/( self.SizeA_t * self.FitA + self.SizeB_t * self.FitB )
#Death time variables
self.DProb_t = (self.Size_t * self.DRate_t) / ( self.BRate_t + self.Size_t * self.DRate_t )
self.DProbA_t = self.SizeA_t / self.Size_t
self.DProbB_t = self.SizeB_t / self.Size_t
#one_run: simulates the process and updates the number of residents and mutants alive at time T
def one_run (self,T):
#Setup initial conditions
self.reset()
while ( self.SizeA_t > 0 and self.SizeB_t > 0 and self.Time < T ):
#Define waiting time until next system update
W = np.random.exponential ( self.BRate + self.Size_t * self.DRate )
#Update time
self.Time = self.Time + W
#Update population size and type
S = np.random.uniform(0,1)
R = np.random.uniform(0,1)
if S < self.BProb_t:
#reproduction
self.Size_t = self.Size_t + 1
if R < self.BProbA_t:
self.SizeA_t = self.SizeA_t + 1
else:
self.SizeB_t = self.SizeB_t + 1
else:
#death
self.Size_t = self.Size_t - 1
if R < self.DProbA_t:
self.SizeA_t = self.SizeA_t - 1
else:
self.SizeB_t = self.SizeB_t - 1
#Update rate? Or is rate a constant?
#Update probabilities
if ( self.Size_t == 0):
print("Simulation: extintion")
break
else:
self.update_prob()
#Save updated time and size
self.TimeEvol = np.append(self.TimeEvol, [self.Time])
self.SizeEvol = np.append(self.SizeEvol, [self.Size_t])
self.SizeAEvol = np.append (self.SizeAEvol, [self.SizeA_t])
self.SizeBEvol = np.append (self.SizeBEvol, [self.SizeB_t])
if (self.SizeA_t == 0):
self.fB = self.fB + 1
def one_run_plot (self):
#Plot
plt.figure(1)
plt.plot(self.TimeEvol, self.SizeEvol, color='brown')
plt.ylabel("Population size")
plt.xlabel("Time")
plt.savefig('One_run_populationsize.eps', format='eps', dpi=1000)
plt.figure(2)
plt.ylabel("Number of mutants")
plt.xlabel("Time")
plt.plot(self.TimeEvol, self.SizeBEvol, color='darkcyan')
plt.savefig('One_run_mutants.eps', format='eps', dpi=1000)
#simulation calls the function one_run(T) repeatedly,
#NumRuns number of times and returns the average number of cells at time T across multiple runs.
def simulation(self, T, NumRuns):
#inicialization
TotalSize = 0
Fa = 0
Fb = 0
counterfa = 0
self.AverageSize = 0
self.FixationProb = 0
TimeEvolAverage = np.array([0])
SizeEvolAverage = np.array([self.Size])
for n in range(0, NumRuns):
self.one_run(T)
counterfa = counterfa + self.fB
TotalSize = TotalSize + self.Size_t
self.AverageSize = TotalSize / NumRuns
self.FixationProb = counterfa / NumRuns
return self.FixationProb
In [14]:
ResInit = 100
MutInit = 1
FRes = 1
FMut = 3
BirthRate = 1
DeathRate = 1
X = Rings(ResInit, MutInit, FRes, FMut, BirthRate, DeathRate)
X.one_run(1000000000)
X.one_run_plot()
In [4]:
FixationProbArray = np.array([0])
PopulationSizeArray = np.array([0])
prob = 0
for n in range(10, 1000, 10):
PopulationSizeArray = np.append(PopulationSizeArray, [n])
prob = X.simulation(100000000000,n)
FixationProbArray = np.append (FixationProbArray, [prob])
In [5]:
plt.plot(PopulationSizeArray, FixationProbArray, color='red')
plt.ylabel("Estimated fixation probability")
plt.xlabel("Number of runs")
plt.savefig('Numberruns_probability.eps', format='eps', dpi=1000)
In [6]:
ResInit = 100
MutInit = 1
FRes = 1
FMut = 3
BirthRate = 1
DeathRate = 1
X = Rings(ResInit, MutInit, FRes, FMut, BirthRate, DeathRate)
X.one_run(1000000000)
X.one_run_plot()
FixationProbArray = np.array([0])
PopulationSizeArray = np.array([0])
prob = 0
for n in range(10, 1000, 50):
ResInit = n
X = Rings(ResInit, MutInit, FRes, FMut, BirthRate, DeathRate)
PopulationSizeArray = np.append(PopulationSizeArray, [n])
prob = X.simulation(100000000000,400)
FixationProbArray = np.append (FixationProbArray, [prob])
In [7]:
print(PopulationSizeArray)
In [8]:
print(FixationProbArray)
In [9]:
plt.plot(PopulationSizeArray, FixationProbArray, color='red')
plt.ylabel("Estimated fixation probability")
plt.xlabel("Population size")
plt.savefig('PopulationSize_probability.eps', format='eps', dpi=1000)