细菌种群模拟

模拟疾病和细菌种群的传播

设计并实现细菌种群的动态随机模拟，并根据模拟的结果得出各种治疗方案如何影响细菌传播的结论。

背景

细菌是单细胞生物，是无性繁殖的，有些细菌会导致疾病，有些是无害的，还有些是有益的。坏细菌可引起如链球菌感染等感染。细菌感染用抗生素治疗，抗生素能杀死坏的细菌细胞。
引起感染的细菌可以通过两种方式抵抗和发展抗生素的恢复能力。

• 自然情况；
• 通过不正当使用抗生素。因此，细菌群体在患者体内可进化。在这个问题集中，我们将利用模拟来探索将抗生素引入细菌群体的效果，并确定如何在模型中最好的解决这些治疗挑战。我们将实施一个高度简化的细菌种群动态模型，该模型具有生物学相关的特征。

Problem 1: Implementing a Simple Simulation

这里我们先考虑一个简单的模型：患者不服用任何抗生素，我们只是简单地模拟细菌种群的增长过程，相当于患者没有得到治疗的情况。
实现SimpleBacteria 类和 Patient 类的功能设计。

SimpleBacteria 实现细节：

init() ：初始化细菌细胞对象
这里重点是两个参数的理解，虽然按代码注释翻译为繁殖率和死亡率比较直接，但是这么理解对理解整个题目表达的意思没有多大的帮助，因此，这里理解为“可能性”一词较为合适。

• birth_prob:细胞发生繁殖的可能性
• death_prob：细胞已经死亡的可能性

# 初始化
def __init__(self, birth_prob, death_prob):self.birth_prob = birth_probself.death_prob = death_prob

is_killed(): 随机地确定该细菌细胞在患者体内是否会在一个时间步骤内被杀死。这里重点是要理解死亡率和随机确定的概率之间的关系，明确判断细菌细胞死亡的判定条件。首先随机产生一个概率随机数 random.random() 这个数作为判断细菌是否死亡的标准。我们可以这么理解，这个数是一个死亡系数，那么按常理，当细菌细胞自身的最大死亡率self.death_prob（理解为细胞已经死亡的可能性）大于这死亡系数时，就可以认为这个细菌细胞是已经死亡了的。因此，代码如下：

# 判断细胞是否是被杀死了
def is_killed(self):return self.death_prob > random.random()

timestep: 意思容易理解，翻译起来却总觉得不合适，这个字面意思是时间步骤，但是思前想后，觉得把它翻译为 周期 比较合适，一个时间段 timestep 就是一个周期。

reproduce():随机确定该细菌细胞是否在一个周期出现繁殖。这个方法会被Patient类和TreatedPatient类中的update（）方法调用。

• pop_density: 种群密度，定义为当前细菌种群数量除以最大种群数量，如果这种细菌细胞繁殖，那么reproduce（）会创建并返回细菌细胞的后代即SimpleBacteria的实例（它具有与其父亲相同的birth_prob和death_prob值）。
• SimpleBacteria：一个新的实例，代表这种细菌细胞的后代（如果细菌繁殖）。孩子应该具有与该细菌相同的birth_prob和death_prob值。

基于上面提到过的理解思想，这里重点依旧是判断细菌细胞是否出现繁殖的条件。首先，代码给出了繁殖可能性reproduce的计算公式,那么当这个可能性大于得到的随机概率时，我们就认为细胞出现了繁殖。代码如下：

# 细胞是否出现繁殖，若繁殖，返回繁殖的对象
def reproduce(self, pop_density):if (self.birth_prob * (1 - pop_density)) >= random.random():return SimpleBacteria(self.birth_prob, self.death_prob)# 不繁殖 抛出异常else:raise NoChildException('Not Reproducible')

Patient实现细节：
init():初始化患者个体

def __init__(self, bacteria, max_pop):self.bacteria = bacteriaself.max_pop = max_pop

get_total_pop():获取患者体内细菌细胞的数量，返回列表长度即可。

# 获取患者体内细菌细胞的数量
def get_total_pop(self):return len(self.bacteria)

update():更新患者体内细菌种群状态的变化。首先得要理解细菌种群变化的是什么？从返回结果也容易知道，其实就是细菌细胞数量的变化，引起细菌细胞数量变化的原因主要有两点：一是细胞死亡，需要减去细胞死亡的数量；而是细胞繁殖，需要增加新生细菌细胞的数量。代码如下：

# 更新患者体内细菌种群状态的变化
def update(self):new_bacteria = []for ind, bac in enumerate(self.bacteria):# 移除已经死亡的细菌细胞if bac.is_killed():self.bacteria.pop(ind)else:pop_density = self.get_total_pop() / float(self.max_pop)try:# 增加新繁殖的后代细菌细胞new_bacteria.append(bac.reproduce(pop_density))except NoChildException:continueself.bacteria = self.bacteria + new_bacteria# 返回更新后患者体内的细菌细胞数量return self.get_total_pop()

Problem 2: Running and Analyzing a Simple Simulation

模拟引入抗生素之前的细菌细胞

calc_pop_avg(): 获取n个周期实验的平均的细菌种群大小，for循环累加求和再求平均数即可。

# 获取经过n个周期实验过程的平均的细菌种群大小
def calc_pop_avg(populations, n):sum_pop = 0for i in range(len(populations)):sum_pop += populations[i][n]return sum_pop / len(populations)

simulation_without_antibiotic(num_bacteria,max_pop,birth_prob,death_prob, num_trials):
各个参数的理解：

• num_bacteria表示患者体内细菌实例的数量，即SampleBacteria对象的数量
• max_pop:患者体内最大的细菌种群大小birth_prob:出现繁殖的最大可能性
• death_prob:细胞已经死亡的最大可能性 num_trails表示实验次数。
• 返回的结果是一个二维列表populations:其中每一项populations[ i ][ j ]的含义，i表示的是实验次数第i次，j表示的是第j个bacteria对象。
  所有的实验时间取 300个周期

这里理解起来比较绕，也花费了不少时间捋其中的关系。我是这么求解的，首先得明确结果列表里每一项是什么，显然列表长度是300，列表的每个元素是一个状态，上面在首先类Patient时有点提示，那就是会用到 update()函数。因此这里的结果列表的每一个元素就是patient.update()，要求patient.update(),那就需要求patient,即需要求patient的两个参数,max_pop已经给出来了，剩下的就是bacterias 这个for循环累加求和即可实现。
代码如下：

def simulation_without_antibiotic(num_bacteria,max_pop,birth_prob,death_prob,num_trials):populations = []for i in range(0, num_trials):bacterias = [SimpleBacteria(birth_prob, death_prob) for j in range(0,num_bacteria)]patient = Patient(bacterias, max_pop)populations.append([patient.update() for i in range(0, 300)])return populations

测试代码

populations = simulation_without_antibiotic(100, 1000, 0.1, 0.025, 50)
# Plottingave_list = []
for i in range(300):ave_list.append(calc_pop_avg(populations, i))
make_one_curve_plot(range(300), ave_list, 'Timestep', 'Average Population', 'Without Antibiotic')

测试结果

Problem 3: Calculating a Confidence Interval

计算 populations 在时间步 t 上的标准差，并且不能使用计算标准差的库函数或者建立在库函数方法基础上的计算方法。

各个参数的含义：
populations：二维列表，populations[ i ][ j ]表示第 i 次实验中 第 j 个时间步上的 bacteria 数量
t: 一个整数，代表时间步长。
不能采用相关库函数，那只能按照公式求出一个个的变量来了，problem3主要就是按公式计算正确的结果并返回。
均值计算公式如下：

标准差计算公式如下：

平均值的标准误差：

代码如下：

# 求 populations 在 t 上的标准差
def calc_pop_std(populations, t):# 先求均值mean = calc_pop_avg(populations, t) std = 0for population in range(len(populations)):std += (populations[population][t] - mean) ** 2# 再求标准差return math.sqrt(std * (1 / len(populations)))

def calc_95_ci(populations, t):mean = calc_pop_avg(populations, t)  # 均值sd = calc_pop_std(populations, t)  # 标准差return (mean, 1.96 * sd / math.sqrt(len(populations)))

执行测试文件，结果ok

Problem 4: Implementing a Simulation with an Antibiotic

在这个问题中，我们考虑了两方面的效果

• 一是给患者施用抗生素的效果
• 二是细胞遗传或突变导致对抗生素产生的抵抗性。

随着细菌群体的繁殖，细菌后代将发生突变。一些细菌细胞获得良性突变，如抗生素耐药性。注意，种群密度越低的细菌发生变异的速度越快。

ResistantBacteria 类的实现，继承自类SimpleBacteria，表示具有抗生素抵抗的细菌类。
涉及到的新参数含义：

• resistant： 一个 bool 值，表示是否对抗生素有抵抗特性
  剩下的参数跟前几个 problem 中同名的参数含义是一样的
  里面包含的方法基本和 父类 SimpleBacteria 中一样，但是计算方法有区别，因此都需要重写。大体思路还是跟实现父类一致的，这里主要给出 细菌繁殖函数 reproduce()的代码设计和解释。
  首先。细胞只有繁殖了，才有可能突变，才有可能对抗生素产生抵抗性，但是繁殖了不是一定会发生突变。因此，需要用到嵌套的条件语句。
  伪代码如下：

if 细胞繁殖：if 细胞自身就已经有抗生素抵抗性 或者 后代细胞突变成有抵抗性的细胞：产生有抵抗性的后代细胞，即返回新对象else:产生无抵抗性的后代细胞
else:未繁殖，不产生后代细胞，返回 False

代码如下：

    # 细菌繁殖def reproduce(self, pop_density):reproduce = self.birth_prob * (1 - pop_density)  # 繁殖率mutate = self.mut_prob * (1 - pop_density)   # 突变率# 细胞繁殖if (reproduce > random.random()):# 产生对抗生素有抵抗性的后代if self.resistant or random.random() > mutate:return ResistantBacteria(self.birth_prob, self.death_prob, self.resistant, self.mut_prob)# 产生对抗生素无抵抗的后代return ResistantBacteria(self.birth_prob, self.death_prob, not self.resistant, self.mut_prob)# 细胞未发生繁殖else:raise NoChildException('Not Reproducible')

TreatedPatient​ 类的设计：
TreatedPatient​ 表示接受治疗的患者，继承自父类 Patient
类中各个方法实现细节：

    def __init__(self, bacteria, max_pop):Patient.__init__(self, bacteria, max_pop)self.on_antibiotic = False# 患者服用抗生素def set_on_antibiotic(self):self.on_antibiotic = True# 获取具有抗生素抵抗性细胞数量def get_resist_pop(self):count = 0for x in self.bacteria:if x.get_resistant():count += 1return count

下面重点讲一下 细菌状态更新的 函数 update() 功能设计.
整体思路跟 Patient 类中的 update 设计思路是一样的 ，减去已经杀死的细胞，增加繁殖后的细胞，但是情况又有点区别。具体看代码注释。

    def update(self):new_bacteria = []for ind, bac in enumerate(self.bacteria):# 减去死亡的细胞if bac.is_killed():self.bacteria.pop(ind)else:# 减去不具有抗生素抵抗性的细胞if self.get_antibiotics() and not bac.get_resistant():self.bacteria.pop(ind)else:pop_density = self.get_total_pop() / float(self.max_pop)try:# 增加繁殖的后代new_bacteria.append(bac.reproduce(pop_density))except NoChildException:continueself.bacteria = self.bacteria + new_bacteria# 返回一个周期后的细胞总数量return self.get_total_pop()

Problem 5: Running and Analyzing a Simulation with an Antibiotic

模拟服用抗生素的患者体内细菌细胞的变化，分析结果。实验时常设置为400个周期，前150个周期不给患者服用抗生素，后250个周期给患者服用抗生素。
跟problem 2的做法类似，列表存储for循环每次实验的结果。

def simulation_with_antibiotic(num_bacteria,max_pop,birth_prob,death_prob,resistant,mut_prob,num_trials):num_step = 400  # 实验时常 400个周期populations = []resistant_pop = []for tri in range(num_trials):# 初始化populations.append([None] * num_step)resistant_pop.append([None] * num_step)bac_list = []for bac in range(num_bacteria):bac_list.append(ResistantBacteria(birth_prob, death_prob, resistant, mut_prob))# 创建一个可接受治疗的患者对象pa = TreatedPatient(bac_list, max_pop)for step in range(num_step):# 150个周期后服用抗生素if step >= 150:pa.set_on_antibiotic()populations[tri][step] = pa.update()resistant_pop[tri][step] = pa.get_resist_pop()# 不服用抗生素else:populations[tri][step] = pa.update()resistant_pop[tri][step] = pa.get_resist_pop()return (populations, resistant_pop)

各个参数的含义：

• num_bacteria: 细菌细胞数量

• max_pop：细菌种群细菌数量最大时的细菌数量

• birth_prob：出生率

• death_prob：死亡率

• resistant：bool值，细菌细胞是否具有抗生素抵抗性

• mut_prob：突变率

• num_trials：实验次数
  测试代码

ave_list1 = []
ave_list2 = []
num_step = 400
for i in range(num_step):ave_list1.append(calc_pop_avg(total_pop, i))
for i in range(num_step):ave_list2.append(calc_pop_avg(resistant_pop, i))# make_two_curve_plot(range(num_step), ave_list1, ave_list2, 'Total', 'Resistant',
#                     'Average Population', 'Without Antibiotic', 'Simulation of A')make_two_curve_plot(range(num_step), ave_list1, ave_list2, 'Total', 'Resistant','Average Population', 'Without Antibiotic', 'Simulation of B')

测试结果

可知，两条曲线最终均收敛到0

Problem 6 Write-up

回答下列问题：

1. What happens to the total population before introducing the antibiotic?

	Sim A: total population 一开始增长，后逐渐收敛Sim B: population 一开始增长，后逐渐收敛

2. What happens to the resistant bacteria population before introducing the antibiotic?

    Sim A: 先增长后下降Sim B: 先增长后下降

3. What happens to the total population after introducing the antibiotic?

    Sim A: 突然下降然后轻微增加，最终收敛Sim B: 突然下降并减少到零

4. What happens to the resistant bacteria population after introducing the antibiotic?

    Sim A: 轻微增加然后收敛Sim B: 急速下降到 0