易感-感染-易感(SIS)模型的实现方式

       连续时间SIS模型与算法

        易感-感染-易感模型,这个模型最早在 1932 年由 Kermack 和 McKendrick 提出。一般用来描述一些治愈后可再次感染的疾病,如流感与肺结核等疾病。 SIS模型的传播过程可以表述为:易感者(S)通过接触其感染者(I)邻居感染疾病,而感染疾病后的患者会成为新的感染源,同时感染者也会随着时间的推移治愈疾病成为一个易感者。其连续时间模型可以见下图:

 其中β为感染速率,μ为恢复速率。对于连续时间的SIS模型,通常采用Gillespie算法,即遗传算法来进行模拟,该算法是基于独立发生的事件,这意味着动力学的每次演变都是由事件决定的。而事件的时间连续性由系统状态的保持时间 𝜏 所分割。具体来说:

(1) 在任意时刻 𝑡,我们根据计算每个个体的转移率 𝜆𝑖(𝑡),即状态为 S 和状态为 𝐼 的个体, 𝜆𝑖(𝑡) 分别等于 βm 和 μ ,m 表示 i 节点周围的感染者邻居数量。

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

(3) 在下一时刻 𝑡 + 𝜏 时,以概率 𝜆𝑖(𝑡)/𝜔(𝑡) 采样选择节点发生状态改变。

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

C++代码实现

        注:下面代码中网络是直接在模拟中生成的,采用的ER随机网络。通常不建议这样使用,因为每次进行模拟,系统的网络结构就发生了变化。一般来说都是将网络单独生成,再通过导入的方法放入动力学程序当中。

#include <iostream>
#include <ctime>
#include <stdlib.h>
#include <fstream>
#include <cmath>
#include "mt19937ar.c"   //从我文章中ER网络实现获取
#include <algorithm>
#include <vector>

using namespace std;


int main()
{
    int Network_Size = 1000;
    double r = 0.4;
    double g = 0.6;
    int I0 = 100;
    double tend = 100;
    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_ER_ext;
    outfile_ER_ext.open("ER_GA_ext.txt");

    double average_degree = 4;
    int edge = 0.5*Network_Size*average_degree;

    for(;;)       //生成ER随机网络
    {
        double judge = 0;
        int random_1 = genrand_int32()%Network_Size;
        int random_2 = genrand_int32()%Network_Size;
        if(random_1 != random_2)
        {
            for(int i = 0; i < adjlist[random_1].size(); i++)
            {
                if(adjlist[random_1][i] == random_2)
                {
                    judge = 1;
                    break;
                }
            }
            if(judge == 1)
            {
                continue;
            }
            else
            {
                adjlist[random_1].push_back(random_2);
                degree[random_1]++;
                adjlist[random_2].push_back(random_1);
                degree[random_2]++;
                edge--;
            }
        }
        if(edge == 0)
            break;
    }

    vector<vector<int> > X_vec;                    //记录所有独立模拟的时间序列
    int n_intervals = tend;
    double interval = tend/(1.0*n_intervals);      //划分取值区间
    vector<int> X_vec_ind(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;
        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] == 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[0] = I;
        t = dt;
        for(;;)
        {
            double u = genrand_real2();
            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[count1] = I;
                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;
            }



            for(int i = 0; i < Network_Size; i++)      //轮盘法选择节点发生状态转变
            {
                if(u < nk[0])
                {
                    if(I_node[0] == 1)
                    {
                        I_node[0] = 0;
                        I--;
                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[0][j]]--;
                        }
                        break;
                    }
                    if(I_node[0] == 0)
                    {
                        I_node[0] = 1;
                        I++;

                        for(int j = 0; j < adjlist[0].size(); j++)
                        {
                            kinf[adjlist[0][j]]++;
                        }

                        break;
                    }
                }

                if(u > nk[i-1] && u < nk[i])
                {
                    if(I_node[i] == 1)
                    {
                        I_node[i] = 0;
                        I--;
                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[i][j]]--;
                        }

                        break;
                    }

                    if(I_node[i] == 0)
                    {
                        I_node[i] = 1;
                        I++;
                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[i][j]]++;
                        }
                        break;
                    }
                }

                if(I == 0)
                {
                    while (count1 <= n_intervals)
                    {
                        X_vec_ind[count1] = 0;
                        count1 ++;
                    }
                    break;
                }

            }

            wt = 0;
            for(int i = 0; i < Network_Size; i++)
            {
                if(I_node[i] == 1)
                {
                    rate[i] = g;
                }
                if(I_node[i] == 0)
                {
                    rate[i] = r*kinf[i];
                }
                wt = wt + rate[i];
            }
            dt = -log(1 - v)/wt;
            t = t + dt;
        }
     X_vec.push_back(X_vec_ind);
    }

        double EXt;
        for(int i=0; i <= n_intervals; i++)        //取时间序列的平均值
        {
            EXt = 0;
            for(int j=0; j < counter; j++)
                EXt += X_vec[j][i];

            EXt /= 1.0*counter*Network_Size;
            outfile_ER_ext << EXt << "\n";
        }
        outfile_ER_ext<< endl;




}

离散时间SIS模型与算法

        而离散时间中,转变概率替代了转变速率,事件的转变时间由时间间隔 Δ𝑡 控制,其模型如下:

 其中β'与μ'分别表示感染概率与恢复概率,值得注意的是,在大多数离散时间模型研究中,人们习惯使用β'=βΔ𝑡,μ'=μΔ𝑡,使用这种方法的一个弊端是连续时间与离散时间模型在参数较大时,两者无法映射(Physical Review E, 2016, 94(5):052125)

        离散时间模型算法通常采用同步更新方法:每个时间步骤中,所有个体都以同步方式更新其状态。状态为 S 的个体以概率 1-(1-β')^m 转变为状态 𝐼,而状态为 𝐼 的个体以概率 μ' 转变为状态 S 。整个过程迭代到系统达到稳定或系统中 𝐼 个体为零时停止。

C++算法实现

#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 = 10000;
    double average_degree = 4;
    double r = 0.2;
    double g = 0.48;
    double dt = 0.1;
    double tend = 100;
    double counter = 100;
    int I0 = 100;

    vector<int> degree(Network_Size);
    unsigned long idum;
    idum=(unsigned long)time(NULL);
    init_genrand(idum);
    vector<vector<int> > adjlist(Network_Size);
    ofstream outfile_ER_S_ext;
    double r1;
    outfile_ER_S_ext.open("ER_S_ext.txt");
    double Sum_degree = 0;
    int edge = 0.5*Network_Size*average_degree;
    for(;;)
    {
        double judge = 0;
        int random_1 = genrand_int32()%Network_Size;
        int random_2 = genrand_int32()%Network_Size;
        if(random_1 != random_2)
        {
            for(int i = 0; i < adjlist[random_1].size(); i++)
            {
                if(adjlist[random_1][i] == random_2)
                {
                    judge = 1;
                    break;
                }
            }
            if(judge == 1)
            {
                continue;
            }
            else
            {
                adjlist[random_1].push_back(random_2);
                degree[random_1]++;
                adjlist[random_2].push_back(random_1);
                degree[random_2]++;
                edge--;
            }
        }
        if(edge == 0)
            break;
    }


    double rs;
    double gs;
    double t = 0;
    rs = r*dt;
    gs = g*dt;
    int I;
    int random_1;    //用来初始化感染者节点
    vector<int> kinf(Network_Size,0);
    vector<int> I_node(Network_Size,0);
    vector<int> kinf_old(Network_Size);


    int n_intervals = tend/dt;
    double interval = tend/(1.0*n_intervals);
    vector<int> X_vec_ind(n_intervals + 1);
    vector<vector<int> > X_vec;
    int count1;


    for(int n = 0; n < counter; n++)
    {

        for(int i = 0; i < Network_Size; i++)
        {
            kinf[i] = 0;
            I_node[i] = 0;
        }
        I = I0;
        for(;;)         //随机取初始感染者节点
        {
            if(I == 0)
            {
                break;
            }
            random_1 = genrand_int32()%Network_Size;
            if(degree[random_1] == 0)
            {
                continue;
            }
            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]++;
                }
            }
        }

        t = 0;
        I = I0;
        X_vec_ind[0] = I;
        count1 = 1;


        for(;;)
        {
            for(int i = 0; i < Network_Size; i++)
            {
                kinf_old[i] = kinf[i];
            }
            t = t + dt;
            while((t > count1*interval)&&(count1 <= n_intervals))
            {
                X_vec_ind[count1] = I;
                count1 ++;
            }

            if(t > tend)
                break;
            if(I == 0)
            {
                while (count1 <= n_intervals)
                {
                    X_vec_ind[count1] = 0;
                    count1 ++;
                }
                break;
            }

            double u;

            for(int i = 0; i < Network_Size; i++)
            {
                u = genrand_real2();

                if(I_node[i] == 1)
                {
                    if(u < gs)
                    {
                        I_node[i] = 0;
                        I--;

                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[i][j]]--;
                        }
                    }

                }
                else if(I_node[i] == 0)
                {
                    if(u < 1 - pow((1 - rs),kinf_old[i]))
                    {
                        I_node[i] = 1;
                        I++;
                        for(int j = 0; j < adjlist[i].size(); j++)
                        {
                            kinf[adjlist[i][j]]++;
                        }
                    }
                }
            }
        }
        X_vec.push_back(X_vec_ind);
    }

        double EXt;
        for(int i=0; i <= n_intervals; i++)
        {
            EXt = 0;
            for(int j=0; j < counter; j++)
                EXt = EXt + X_vec[j][i];

            EXt = EXt/(counter*Network_Size);
            outfile_ER_S_ext  << EXt << "\n";
        }
        outfile_ER_S_ext << endl;




    return 0;
}

结语

        对于各种网络的实现方法大同小异,但对于不同的网络需要单独优化程序,例如在全连接网络中,每个节点的感染者邻居数是相同的,因此可以极大的简化抽样的步骤。上面这两个程序是我刚入门时参考论文中大佬给出的程序写的,在后面的独立研究中也很少使用到,如有错误,请指正,谢谢。

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值