Complex vs Simple Gillespie in Non-pharmaceutical intervention situations

[ConstantinosM]

Thanks a lot for this very useful module!

I have started working with it and I was wondering if you could possibly help me deciding which simulation function is the best for what I want to do. I have few questions also regarding visualizations.

I am trying to build a high fidelity model by using a SIR/SEIR model along with a custom model that generates beliefs over e.g. social distancing, mask-wearing etc. In other words each node will have a property (e.g. belief that masks are good, thus high probability for wearing a mask which results in lower transmission rate) that will get updated according to a behavioral model (e.g. opinion model, Bayesian model, etc) and it will use information of the neighboring nodes. Then the simulation will return the counts of S,I,R, and the number of masked/unmasked people.

I would like to be able to use node properties to handle when someone wears mask or not - e.g. if majority of surrounding nodes believe in masks). I would also like to use this external behavioral model to define when someone's belief changes at the same timestep that the Gillespie model runs the simulation.

From what I have understood, the most flexible function for running simulations on networks that your module has is the Gillespie_complex_contagion. Do you think that I need this function or the simple version of it? Shall I assign the 'Mask wearing' as a node property or status (such as S,I,R)?

Also I am not 100% sure about their differences (complex vs simple) especially these parts:

Initially intended for a complex contagion.  However, this can allow influence
from any nodes, not just immediate neighbors.

and

This does not handle complex contagions. It assumes that when an individual changes status either he/she has received a "transmission" from a single neighbor or he/she is changing status independently of any neighbors

I see from the examples also that you are adding additional networks (specifying self transitions, or node-to-node transitions along with custom functions when you are using the simple version but not when you are using the complex version)

If you could elaborate more on how actually you calculate transmission in each function would be great (if you need to use math to show me how a transmission is estimated in the two different cases please do!). I read a notebook with another example of yours in DynamicProcessesOnNetworks repo but still it is not clear to me.

It seems to me that the problem I am dealing with is like having 2 "diseases" (or 2 networks), one is the opinion and how it propagates and changes the mask wearing probability and the other is the actual infection which is affected by the transmission rate of the nodes. I have already looked at your 2 diseases example but I am not sure if this is the best approach. Mainly I am thinking to use either the simple or complex Gillespie function with some custom functions. Any help would be great!

For graphics: I tried to use **nx_kwargs (e.g. node_size, edgecolors='black') and although the node size changes nodes appear always without outline (strangely as from the code it seems that the animation function draws nodes and edges separately). I would like to have different colors or even alphas (for example indicating the probability that a node is believing in mask wearing).

I have seen the effort that you have put in the module and documentation and in any case that this leads to any publication or presentation I will cite your work as I really appreciate it.

[joelmiller]

Hi -

I think the best way to do this would be to have both disease and mask wearing be considered as a "status". So we'd have statuses like 'SM', 'IM', 'RM', for the nodes wearing masks and 'SX', 'IX', 'RX' for those without.

Then your mask wearing is related to what the majority of people around the individual are doing. So people would go from the non-wearing to wearing status (or back) based on the neighbors (and perhaps based on their own properties).

For that, we'd need to use the complex contagion. You'd have to write a custom function which would be able to calculate the rate at which the node is changing its status.

If the connections that transmit disease and the connections that affect behavior are different you'd probably have a single network, but put some attribute on the edges to denote which kind of contact it is. Then that custom function you're writing would need to check edge attributes.

Let me know when/if you run into trouble on this. I'm a bit busy right now (what infectious disease modeler isn't?), but I can hopefully help.

[joelmiller]

I'll try to come back to this in a day or so with an example of a complex contagion. I thought I'd created one, but I don't see it in the documentation.

[ConstantinosM]

Thank you so much for the prompt response!Out of coincidence I ended up coming up with your suggestion as well. I assumed that in the SIR model won't make any effect if the R is masked or unmasked.

As a first stage I would like to do something as you suggest (as you see with the code below). Ultimately though I would like to have an external model that will take the neighboring nodes' beliefs, do some processing, and output the new belief of the current node. Then if that belief is lets say above 0.7 the node becomes masked (MS,MI) else it is unmasked (US,UI). The whole simulation loop though is happening in the Gillespie complex function so I am not sure what's the best place to do that update.

I am still going back and forth between the complex and the simple approach as I am trying to understand which one is the best approach. The reason for this is the following (I use as basis one simple example that you have commented in the simulation.py function):

def rate_function(G, node, status, parameters):    
    #This function needs to return the rate at which node changes status.
    #
    tau,gamma = parameters                       
    if status[node] == 'UI') or (status[node] == 'MI'):
        return gamma  # Recovery rate not affected by masked neighbors      
    elif status[node] == 'US': 
      # we discount the transmission rate by the amount of Masked Neighbors/Total Neighbors   or we can use a probability attribute and do weighted average.                                      
        return tau*(len([nbr for nbr in G.neighbors(node) if status[nbr] == 'MI'])/len([nbr for nbr in G.neighbors(node)]))
    elif status[node] == 'MS': 
      # we discount the transmission rate as above but we add an extra 0.7 discount. So if none is umasked still you get a discount                                         
        return 0.7*tau*(len([nbr for nbr in G.neighbors(node) if status[nbr] == 'MI'])/len([nbr for nbr in G.neighbors(node)]))

    else:
        return 0

But in the transition choice I am not sure how to handle it. For example in the simple version I can specify a digraph and then write a function only for that transition (but I have no access to the status[node] as in the complex case, at least from what I know). In the complex example the transition_choice function would need the event or something I guess to indicate which transition to select:

def transition_choice(G, node, status, parameters):
    #this function needs to return the new status of node.  We already
    #know it is changing status.
    #
    #this function could be more elaborate if there were different
    #possible transitions that could happen.
    if (status[node] == 'UI') or (status[node] == 'MI'):                       
        return 'R'
    elif status[node] == 'US':
        return 'UI'
    elif status[node] == 'MS':
        return 'MI'

# Assign mask belief attribute:
masks = {node: np.random.random_sample() for node in G} # better use a beta distr
nx.set_node_attributes(G, values=masks, name = 'masked')

# Rest of your code as is

How can I add transitions when the status is the same? For example: UI (unmasked infected) changes to MI (masked infected) or R (Recovered). Also in which of the 3 functions, that the complex contagion uses, is the best in terms of simulation timing in order to update node/edge attributes as you suggest. I guess they have a specific temporal order that are being called.

what do you exactly mean with "put some attribute on the edges to denote which kind of contact it is"? Shall I use an attribute such as contact type with two possible values: opinion, infection? I am not sure how the simulation will iterate over different edge attributes. Which function should specify for example when an US (unmasked susceptible) converts to a MS (masked susceptible) or UI (unmasked infected)?

Any example with a complex contagion would be great!

[ConstantinosM]

So this is what I ended up doing assuming the following transitions: and the code:

N = 100
tmax = 30
init_infected = 10
beta = 1.5
gamma = 1.0
alpha = 0.4 # how much a node is affected by its neighbors prob of masking
mu = 1 #  scaling down masks
mask_sc = 0.1

G = nx.grid_2d_graph(N,N)

a=0.0
b=0.2
masks = {node: (b - a) * np.random.random_sample() + a for node in G}
nx.set_node_attributes(G, values = masks, name = 'masked')

print('Initial masked probability:', np.mean([G.nodes[nbr]['masked'] for nbr in G]))

# The idea is that if rate specifies when something happens then I deal only with one transition and not multiple possible ones
# and leave the algorithm to decide when something should happen.
# Transitions:
# US --> UI --> R
# |      |
# MS --> MI --> R

belief_list = []

def mask_weighting(G, source, **kwargs): # US or UI is converted to MS or MI
    global belief_list
    n_neighbors = len([nbr for nbr in G.neighbors(source)])
    mask_neigbors_prob =np.sum([G.nodes[nbr]['masked'] for nbr in G.neighbors(source)])/n_neighbors
    new_belief = G.nodes[source]['masked'] + kwargs['alpha']*mask_neigbors_prob
    G.nodes[source]['masked'] = np.minimum(1.0,new_belief) # new belief cant go beyond 1.0
    scale = G.nodes[source]['masked']
    return scale # The lower the mask prob the lowest the prob for an unmasked person to become masked

H = nx.DiGraph()  #DiGraph showing possible transitions that don't require an interaction
H.add_edge('UI', 'R', rate = gamma)   #I->R
H.add_edge('MI', 'R', rate = gamma)   #I->R
H.add_edge('UI', 'MI', rate = mu, rate_function=mask_weighting)   #UI->MI
H.add_edge('US', 'MS', rate = mu, rate_function=mask_weighting)   #US->MS

J = nx.DiGraph()    #DiGraph showing transition that does require an interaction.
J.add_edge(('UI', 'US'), ('UI', 'UI'), rate = beta) # unmasked infected meets unmasked susceptible
J.add_edge(('UI', 'MS'), ('UI', 'MI'), rate = mask_sc*beta) # unmasked infected meets masked susceptible
J.add_edge(('MI', 'MS'), ('MI', 'MI'), rate = (mask_sc/2)*beta) # small rate of a masked infected to infect a masked susceptible susceptible
J.add_edge(('MI', 'US'), ('MI', 'MI'), rate = mask_sc*beta) # higher when susceptible doesn't wear a mask

IC = defaultdict(lambda: 'US')
for node in init_infected:
    IC[node] = 'UI'

return_statuses = ['US', 'UI', 'R', 'MS', 'MI']

# ANIMATION
color_dict = {'US': '#009a80', 'UI':'#ff2000', 'R':'gray', 'MS': '#5AB3E6', 'MI':'#a15f5a'}#,'Vac': '#5AB3E6'}
pos = {node:node for node in G}
tex = False
sim_kwargs = {'color_dict':color_dict, 'pos':pos}

sim = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses, tmax = tmax,
                            return_full_data=True, sim_kwargs=sim_kwargs,
                            spont_kwargs = {'alpha':alpha})

times, D = sim.summary()

newD = {'US': D['US'], 'UI': D['UI'], 'R': D['R'], 'MS': D['MS'], 'MI': D['MI'], 'Total Inf':D['UI']+D['MI']}

new_timeseries = (times, newD)
sim.add_timeseries(new_timeseries, color_dict={'US':'#009a80', 'UI':'#ff2020', 'R':'gray',
                                               'MS': '#5AB3E6', 'MI':'#a15f5a', 'Total Inf':'#00c3ff'})

ani=sim.animate(ts_plots=[['UI','MI'], ['MS', 'US'], ['Total Inf']], node_size = 4)

And I got the effect of having masked people containing the virus spread:

Whenever you have some time, if you could help me convert this to the complex_contagion scenario with edge attributes would be really great!

[joelmiller]

Hi,

Sorry for the delay getting back to you. I'm a bit busy, but hoping to do an example out of what you're doing. However, in the mean-time, I just discovered that I did write up an example here: https://arxiv.org/pdf/2001.02436.pdf starting on page 17. Take a look and see if it helps you.

Otherwise, can you respond telling me what you'd like the rule for masks to be, and it might take time, but I can try to get a good example out to you.

[ConstantinosM]

It is fine Joel! I really appreciate your help. Yes please if you could get an example with non trivial state transitions would be great. I have seen the example you mentioned and another example that you have with complex_contagion but they do not account for complex transitions as my example above (e.g. from state A you can get to state B or C with some probability).

As a start I would like to have a rule similar to the one I have in my example. So someone will update his current probability of wearing a mask (node attribute) according to the weighted average of probabilities of wearing a mask of his close neighbors. I know that this is a monotonic function that eventually will convert everyone into a mask wearer but for now it is fine. I can try different rules later if I have an initial model that can work on. In my example I initialized the mask probability attribute from a uniform and this also can be changed as well to reflect various other characteristics (e.g. education, political affiliation etc)

Thank you again so much!