SIRS 模型与 NetworkX 上的代理

Question

我想在 NetworkX 图上模拟和可视化 SIRS 模型（恢复/删除的节点可能再次变得容易受到影响）。每个节点也作为一个代理运行，并且可以在每个时间步选择以概率 p 进行隔离，因此无法在该时间步内被感染。

我已经用 EoN.Gillespie_simple_contagion 模拟了一个 SIRS 模型，我试图了解我是否可以修改这些方法以包含代理行为，或者我是否需要编写一个定制方法。

我已经看到可以修改一些 EoN 方法：Modified SIR model

这是我正在使用的代码：

import networkx as nx
import matplotlib.pyplot as plt 
import EoN
from collections import defaultdict

# parameters required for the SIRS model
a = 0.1
b = 0.01
y = 0.001
d = 0.001

# Simple contagions
# the below is based on an example of a SEIR disease (there is an exposed state before becoming infectious)
# from https://arxiv.org/pdf/2001.02436.pdf

Gnp = nx.gnp_random_graph(500, 0.005)

H = nx.DiGraph() #For the spontaneous transitions
H.add_edge('I', 'R', rate = b)  # an infected node can be recovered/removed
H.add_edge('I', 'S', rate = y)  # an infected node can become susceptible again
H.add_edge('R', 'S', rate = d)  # a recovered node can become suscepticle again

J = nx.DiGraph() #for the induced transitions
J.add_edge(('I', 'S'),('I', 'I'), rate = a)  # a susceptible node can become infected from a neighbour
IC = defaultdict(lambda: 'S')

# set all statuses except one to susceptible. only one node shall begin infected
for node in range(500):
    IC[node] = 'S'
IC[0] = 'I'

return_statuses = ('S', 'I', 'R')
print('doing Gillespie simulation')

t, S, I, R = EoN.Gillespie_simple_contagion(Gnp, H, J, IC, return_statuses, tmax = 500)

print('done with simulation, now plotting')
plt.plot(t, S, label = 'Susceptible')
plt.plot(t, I, label = 'Infected')
plt.plot(t, R, label = 'Recovered')
plt.xlabel('$t$')
plt.ylabel('Number of nodes')
plt.legend()
plt.show()

Answer 1

可以为所需功能编写定制方法：

def SIRS_simulation_with_quarantine(graph, params, tmax):
    a, b, y, d, q = params
    t=0
    graph_current_timestep = graph.copy()
    graph_previous_timestep = graph.copy()
    Infected_nodes = [x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='I']
    S = []
    I = []
    R = []
    T = []
    
    while(t<tmax):
        
        for node in graph.nodes():
            
            # initialise random probablilty for that node
            p = np.random.rand()
            # determine neighbours for that node
            neighbors = [n for n in graph_previous_timestep.neighbors(node)]
            
            # depending on the state, check for any state changes and update node attributes
            if graph_previous_timestep.nodes[node]['state'][0] == 'S':
                if (q >  np.random.rand()):                                             # if the node didnt quarantine...
                    if len(set.intersection(set(neighbors), set(Infected_nodes)))>0:    # if any neighbors are infected...
                        if (p < a):                                                     
                            attrs = {node: {'state': 'I'}}
                            nx.set_node_attributes(graph_current_timestep, attrs)                   
            elif graph_previous_timestep.nodes[node]['state'][0] == 'I':                
                if (p < b):
                    attrs = {node: {'state': 'R'}}
                    nx.set_node_attributes(graph_current_timestep, attrs)                    
                elif (p < y):
                    attrs = {node: {'state': 'S'}}
                    nx.set_node_attributes(graph_current_timestep, attrs)

            elif graph_previous_timestep.nodes[node]['state'][0] == 'R':                
                if (p < d):
                    attrs = {node: {'state': 'S'}}
                    nx.set_node_attributes(graph_current_timestep, attrs)

        # after the changes to states have occured for all nodes, record number of nodes at each state, increment timestep
        Infected_nodes = [x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='I']
        S.append(len([x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='S']))
        I.append(len(Infected_nodes))
        R.append(len([x for x in graph_current_timestep.nodes() if graph_current_timestep.nodes[x]['state'][0]=='R']))
        T.append(t)      
        graph_previous_timestep = graph_current_timestep.copy()
        t = t+1
    
    return T, S, I, R

Answer 2

为此，我们将引入一个新的节点状态'X' 来表示易感隔离节点。我在'H' 中添加了一条从'S' 到'X' 的边，以及从'X' 到'S' 的另一条边。我已将'X' 添加到return_statuses，并绘制了它。

为了简化您的代码，我删除了最初将所有内容分配为 'S' 的循环。因为它使用了default_dict，所以一切都隐含地以'S' 开头。

您提到了下一个“时间步长”并以概率p 移动。这不是这段代码中会发生的事情。该代码在连续时间内工作。所以它需要一个费率。它没有具体说明一个节点保持隔离的时间，而是给出了它在任何时候离开的可能性有多大。如果这对您没有意义，您应该阅读一下指数分布。

虽然我没有在这里使用它，但您可能想知道 EoN 附带的工具可以让您animate the simulation。

import networkx as nx
import matplotlib.pyplot as plt 
import EoN
from collections import defaultdict

a = 0.1
b = 0.01
y = 0.001
d = 0.001
to_isolation_rate = 0.05
from_isolation_rate = 1

# Simple contagions
# the below is based on an example of a SEIR disease (there is an exposed state before becoming infectious)
# from https://arxiv.org/pdf/2001.02436.pdf

Gnp = nx.gnp_random_graph(500, 0.005)

H = nx.DiGraph() #For the spontaneous transitions
H.add_edge('I', 'R', rate = b)  # an infected node can be recovered/removed
H.add_edge('I', 'S', rate = y)  # an infected node can become susceptible again
H.add_edge('R', 'S', rate = d)  # a recovered node can become suscepticle again
H.add_edge('S', 'X', rate = to_isolation_rate)
H.add_edge('X', 'S', rate = from_isolation_rate)

J = nx.DiGraph() #for the induced transitions
J.add_edge(('I', 'S'),('I', 'I'), rate = a)  # a susceptible node can become infected from a neighbour
IC = defaultdict(lambda: 'S')

IC[0] = 'I'

return_statuses = ('S', 'I', 'R', 'X')
print('doing Gillespie simulation')

t, S, I, R, X = EoN.Gillespie_simple_contagion(Gnp, H, J, IC, return_statuses, tmax = 500)

print('done with simulation, now plotting')
plt.plot(t, S, label = 'Susceptible')
plt.plot(t, I, label = 'Infected')
plt.plot(t, R, label = 'Recovered')
plt.plot(t, X, label = 'Isolating')
plt.xlabel('$t$')
plt.ylabel('Number of nodes')
plt.legend()
plt.show()