<p>To do this we will introduce a new node status <code>'X'</code> to denote the susceptible isolating node.  I've added an edge from <code>'S'</code> to <code>'X'</code> and another from <code>'X'</code> to <code>'S'</code> in <code>'H'</code>.  I've added <code>'X'</code> to the <code>return_statuses</code>, and I've plotted it.</p>
<p>To simplify your code, I have removed the loop that assigned everything to initially be <code>'S'</code>.  Because it uses a <code>default_dict</code> it all starts out as <code>'S'</code> implicitly.</p>
<p>You mention the next &quot;time step&quot; and shifting with probability <code>p</code>.  This isn't what will happen in this code.  The code works in continuous time.  So it will need a rate.  It doesn't specify exactly how long a node remains isolated, but rather gives an idea of how likely it is to leave at any moment.  If this doesn't make sense to you, you should read up a little bit on the exponential distribution.</p>
<p>Although I haven't used it here, you may want to be aware that EoN comes with tools that allow you to <a href="https://epidemicsonnetworks.readthedocs.io/en/latest/Examples.html#visualizing-or-animating-disease-spread" rel="nofollow noreferrer">animate the simulation</a>.</p>
<pre><code>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()
</code></pre>
SIRS model with Agents on NetworkX
<p>I want to simulate and visualize a SIRS model on a NetworkX graph (where a recovered/removed node can again become susceptible). Each node also operates as an agent, and can at each time step, choose to isolate with a probability p, and thus be unable to become infected in that timestep.</p>
<p>I have simulated a SIRS model with EoN.Gillespie_simple_contagion and I am trying to understand if I can modify these methods to include the agent behavior, or if I need to write a bespoke method.</p>
<p>Ive seen its possible to modify some EoN methods: <a href="https://stackoverflow.com/questions/56927079/modified-sir-model">Modified SIR model</a></p>
<p>Here is the code I am using:</p>
<pre><code>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()
</code></pre>