数学生物学-1-指数增长与衰退(Exponential Growth and Decay)

我们注意到数学专业的学生往往没有生物学系统的思维，而生物专业的学生则常常对数学感到恐惧。许多生物专业的学生可能一直对科学感兴趣，部分原因是他们认为生物学是一个不需要任何数学技能或背景就可以研究的科学领域。我认为这是不对的思维。

虽然大多数生物学本科课程中可能包含统计学，但通常不包括微积分及以上的数学内容。我们认为，生物学教育中普遍存在的情况至少在一定程度上忽视了数学。因此，学生在研究生阶段以及学术或其他研究工作中学到的数学知识很少。

因此，本专栏的目标是以一种不会让生物学本科生过度恐惧的方式，介绍这个跨学科领域，并强调通过应用基于数学的方法可以了解生物学系统的内容。

本专栏参考An Invitation to Mathematical Biology

1.简介

2.指数增长与衰退

(1)指数增长

[1]例一：理想指数增长模型

[2]例2：承载能力模拟模型

（2）指数衰退

[1]例1：理想指数衰减模型

[2]感染康复模型

1.简介

用数学去解释生物问题往往可以从宏观的和整体角度解决，可以有解析解与数值解解法，在计算生物学中，我们通常开发新算法然后应用到计算机软件中，对一些问题的数值解进行计算，快速浏览几本数学生物学期刊，你会发现研究论文通常以计算方法为主。

数学与生物学之间加强联系的重要性在NSF赞助的题为“数学与生物学：接口、挑战与机遇”的研讨会上得到了清晰的阐述，作者们描述了数学有可能像牛顿微积分对物理学那样彻底改变生物学（Levin 1992）。这种互动将使双方受益：正如过去动力学系统和混沌理论在很大程度上是为了解决生物问题而发展起来的一样，这种互动也将丰富数学。

2.指数增长与衰退

许多物理现象遵循（近似）指数增长或指数衰减的规律，显然生物现象也不例外，正如我们接下来将要看到的。这种增长和衰减情况的典型例子相当精确，如细菌生长、传染病康复和放射性衰减

(1)指数增长

[1]例一：理想指数增长模型

在合适营养环境中，细菌在一定的时间间隔Δt内体积会增加。假设我们在培养基中放置了60个大肠杆菌细胞，每个细胞每20分钟分裂成两个。那么在第一个20分钟后，就有120个细胞，再过20分钟就有240个细胞，然后在60分钟后有480个细胞，依此类推，只要培养基中有食物供应，细胞数量就会持续翻倍。如果在180分钟内（即9个20分钟的时间段），培养基中的食物充足，那么最终会有60×2^9=30,720个细胞。

如果说时间趋向∞，那么细胞数目也就会趋向∞，这时我们考虑到培养基营养不是无限的及细胞寿命不是无限的，因此会有衰减的情况，这样想是对的，但是我更偏向于模型是一步一步完善的，我们先考虑最简单的模型。

这种增长模式可以用离散时间方程来表示，即：

$N_{t+1} = N_{t}+kN_{t}$

为了计算每分钟的增长率，我们需要将20分钟内的翻倍效应转换为每分钟的增长率：

$k = 2^{\frac{1}{20}}-1$

具体来说，2^{1/20} 表示每分钟的增长因子，它告诉我们每分钟细胞数量相对于前一分钟增加了多少。然后，我们从这个增长因子中减去1，得到的实际增长率 k 就是每分钟细胞数量相对于前一分钟的净增长比例。

我们通过python可视化看一下：

import matplotlib.pyplot as plt

# 定义初始细胞数量和时间间隔（分钟）
N0 = 60
interval = 20

# 定义时间点（分钟）
time_points = [i*interval for i in range(10)]  # 从0到180分钟，每20分钟一个点

# 计算每个时间点的细胞数量（指数增长）
cell_counts = [N0 * 2**(i/interval) for i in time_points]

# 创建图表
plt.figure(figsize=(10, 6))
plt.plot(time_points, cell_counts, marker='o', linestyle='-')

# 添加标题和标签
plt.title('Escherichia coli Cell Growth Over Time (Exponential)')
plt.xlabel('Time (min)')
plt.ylabel('Number of Cells')

# 显示网格
plt.grid(True)

# 显示图表
plt.show()

这个例子中，如果时间低于一个小时，那还是比较可信的，但是如果用于估计长时间的问题就变得离谱了。

[2]例2：承载能力模拟模型

前面的指数增长例子似乎表明细胞数量会无限增加。然而，我们可以想象，对于许多实际的生物种群，存在许多限制因素阻止了无限制的增长。在这里，我们将考虑在实验室条件下生长的细菌种群。数据由内华达大学拉斯维加斯分校化学与生物化学系的Ernesto Abel-Santos提供。实验测量了某种培养基中炭疽杆菌在不同小时数的生长情况。从表中可以看出，最初有一个快速增长（如上面例子所示），但在大约10小时后，增长放缓，最终趋于平稳，接近所谓的培养基的“承载能力”。这种模式有时被称为S曲线或逻辑曲线。

我们可以用一条遵循逻辑函数的曲线来拟合这组数据。在这种情况下，我们使用一个函数，该函数最小化了数据点和拟合曲线之间的误差，这种方法被称为非线性最小二乘法。拟合曲线似乎合理地近似了实际观测到的数据。曲线拟合过程还为我们提供了两个重要参数的估计：一个增长率和一个承载能力。我们通过R语言进行拟合,分别使用二次多项式拟合和非线性最小二乘法拟合：

# 准备数据
time <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15)
population <- c(0.0031, 0.0046, 0.0094, 0.0218, 0.0469, 0.0759, 0.0989, 0.1172, 0.1310, 0.1407, 0.1481, 0.1532, 0.1551, 0.1575, 0.1580, 0.1590)

# 使用多项式回归进行拟合，这里以二次多项式为例
fit <- lm(population ~ poly(time, 2), data = data.frame(time, population))

# 打印模型的系数
coefficients(fit)

# 绘制原始数据和拟合曲线
plot(time, population, type = "b", col = "blue", main = "Population Size Over Time")
lines(time, fitted(fit), col = "red", lwd = 2)
legend("topright", legend = c("Observed", "Fitted"), col = c("blue", "red"), lty = 1)

# 使用nls()函数进行非线性最小二乘拟合
# 假设初始参数K=0.16, r=0.5, t0=5
model <- nls(population ~ K / (1 + exp(-r * (time - t0))), 
             start = list(K = 0.16, r = 0.5, t0 = 5))

# 查看拟合结果
summary(model)

# 获取拟合参数
params <- coef(model)
K_fit <- params['K']
r_fit <- params['r']
t0_fit <- params['t0']

# 绘制原始数据和拟合曲线
plot(time, population, type = "b", col = "blue", main = "Population Size Over Time")
curve(K_fit / (1 + exp(-r_fit * (x - t0_fit))), add = TRUE, col = "red", lwd = 2)
legend("topleft", legend = c("Observed", "Fitted"), col = c("blue", "red"), lty = 1)

这是二次多项式拟合的结果，其实还能接受，但是不准确，我们看看非线性最小二乘拟合：

这个是已知结果去求方程，也就是数值拟合！这对于我们理解表观数据很有用，另外，我们可以先定义模型，然后用模型求解，看是否符合数据，这是数值模拟的思维，我们用python实现：

我们关注到这是一个logistics模型：

# 导入所需的库
import numpy as np  # 导入NumPy库，用于数值计算
import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘图

# 定义模型参数
r = 4.41708  # 增长率
K = 0.15618  # 环境承载量
steps = 15  # 模拟的步数

# 初始化种群比例和环境承载量比例的数组
p = np.zeros(steps+1)  # 种群比例数组，多出一个元素用于存储最终结果
n = np.zeros(steps+1)  # 环境承载量比例数组，同上
p[0] = 0.0031  # 初始种群比例
n[0] = 0.0  # 初始环境承载量比例

# 进行模拟
for i in range(0, steps):
    p[i+1] = p[i] + r*p[i]*(K- p[i])  # 根据离散逻辑斯蒂模型更新种群比例
    n[i+1] = n[i] + 1.0  # 更新时间步

# 打印模拟数据
print(np.transpose([n,p]))  # 打印时间和种群比例的转置矩阵

# 绘制结果图表
plt.plot(n,p, 'b.', fillstyle='full')  # 绘制时间与种群比例的图表，点以蓝色填充
plt.xlabel('Time')  # 设置横轴标签为“Time”（时间）
plt.ylabel('p')  # 设置纵轴标签为“p”（种群比例）
plt.show()  # 显示图表

回顾之前的表格和图表，我们可以注意到承载能力的值，通过曲线拟合估计为0.1562。事实上，很自然地推断出变化Δpn= pn+1- pn与乘积pn (0.1562- pn)成正比，因为当(0.1562- pn) = 0时，Δpn大约为零。换句话说，我们推断出

其中r > 0是某个常数。实际上，前面描述的逻辑曲线拟合给出了r = 4.417的值。因此，我们得到公式：

可以从p0 = 0.0031数值解得p1 = 0.0052，然后从p1得到p2，依此类推。

（2）指数衰退

[1]例1：理想指数衰减模型

现实的“指数衰减”模型支配着放射性物质质量m的自发衰减（Nagy 2011；Sharon和Sharon 2021）。在这种情况下，从一个特定的初始质量m0开始，每过一段时间H，其质量就会减半。这里，H被称为物质的半衰期[每种放射性物质的一个特性，表示其质量减半所需的时间，以年为单位]。

因此，初始质量m0在时间t = H时减半至m0/2，然后在时间t = 2H时再次减半至m0/4，如此“无限循环”。实际上，这种指数衰减是相当精确的，放射性物质的半衰期可以从几千年到几百万年不等。例如，用于估计文物年龄的碳-14的半衰期约为5730年，而铀-235的半衰期约为7.1亿年！

正如你可以从之前的例子中推断出的那样，放射性物质相当精确的指数衰减过程是由以下方程推出：

当时间以年为单位表示时，我们可以通过以下方式计算放射性物质的质量衰减：

如此一来，我们也可以表示出指数衰减的离散时间方程：

k一样可以计算为：

我们使用python进行可视化模拟：

import matplotlib.pyplot as plt
import numpy as np

# 定义初始质量和半衰期
m0 = 100  # 初始质量
H = 5730  # 半衰期（例如，碳-14的半衰期）

# 生成时间序列（例如，接下来的10个半衰期）
t = np.arange(0, 10 * H, H)

# 计算每个时间点上的质量
m = m0 * (1/2) ** (t / H)

# 绘制指数衰减曲线
plt.plot(t, m)
plt.xlabel('Time (years)')
plt.ylabel('Mass')
plt.title('Exponential Decay Curve')
plt.show()

[2]感染康复模型

现在我们从生物学角度考虑另一个指数衰减的例子，涉及人群中感染者的康复情况。I 表示感染者的数量。在感染者群体中，每天有 1/4 的人康复，我们将这个值称为 λ = 0.25。因此，我们的时间单位是一天，我们的感染者康复模型变为：

我们使用R进行可视化：

# 设置初始参数
I0 <- 100 # 初始感染者数量
lambda <- 0.25 # 每天康复的比例
days <- 0:10 # 天数范围

# 计算每天的感染者数量
infected <- numeric(length(days))
infected[1] <- I0
for (i in 2:length(days)) {
  infected[i] <- infected[i-1] * (1 - lambda)
}

# 绘制图表
plot(days, infected, type = "l", xlab = "Days", ylab = "Number of Infected", main = "Exponential Decay of Infected Individuals", col = "blue", lwd = 2)

# 添加图例
legend("topright", legend = "Infected", col = "blue", lty = 1, lwd = 2)