如何通过Matlab对新型冠状病毒肺炎的爆发过程进行系统建模

目前全世界已经发生了四起大规模的冠状病毒大暴发，上述爆发包括 2003 年在中国大陆爆发的“严重急性呼吸系统综合症”（SARS），2012 年在沙特阿拉伯的“中东呼吸综合症”（MERS）爆发，2015 年在韩国发生的 MERS 爆发。这些暴发分别导致 8,000 例和 2,200 例确诊的 SARS 和 MERS 病例。最近，在湖北省省会武汉市发生了第四次大暴发。

在病毒连续爆发过程中，医学研究者获得了病毒大量的传播参数。而数学家则通过上述参数对病毒的传播模式进行了大量分析。在病毒传播模型的分析过程中，利用常微分方程综合考虑病毒传播的各种参数可以有效的对病毒传播过程进行模拟。
建立一套方法让一线医学研究者快速理解和分析数学工作者的工作结果至关重要。很多时候，需要人通过一定工作来完成数学模型---方程组----计算机代码之间的鸿沟，我们以最近发表的论文：
Tang B , Wang X , Li Q , et al. Estimation of the Transmission Risk of 2019-nCov and Its Implication for Public Health Interventions[J]. Social Science Electronic Publishing.
为例，用Matlab的常微分方程工具箱（ODE toolbox）为载体,看看如何利用matlab的实时代码功能快速获得数学家研究结果的模拟
首先：该论文利用了通用的 SEIR 型流行病学模型对一次传染性流行病的大爆发进行了研究，模型主要通过对流行病的自然传播过程进行分析，然后结合诸如检疫，隔离和治疗等干预措施，对整个一次流行病爆发进行系统模拟。
为了理解他们的工作，我们现将模型要点分析如下：首先，我们将爆发过程中的武汉居民按照分为易感者（S），暴露者（E），无症状传播者（A），有症状感染者（I），住院者（H）和康复者（R）几个仓室，这些舱室之间的关系如下：

并进一步对人群进行分层包括易感者隔离区（S_q），暴露者隔离区（E_q）和感染者隔离区（I_q）
在政府的持续干预下，比例为q的患者不断被识别出来，并且被干预。被隔离的人进入隔离室S_q和E_q.
其他1-q的患者持续的暴露在病毒当中被感染。最后感染者死亡或者康复，获得该次流行病爆发的免疫力。
详细的参数设定可以看论文的参数列表：

总体上说：
可以用微分方程描述各舱室之间关系如下：

论文中的方程组描述的是各个城市当中各种状态的人数转化的情况，比如第一个公式表示易感者的变化率β与感染者接触c，每次感染传播概率q，以及接触者的隔离率A，还有已经隔离但因为确定无症状，重新回归社区的速率(λ)都有关系。
利用matlab可以方便的模拟上述方程组：
%这是一个模拟1008万人口的城市在疾病爆发的情况下感染人数的爆发模型
%使用的是SEIR模型
%假设从t0时间开始，城市出现感染，到100天的时候结束🔚

to = 0;
tf = 100;
tspan = [to tf];
y0 = [11108100 105.1 27.679 53.839 739 1.1642 1 2];
[t,S] = ode45(@SEIRODE, tspan, y0);
%ode45参数：1.@函数句柄 2.t的取值代表计算周期， 3. 7个初始值
%ode45是用来求解常微分函数的方法

运行上述代码，再利用下面的代码绘图，可以获得以下结果：

plot(t,S)
title('Human Population Without Control')
xlabel('Time')
ylabel('Susceptible, Exposed, Infected, Recovered')
legend('易感者', '暴露人数', '感染人数', '无症状感染者','初步隔离者','隔离暴露者','隔离感染者','康复者')
可以看到，在控制感染者出行的情况下易感者数量在14-16天左右开始下降。而因为有接触被隔离的人的数量还在快速上升。

为了更好的看到上图下方的不明显曲线，我们取消了了上面曲线中数量最为巨大的易感者和隔离者人数，此时曲线变为
figure()
plot(t,S(:,[2,3,4,6,7]))
legend( '暴露人数', '感染人数', '无症状感染者','隔离暴露者','隔离感染者')

最后，我们通过代码获得实现本次模拟的一些结果：
max_p=num2str(max(S(:,2)));
[~,temp]=max(S(:,2));
max_p_data=num2str(t(ceil(temp)));
max_g=num2str(max(S(:,3)));
[~,temp]=max(S(:,3));
max_g_data=num2str(t(ceil(temp)));
max_gl=num2str(max(S(:,7)));
[~,temp]=max(S(:,7));
max_gl_data=num2str(t(ceil(temp)));
formatSpec = "最大暴露人数为: %s 出现时间为第 %s 天 \n " + ...
"最大感染人数为: %s 出现时间为第 %s 天 \n" + ...
"最大隔离感染者人数为： %s 出现时间为第 %s 天";
sprintf(formatSpec,max_p,max_p_data,max_g,max_g_data,max_gl,max_gl_data)

方程组具体参数设置以论文为准。