为此，我们将引入一个新的节点状态 'X' 来表示易受影响的隔离节点。我在 'S' 中添加了从 'X' 到 'X' 的边以及从 'S' 到 'H' 的另一个边。我已将 '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',rate = d)  # a recovered node can become suscepticle again
H.add_edge('S','X',rate = to_isolation_rate)
H.add_edge('X',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','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,label = 'Susceptible')
plt.plot(t,label = 'Infected')
plt.plot(t,label = 'Recovered')
plt.plot(t,X,label = 'Isolating')
plt.xlabel('$t$')
plt.ylabel('Number of nodes')
plt.legend()
plt.show()

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

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,R