# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import random
import math
import numpy as np
from numpy import searchsorted as mychoose
from scipy.signal import fftconvolve
plt.rcParams.update({'font.size': 13})
plt.rcParams.update({'axes.labelsize': 'large'})
##########################################################
'''a stochastic simulation of many infected cells at once
for the two types of spores model'''
#the variables are
#s = number of intracellular spores
#n = number of intracellular newly germinated bacteria
#b = number of intracellular vegetative bacteria
#c = 1 if the cell is alive, and c = 0 otherwise
#define the events that can occur inside an infected cell
def germination(s,n,b,c):
return (s-1,n+1,b,c)
def maturation(s,n,b,c):
return (s,n-1,b+1,c)
def death_ngb(s,n,b,c):
return (s,n-1,b,c)
def birth(s,n,b,c):
return (s,n,b+1,c)
def death(s,n,b,c):
return (s,n,b-1,c)
def rupture(s,n,b,c):
return (0,0,0,0)
#time-course for the simulation
sim_times = np.array([0,2,4,6,8,10,15,20,25,30,35,40,45,50])
#function to run a simulation for an intial number of infected cells/ intracellular spores given by 'spores'
#outputs arrays for the total number of spores, newly germinated bacteria, and vegetative bacteria inside cells at each of the time points in 'sim_times'
#arising from each of the two types of spores.
#also outputs the list of times when a cell ruptures
def simulation(e,gA,gB,mut,lamb,mu,gamma,spores):
#initial number of cells, each containing one spore
x0 = spores
#initialise time-courses for number of spores, ngb, and vegetative bacteria
#for each of the two types of spores
spAstar = np.zeros(len(sim_times))
spBstar = np.zeros(len(sim_times))
ngbAstar = np.zeros(len(sim_times))
ngbBstar = np.zeros(len(sim_times))
bacAstar = np.zeros(len(sim_times))
bacBstar = np.zeros(len(sim_times))
#initialise a list for the rupture times
rupture_times = [0]
events = {0:germination, 1:maturation, 2:death_ngb, 3:birth, 4:death, 5:rupture}
#simulate each of the initial x0 cells
for _ in range(x0):
#germination rate is gA with probability e and gB otherwise
g = random.choices(population=[gA,gB],weights=[e,1-e],k=1)[0]
#set initial time to 0
t = 0
#initial number of intracellular spores
s = 1
#initial number of intracellular newly germinated bacteria
n = 0
#initial number of intracellular vegetative bacteria
b = 0
#cell is initially alive
c = 1
#initialise lists for the number in each population
#and the transition times
S = []
N = []
B = []
T = []
S.append(s)
N.append(n)
B.append(b)
T.append(t)
#maximum number of hours to run the simulation for
tmax = max(sim_times)
#run the simulation until the cell has either ruptured, recovered, or the time is over tmax
while (s+n+b)!=0 and t