SEIR流行病动力学模型的C++实现

SEIR简要描述

        对于SEIR动力学，即易感-潜伏-感染-康复动力学过程，与易感-感染-康复(SIR)动力学不同的是，该动力学中多了潜伏态过程，即个体被感染疾病后，需要经过一定的时间才能表现疾病症状，在此期间我们称为潜伏态。

        该动力学过程可以描述为：易感态个体接触周围患病个体后，通过感染速率β感染疾病，并且进入疾病潜伏态，随后通过速率λ转变为感染态个体，而感染个体会以μ的康复速率恢复为康复者，见下图。

        SEIR的疾病过程可以用来描述很多带有潜伏期的疾病，如梅毒、狂犬病以及最近肆虐全球的COVID-19等。由于其及其模型变体在预测新冠传播方面有非常多的应用，因此这个模型也成为了现今的研究热点。我最近在与医科同学交流时，了解到狂犬病的症状与SEIR模型较为相似，因此想来做一些这方面的工作，所以写下了这个程序供大家参考。对于SEIR模型的性质，及数学模型，大家可以在很多文献或书籍中得到一个更全面的认知，我也不做介绍，直接上代码。

 C++代码实现

        采用Gillespie算法进行实现，具体来说：

        (1)在任意时刻 𝑡，我们算每个个体的转移率 𝜆𝑖(𝑡)，即状态分别为 0 ，1，2，3(易感，感染，康复，潜伏)的个体， 速率𝜆𝑖(𝑡) 分别为βkinf(i)，μ，0，λ，其中kinf(i)为节点的感染者邻居数。

        (2) 然后计算出系统的总速率 𝜔(𝑡) = ∑𝑖 𝜆𝑖(𝑡)，其中系统状态保持时间 𝜏 = 1/𝜔(𝑡)。

        (3) 在下一时刻 𝑡 + 𝜏 时，以概率 𝜆𝑖(𝑡)/𝜔(𝑡) 采样选择节点发生状态改变，选择节点为易感态S时，转变为潜伏态E；选择节点为潜伏态E时，转变为感染态I；选择节点为感染态I时，转变为康复态R；选择康复态R时，跳过该节点，选择其他节点。

        (4) 重复 (1) 到 (3) 过程，直到整个系统达到稳定状态或系统中 𝐼 个体为零时停止迭代。

#include <iostream>
#include <ctime>
#include <stdlib.h>
#include <fstream>
#include <cmath>
#include "mt19937ar.c"
#include<algorithm>
#include <vector>

using namespace std;


int main()
{
    int Network_Size = 100000;
    double r = 0.4;                //感染速率
    double g = 0.2;               //恢复速率
    double k = 0.4;              //潜伏态转变速率
    int I0 = 100;
    double tend = 50;
    double counter = 10;
    vector<int> degree(Network_Size,0);
    vector<vector<int> > adjlist(Network_Size);

    unsigned long idum;
    idum=(unsigned long)time(NULL);
    init_genrand(idum);
    ofstream outfile_S;
    outfile_S.open("SEIR_S.txt");
    ofstream outfile_E;
    outfile_E.open("SEIR_E.txt");
    ofstream outfile_I;
    outfile_I.open("SEIR_I.txt");
    ofstream outfile_R;
    outfile_R.open("SEIR_R.txt");
    ifstream infile;
    infile.open("ER_degree1.txt");
    ifstream infile2;
    infile2.open("ER_adjlist1.txt");
    double Sum_degree = 0;

    for(int i = 0; i < Network_Size; i++)          //导入节点度与网络结构
    {
        infile >> degree[i];
        Sum_degree = Sum_degree + degree[i];
        adjlist[i].resize(degree[i]);
    }
    for(int i = 0; i < Network_Size; i++)
    {
        for(int j = 0; j < degree[i]; j++)
        {
            infile2 >> adjlist[i][j];
        }
    }


    vector<vector<int> > X_vec_S;
    vector<vector<int> > X_vec_E;
    vector<vector<int> > X_vec_I;
    vector<vector<int> > X_vec_R;
    int n_intervals = tend;
    double interval = tend/(1.0*n_intervals);
    vector<int> X_vec_ind_S(n_intervals + 1);
    vector<int> X_vec_ind_E(n_intervals + 1);
    vector<int> X_vec_ind_I(n_intervals + 1);
    vector<int> X_vec_ind_R(n_intervals + 1);
    vector<double> rate(Network_Size,0);
    vector<int> I_node(Network_Size,0);
    vector<double> sum(Network_Size,0);
    vector<double> nk(Network_Size,0);
    vector<int> kinf(Network_Size,0);

    for(int n = 0; n < counter; n++)
    {
        double wt = 0;
        double dt;
        double t = 0;
        double I;
        double R = 0;
        double E = 0;
        int random_1;    //用来初始化感染者节点
        int count1;
        for(int i = 0; i < Network_Size; i++)
        {
            rate[i] = 0;
            I_node[i] = 0;
            kinf[i] = 0;
        }
        I = I0;
        for(;;)         //随机取初始感染者节点
        {
            if(I == 0)
            {
                break;
            }
            random_1 = genrand_int32()%Network_Size;
            if(I_node[random_1] == 0)
            {
                I_node[random_1] = 1;
                I--;
            }
            else
                continue;
        }
        for(int i = 0; i < Network_Size; i++)
        {
            for(int j = 0; j < adjlist[i].size(); j++)
            {
                if(I_node[adjlist[i][j]] == 1)
                {
                    kinf[i]++;
                }
            }
        }
        for(int i = 0; i < Network_Size; i++)
        {
            if(I_node[i] == 1)
            {
                rate[i] = g;
            }
            if(I_node[i] == 2)
            {
                rate[i] = 0;
            }
            if(I_node[i] == 3)
            {
                rate[i] = k;
            }
            if(I_node[i] == 0)
            {
                rate[i] = r*kinf[i];
            }
            wt = wt + rate[i];
        }
        double v = genrand_real2();
        dt = -log(1 - v)/wt;
        I = I0;
        count1 = 1;
        X_vec_ind_I[0] = I;
        X_vec_ind_S[0] = Network_Size-E-I-R;
        X_vec_ind_E[0] = E;
        X_vec_ind_R[0] = R;
        t = dt;
        for(;;)
        {
            //cout << t << '\t' << E << '\t' << I << '\t' << R << endl;
            double u;
            v = genrand_real2();
            for(int i = 0; i < Network_Size; i++)
            {
                nk[i] = 0;
                sum[i] = 0;
            }
            while((t > count1*interval)&&(count1 <= n_intervals))
            {
                X_vec_ind_E[count1] = E;
                X_vec_ind_I[count1] = I;
                X_vec_ind_R[count1] = R;
                X_vec_ind_S[count1] = Network_Size-I-R-E;
                count1 ++;
            }
            if(t > tend)
            {
                break;
            }

            for(int i = 0; i < Network_Size; i++)
            {
                if(i == 0)
                {
                    sum[i] = rate[i];
                }
                else
                {
                    sum[i] = sum[i-1] + rate[i];
                }
                nk[i] = sum[i]/wt;
            }
            u = genrand_real2();
            for(int i = 0; i < Network_Size; i++)
            {

                if( u < nk[i])
                {
                    if(I_node[i] == 1)
                    {
                        I_node[i] = 2;
                        I--;
                        R++;
                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[i][j]]--;
                        }
                        break;
                    }
                    if(I_node[i] == 2)
                    {
                        continue;
                    }
                    if(I_node[i] == 3)
                    {
                        I_node[i] = 1;
                        E--;
                        I++;
                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[i][j]]++;
                        }
                    }
                    if(I_node[i] == 0)
                    {
                        I_node[i] = 3;
                        E++;
                        break;
                    }
                }
                if(I == 0)
                {
                    while (count1 <= n_intervals)
                    {
                        X_vec_ind_I[count1] = 0;
                        X_vec_ind_E[count1] = E;
                        X_vec_ind_R[count1] = R;
                        X_vec_ind_S[count1] = Network_Size-E-I-R;
                        count1 ++;
                    }
                    break;
                }

            }

            wt = 0;
            for(int i = 0; i < Network_Size; i++)
            {
                if(I_node[i] == 1)
                {
                    rate[i] = g;
                }
                if(I_node[i] == 2)
                {
                    rate[i] = 0;
                }
                if(I_node[i] == 3)
                {
                    rate[i] = k;
                }
                if(I_node[i] == 0)
                {
                    rate[i] = r*kinf[i];
                }
                wt = wt + rate[i];
            }
            dt = -log(1 - v)/wt;
            t = t + dt;
        }
     X_vec_E.push_back(X_vec_ind_E);
     X_vec_I.push_back(X_vec_ind_I);
     X_vec_R.push_back(X_vec_ind_R);
     X_vec_S.push_back(X_vec_ind_S);
    }

    double EXt1,EXt2,EXt3,EXt4;
    for(int i=0; i <= n_intervals; i++)
    {
        EXt1 = 0;
        EXt2 = 0;
        EXt3 = 0;
        EXt4 = 0;
        for(int j=0; j < counter; j++)
        {
            EXt1 += X_vec_S[j][i];
            EXt2 += X_vec_E[j][i];
            EXt3 += X_vec_I[j][i];
            EXt4 += X_vec_R[j][i];
        }

        EXt1 /= 1.0*counter*Network_Size;
        EXt2 /= 1.0*counter*Network_Size;
        EXt3 /= 1.0*counter*Network_Size;
        EXt4 /= 1.0*counter*Network_Size;
        outfile_S << EXt1 << "\n";
        outfile_E << EXt2 << "\n";
        outfile_I << EXt3 << "\n";
        outfile_R << EXt4 << "\n";
    }
    outfile_S << endl;
    outfile_E << endl;
    outfile_I << endl;
    outfile_R << endl;


}

结果

        模拟采用平均度为4的ER随机网络，网络规模为100000，其中感染速率β=0.4，恢复速率μ=0.2，潜伏转变速率λ=0.2，初始感染总数为100。

        若程序或动力学理解有错，希望有大佬能指正，感谢！