出一个微分方程组模型（动力学模型）的推送，最具代表性的就是传染病SIR模型，同时万一美赛用上了不会翻译，这里尽量保证给出原汁原味的英文，先来看传染病模型叭：

传染病SIR模型

一个最简单的设定， 代表易感者(susceptibles)， 代表患者(infectives)， 代表康复且不会再得的人群(recovered individuals)

首先我们假设如果没有康复者会咋样：

这个情况就非常简单了，易感者减少的人群数量等于感染者增加人群数量，且由于传染需要易感者和感染者接触，所以和两者人数都有关系：

由于总人数不变，设置 ，对第二个微分方程进行变量替换并求解微分方程得：

假设 时，，可以求解出 ：

因此直接解析式就求完了，简简单单：

解析式都有了，说明这个 SI 问题过于简单，随便带入一下数据就能求出来各个参数从而得到完整的变化曲线，这里不多赘述，我们再来看 SIR 模型，这个就很难直接得到解析式，而要分析不同时间各个类型人群数量的变化趋势啦：

依旧是感染者增加速度和感染者人数和易感者人数都有关系，而感染者痊愈的人数变化速度只和感染者人数有关，

由于：

所以总人数是保持不变的(is closed)（暂时不考虑死亡），因此可以设置 ,对原来公式进行无量纲化(Non-dimensionalise)(就是去除单位，比如原本是90感染者，变为90%感染者)，对其进行如下变量替换：

可以得到：

其中是基本繁殖率(basic reproductive ratio)是记录传染病传染能力的重要属性，

由于前两个微分方程与 没有关系，我们只研究 ，可以通过 进行计算，注：有的时候传染病模型是考虑复感者的，比如有些传染病SIR模型是这样的， 代表复感率(echoing ratio)：

这种情况就可以使用把第二个式子中的进行替代，这个过程叫做解耦 (decouple)

但这篇文章就先不考虑那么多，这里考虑无复感者：

然后我们可以通过分析绘制相图(phase portrait)

相图绘制分析过程

为了达到稳态(steady states), 所有变化率都需要为 0:

如果 :

如果 :

我们发现只要 ，无论 取何值都是满足条件的

而对于区域(为了保证物理意义，人数必须是正数)(physical sense) , 发现总有

对于区域 , 总有 , 而区域 , 有 , 同时, 如果 , 有

如果 , , 如果 ,  

因此, 对于 和 , 可分别绘制相图:

这里相图范围是三角形是因为

可以用MATLAB试一下，使用NAN填充三角形外面区域：

R0 = 3;
[S,I] = meshgrid(0:.005:1);
dS = - R0.*S.*I;
dI = R0.*S.*I - I;

dS(S+I>1) = 0;
dI(S+I>1) = 0;

S = streamslice(S,I,dS,dI,5);
axis([0,1, 0,1])

我们想知道如果初始 ，且大部分人都是易感者 (相图里初始点在 那条直线上最右侧即(0,1)点附近) 那么疫情最后有多少人会被感染过(size of the epidemic, i.e. how many will suffer from the disease)

从相图来看如果那么最终有人不会被感染，即最后:

那么事实是这样嘛，我们要用公式来确认，首先我们有：

我们对式子1进行常微分方程求解：

带入初始条件 可得：

由于这个模型不考虑复感，因此最终感染者一定会被清零，假设最终痊愈者数量为 ，我们有：

带入前式有：

我们把等式两边的曲线画在同一图像里：左侧是，右侧是 情况:

我们发现 的时候有解 (When R0 > 1 there is a single positive root, so some susceptibles can remain in the population after the epidemic has passed, with the total (dimensional) number in the removed class is NR∞, where R∞ is found as the positive root above.)

写段MATLAB代码绘图试试（这篇写的太赶了来不及修饰图像了）:

tspan = 0:.001:10;
SIR0 = [1-1e-3; 1e-3; 0];
[t, SIR] = ode45(@(t, SIR) funSIR(t, SIR), tspan, SIR0);

hold on
plot(tspan, SIR(:,1), 'LineWidth',2)
plot(tspan, SIR(:,2), 'LineWidth',2)
plot(tspan, SIR(:,3), 'LineWidth',2)
legend({'S','I','R'})

function dSIR = funSIR(~, SIR)
R0 = 3;
S  = SIR(1);
I  = SIR(2);
R  = SIR(3);
dS = - R0.*S.*I;
dI = R0.*S.*I - I;
dR = I;
dSIR = [dS; dI; dR];
end

可以看到最后 也没变为0.

平流反应扩散模型(Reaction-advection-diffusion equation)的行波解(traveling waves)

这题主要目的就是将一个和位置和时间有关的方程组变成只和 有关的解，并绘制一个稳态到另一个稳态之间的连线。 用来分析各个稳态的稳定性和状态变化趋势。这里实际不是一个完全的平流反应扩散模型，因为这里没有考虑平流。

我们采用一个简单的模型来表示酵母在葡萄糖上的生长，假设这是一个一维空间模型。令表示位置和时间时刻酵母的浓度，表示葡萄糖的浓度。我们将取正的常数、和分别表示葡萄糖的临界浓度、葡萄糖的扩散系数以及酵母的生长速率常数。我们假设酵母的运动可以忽略不计，因此只有葡萄糖能够扩散。描述酵母生长及相应动态的方程为

• 表示能量转换率的倒数，它决定了酵母消耗葡萄糖的效率与其生长速率之间的关系。它表示酵母浓度增加一个单位所需的葡萄糖浓度。
• 只有当 大于 时，酵母才能生长。
• 表示葡萄糖和酵母数量共同影响的酵母生长速率。
• 表示葡萄糖的消耗速率。
• 表示葡萄糖的扩散。

行波解假设及变量替换

将行波假设 和 代入，其中 且波速 ，从而得到关于 和 的二阶常微分方程（ODE）系统，其中 是自变量。(Substitute the travelling wave ansatz and with and constant wave speed to obtain a second-order system of ordinary differential equations (ODEs) for and in terms of )

以下请忽略我的垃圾英语，实在是没空把我小作业里的英文都改了，但是这些因果叙述关系应该还是很好懂的：

It is easy to obtain:

Substitution variables:

假设有初值条件 和 ，将其转换为一阶常微分系统(first-order system of ODEs)：

Consider:

Intergal from each side of equation:

Take the initial condition and :

Therefore:

Which means:

解的稳定性分析

通过变化率为0求出稳态，通过计算雅可比矩阵特征值判断稳定性，两个特征值都是负数说明每次点有往外走的趋势，就会施加一个反方向的“力”，导致点再走回去，这就是特征值都是负数这个稳态是稳定的原因：

When the system reaches its steady states, the derivative terms turn into 0:

To satisfy the first equation, we have

• When

• When

So the steady states are and

Consider the jacobian matrix:

• When  

This is a lower triangular matrix, so the characteristics are equal to the corner elements :

As , then which is stable.

• When  

We have the characteristic function:

easy to obtain:

As , then , which means , which is unstable.

Thus the stable state is unstalbe, the stable state is stalbe.

绘制相图

If we take , so we have

If :

If :

If :

So on this line,

For the area , , for the area , .

Then we can draw the phase plane portrait:

MATLAB验证一下：

k = 1; s = 1; a = 1; D = 1; C0 = .5;
[C,G] = meshgrid(0:.005:2);
dC = - k./s.*C.*G;
dG = - s./D.*G - a.*s./D.*(C - C0);

S = streamslice(C,G,dC,dG,25);
axis([0,.6, 0,.6])

hold on
plot([C0,0], [0,a*C0],'LineWidth',1,'Color','k')

酶动力学(Enzyme kinetics)激活与抑制模型(Activation and inhibition)及协同效应(Cooperative phenomena)

我们首先介绍数学生物学中的一个基本原理——质量作用定律(Law of mass action)。该定律指出，如果物质A（例如个体或化学物质）与物质B反应生成物质C，反应过程如下所示：

即A和B结合形成C，那么反应速率由 给出（其中和分别表示相应的浓度或种群密度），因此有：

其中是反应的速率常数。通常，反应是可逆的(reversible)，因此我们写作：

从而得到以下方程：

酶动力学(Enzyme kinetics)基础

考虑一种酶 (enzyme)，它催化底物 (substrate)转化为产物 (product)。该过程可表示为：

即酶 和底物 结合生成中间产物 , 再转化为酶的最初状态和产物 的过程， 根据质量作用定律，我们得到以下方程组：

然后可以注意到，以下关系成立：

因此，可以设置：

表示酶的守恒(conservation of enzyme)，其中是酶的总量(total amount of enzyme)。类似地，

因此，

表示底物的守恒(conservation of substrate)，其中是最初可用的底物总量。

我们希望通过可用的底物浓度来确定反应速率 。通常，底物的量远大于酶的量（），因此我们可以预期酶将以最大效率工作，假设符合酶动力学中的准稳态假设(Quasi-Steady-State Assumption，QSS)(该假设的核心思想是，反应中某些中间产物（通常是酶-底物复合物）的浓度在时间上变化非常缓慢，因此可以认为它们在反应过程中保持近似恒定)，因此有：

因此，

其中。

然后，由于

因此，我们得到

我们就得到了反应的速率，即产物的生成速度。这被称为米哈利斯-门腾动力学(Michaelis-Menten kinetics)，其中是米哈利斯常数(Michaelis constant)，只和酶转换为中间产物的速率常数和转化回酶原来状态的速率常数有关，有些人可能会想，嗯？？？？？最终方程和酶相关的部分只和酶还有中间体总量有关，和各自的量没有关系吗？？？？实际上由 可以计算出各自的量，因此各自的量的信息实际已经包含在公式里了。

激活与抑制模型(Activation and inhibition)

考虑一种酶 ，它催化底物 转化为产物 的反应：

但它也可以与修饰物 (modifier)反应：

我们还假设 和 可以同时与酶结合，且可能具有不同的速率常数：

我们使用质量作用定律写出一系列常微分方程（ODEs）来描述这些反应：

我们注意到与其他变量解耦，因此我们得到：

即酶的守恒，

即底物的守恒，

即修饰物的守恒。

如果某个结合分子（称为配体）(ligand)使得另一个分子无法与酶结合，即我们设置 ，，，这种现象被称为竞争性抑制(competitive inhibition)。我们希望了解修饰物对反应速率的影响。此时：

同时有：

其中， 和 分别表示酶和底物的初始量。由于酶很少，假设 ，由酶动力学中的准稳态假设(Quasi-Steady-State Assumption，QSS)，将以下式子设置为零：

由之前的微分方程设置导数项为0得 , 得到：

我们可以将其写为：

其中，

是底物的米哈利斯常数(Michaelis constant)，

是修饰物的平衡常数(equilibrium constant)。得到，

因此，我们得到：

以及

因此，

而我们之前不存在修饰物时候的反应速率为：

即修饰物改变了酶的米哈利斯常数，将其从 改为 ，在给定底物浓度 的情况下降低了反应速率，同时不改变最大反应速率。

协同效应(Cooperative phenomena)

我们可能遇到一种协同酶(cooperative enzyme)，其中配体(ligand)可以结合到酶的多个位点(sites)上，并且其他已结合的配体的存在可能会改变反应速率常数。

对于具有两个结合位点的酶，我们可以依次结合两个分子，即：

其中因我们有两个位点进行结合/解离/生成(bind/unbind/produce)产物，因此反应速率常数前有一个因子 2， 和 分别表示一个和两个位点被占据的复合物。如果酶不是协同的，那么无论另一个位点是否被占据，反应速率都是不会发生改变，即 ，，。

使用质量作用定律，我们可以得到以下一系列常微分方程：

这给出了酶和底物的守恒定律：

假设 ，我们使用准稳态假设，设置：

还是把之前微分方程对应左侧设置为0我们得到：

这些可以写为：

其中：

是米哈利斯常数。

因此有：

所以：

得到：

非协同(non-cooperative case)

在非协同情况下，我们有 和 ，得到：

这和标准的米哈利斯-门腾动力学公式很像，只不过反应速率是原来的两倍。

高度协同(highly cooperative limiting case)

在高度协同的极限情况下，当第二个分子的结合受到第一个分子存在的极大帮助时， 的浓度会非常低，而 和 的浓度则较高。考虑以下方程：

这两个方程可以结合得到：

如果我们假设 和 的浓度相同，而 的浓度较低，那么根据方程我们得到 ，根据方程 得到 ，并且根据方程 我们可以推断 和 的大小相同。因此，这给出了以下极限：

我们发现和上面的无协同公式：

非常像，我们对这俩公式做一个推广，用于模拟协同酶的行为，就一个是1一个是2嘛，我们直接改成n：

这被称为希尔方程（Hill function）。这里的n指的就是希尔系数，有以下解释：

• ．正协同结合：一旦一个配体分子与酶结合，其对其他配体分子的亲和力就会增加。例如，氧气与血红蛋白结合的希尔系数 （正协同作用的一个例子）在 1．7－3．2 范围内。
• ．负协同结合：一旦一个配体分子与酶结合，它对其他配体分子的亲和力就会降低。
• ．非协同（完全独立）结合：酶对配体分子的亲和力不依赖于其他配体分子是否已经结合。当 时，我们得到一个可以用 Michaelis－Menten 动力学建模的模型。

结

本文讲了三个微分方程组构成的模型，希望能对大家美赛有帮助，整理不易点个赞叭，~~