Many contagion processes evolving on populations do so simultaneously, interacting over time. Examples are co-evolution of social processes and diseases, such as increased hand washing during the influenza (flu) season and mask wearing during pandemic outbreaks, e.g., [9]; and co-transmission of information and evacuation decisions during natural disasters [29]. Commensurately, multi-contagion agent-based simulations (ABSs) that represent populations as networks in order to capture interactions between pairs of nodes (e.g., people are represented as nodes and interactions as edges in a network) are becoming more popular.
Currently, several simulation systems can compute interacting contagion dynamics, but they do so by incorporating data that must be preprocessed for a particular simulation, or additional software must be written. For example, consider each of two interacting contagions, c 1 and c 2 , spreading on separate networks, G 1 and G 2 , respectively. Generally, there are common nodes and/or edges in the graphs so that there is interaction between contagions. Simulators are typically set up to store one graph, so that they require that a union graph G be formed (G = G 1 ∪ G 2 ) where labels on edges designate whether an edge is in G 1 , G 2 , or both. Graph nodes are handled analogously. Consequently, if there are three graphs and simulations need to be performed with all combinations of two graphs, then three union graphs need be generated. With n g of distinct graphs and n inst graphs used in each simulation, the number of required union graphs is the binomial coefficient C(n g , n inst ), which can be large.
Also, software must be added to frameworks to manage contagions. As an example for discrete time simulation, one might execute one contagion at each time step, so that if there are n c contagions and t max simulation time steps (e.g., days), a simulator might execute n c × t max time steps. But this approach adversely affects down-stream post-processing because "time" must be manipulated accordingly. For other frameworks, code has to be written to coordinate the execution of contagions at each time (see Section 1.3). Our system obviates the need for these changes.
We present a new ABS framework for modeling the spread of any number of contagions on any number of networks, using a structured and principled approach in both the software system and its use in applications. We also provide a formal model for our simulation system and case studies to illustrate its use in applications.
1. Multi-Contagion Graph Dynamical Systems Formalism. We start with a graph dynamical systems (GDS) formalism [1] for single-contagion discrete dynamical systems. We extend that formalism to a multiple contagion GDS, MCGDS, operating over multiple graph instances. Any number n c ≥ 1 of contagions can operate on any number n g of graphs, where 1 ≤ n g ≤ n c , because while each contagion spreads on one graph, two or more contagions may spread on a single graph. The formal model is provided in Section 2.
The software framework, called CSonNet (Contagion Simulation ON NETworks) is written in Python. Our simulation system implements concurrency through multi-processing. The system can simulate any number n c ≥ 1 of contagions operating over any number n g of graphs, where 1 ≤ n g ≤ n c . The graph instances, in general, can be disjoint, the same graph, or have common subsets of nodes and edges. Contagions can evolve independently or can interact as they evolve. Contagion models (e.g., SIR epidemic models, various threshold models for social contagions) are added to the system in one of two ways: (i) Python code that conforms to a specified interface is written and added to the system, or (ii) a user composes models through input files, using a set of rules provided by the system. Option (i) is the only case where new code is written; all other inputs (e.g., specification of graphs, the mapping of contagions to the networks on which they evolve, rule-based models) are given through input files for a simulation. The system is overviewed in Section 3 and its strong scaling performance is evaluated in Section 4.
3. Case Studies. We provide two case studies, each with two non-interacting and two interacting contagion simulations: (i) SIR 1 -SIR 2 epidemic contagions, and (ii) Threshold-SIR mixed social and epidemic contagions. Social networks with up to 75,877 nodes are used. While the system can compute contagions on much larger graphs, we choose graphs where fairly demanding multi-contagion simulations can be completed in less than 300 seconds of computational time, since in practice parametric simulation studies over large parameter spaces are typically required. These illustrate several features of CSonNet (see Section 5).
This simulation system significantly extends the single-contagion system described in [23]; that system does not implement the MCGDS of Contribution 1. Our implementation here (Contribution 2) is a large extension of the previous system in terms of capabilities (e.g., software added and software modified). None of the case studies here (Contribution 3) could be executed with the system in [23] without the workarounds specified in Section 1.1.
Multiple contagion simulation systems. To the best of our knowledge, no simulator has multi-contagion capabilities "out of the box" as is the case for our framework. However, because several frameworks enable models to be added to them (as does ours), it is possible to use these simulators for multiple-contagion scenarios, as described in Section 1.1. But these workarounds come with the price of requiring additional code and/or requiring more preprocessing of graphs before running simulations. Among the frameworks for simulations on networks that fit this description, across a range of capabilities, are: NetLogo [24], NDlib [25], MASON [15], and Repast HPC [10].
Studies of multiple contagion systems. Simulations of information spreading and evacuation decision-making in the context of hurricane evacuation modeling are discussed in [29]. Simulations of agents with limited attention spans for which multiple contagions (e.g., ideas) must compete are discussed in [28]. Agent-based models for competing languages are summarized in [21].
A number of papers have addressed models for competing and/or cooperating contagions (such as epidemics, product information, and misinformation in social media) over networks (see e.g., [17,19,28]). Polarization and consensus are studied with competing contagions [27]. One contagion that overtakes and halts a second contagion is studied in [7,20]. A model for a simple contagion producing a complex contagion is given in [17]. Game theory is used to analyze competing contagions in [11]. Many publications have studied optimization problems (e.g., seeding methods for influence maximization, minimizing the spread of contagions) in the context of multiple contagions (see e.g., [6,7]).
The underlying model implemented in our simulation system CSonNet is called a graph dynamical system (GDS). The formalism [1,18] addresses a single contagion, but naturally extends to multiple contagions, which is done here. A GDS can simulate Turing machines for specific complexity classes [1,4,5].
, where G j is the graph on which the contagion c j propagates, V j and E j denote its vertex (i.e., node) and edge sets, 1 ≤ j ≤ n c ; (ii) F c : a sequence of sequences F j , where each F j is the sequence of local functions for contagion c j ; (iii) K c : a sequence of vertex (node) state sets K j , where each K j is the set of admissible vertex states for contagion c j ; and (iv) R c : a sequence of specifications R j , where each R j is the order in which local functions are executed for c j . Each of these elements is defined below.
The graph G j (V j , E j ) represents the interaction network for c j . Let
Agents are vertices and pairwise agent interactions are edges in G j . In general, edges may be directed or undirected. In this paper, for simplicity, we will consider graphs with undirected edges. Vertex (resp., edge) sets across contagions may be the same, disjoint, or have some common subset of nodes (resp., edges).
Each agent v i ∈ V 1 ∪ . . . ∪ V nc has an agent state or vertex state s i which is a sequence s i = (s i,c1 , s i,c2 , . . . , s i,cn c ) of n c elements. Each s i,cj ∈ K j . The system state (also called a configuration) s = (s 1 , s 2 , . . . , s n ) is the sequence of n vertex states. Let the local states, denoted by s[v i ], represent the sequence of states of vertex v i and all of its distance-1 neighbors in each
, where s cj [v i ] is the sequence of length d j (v i ) + 1 of the states of v i and each of its distance-1 neighbors in G j for contagion c j , and d j (v i ) is the degree of v i in G j . These quantities, at time t, are denoted by s t i , s t i,cj , s t , and s t [v i ] respectively. Each agent v i has a local function f i,cj ∈ F j that determines its state transitions for contagion c j . The state of agent v i at time t + 1 for c j is given by
Thus, the next state of v i , with respect to each contagion c j is a function of the current state of v i and those of its distance-1 neighbors in every G j over all contagions. That is, the argument on the right hand side of Equation ( 1) is the same for each f i,cj , for a fixed v i . This expression-and specifically the local functions for each contagion-indicates explicitly how each contagion affects (i.e., interacts with) every other contagion. We assume that specification of update order R j for the local functions for c j corresponds to the synchronous update scheme. That is, all agents v i evaluate their local functions f i,cj , 1 ≤ i ≤ n, for c j and update their states s t+1 i,cj simultaneously at each time step. Furthermore, we assume the updates across contagions are also done in parallel, i.e., all R j are parallel with each other. Hence, f i,cj are computed in parallel for all 1 ≤ i ≤ n and for all 1 ≤ j ≤ n c in Equation ( 1). This enables greater parallelization of simulation computations, leading to more efficient calculations. However, other update orders R j can be used. One example is a sequence W of the node IDs in G j (of length n j ) and the local function for each node is computed in the order specified by W .
Two contagions propagate on the single network at the left in Figure 1, where nodes are people and edges are interactions. The two contagions are: a social contagion c 1 on mask wearing and a disease contagion c 2 . Contagion c 1 uses a threshold model [12], with K 1 = B = {0, 1}, where state 0 (resp., state 1) means that a person (node) is not (resp., is) wearing a mask. The local function f i,c1 for all i ∈ {1, 2, 3, 4} and c 1 is as follows. A node v i changes from state 0 to 1 if at least a threshold θ = 1 number of its neighbors are in state 1. Contagion c 2 is a disease model with states K 2 = {S, I, R}, where the states are susceptible (S), infectious (I), and recovered (R). The local function f i,c2 for all i ∈ {1, 2, 3, 4} and c 2 is as follows. If v i is in state S, then v i changes to state I with probability p = p base • m I • m S for each neighbor in state I, where p base = 0.008, m I = 0.2 if the infected neighbor is wearing a mask and m I = 1 otherwise, and m S = 0.2 if v i is wearing a mask and m S = 1 otherwise. If v i is in state I and transitioned to state I at time t I , then it transitions to state R at time t = t I + t inf , where t inf = 3 is node v i 's infectious duration. If v i is in state R, it remains in state R. Contagion c 1 affects c 2 in that people that are wearing masks have a lesser probability of contracting and transmitting the disease.
The system states s at t = 0, 1, and 2 are given in Figure 1. Contagion c 1 states are under "Social" and c 2 states are under "SIR." At t = 1, c 1 spreads to each of v 2 and v 4 because these latter two nodes have v 1 as a neighbor, s 0 1,1 = 1, and the threshold for all nodes is θ = 1. So at the end of t = 1, three of the four nodes are wearing masks. Also at t = 1, for contagion c 2 , v 3 is initially infected, but the Bernoulli trials for nodes v 2 and v 4 do not result in contagion spread, so the states remain s 1 2,2 = s 1 4,2 = S. (Note that at t = 1 for c 2 and node v 3 , t < t I + t inf = 0 + 3, so f 3,c2 returns the next state as s 3,c2 = I, which is the current state. The local function for v 3 will return the next state as I until t = 3, at which time s 3,c2 = R.) At t = 2 and for c 1 , v 3 changes to s 2 3,1 = 1 due to simple contagion spread from v 2 and v 4 . For c 2 at t = 2, the edge probability p = p base •m I •m S = (0.008)(1)(0.2) = 0.0016 for the edge from v 3 to v 2 and from v 3 to v 4 because v 2 and v 4 are wearing masks at the end of t = 1 but v 3 is not. The random number in [0, 1] generated for the first edge is 0.0013 < p = 0.0016 (i.e., the Bernoulli trial is successful) and so v 3 infects v 2 and s 2 2,2 = I. However, the Bernoulli trial for v 3 to infect v 4 is not successful because the drawn random number is 0.732 > p, so s 2 4,2 = S and therefore the states at t = 2 are as shown in Figure 1. A case study using this two-contagion model is given in Section 5.1. 3 CSonNet Modeling and Simulation Software System
CSonNet is a discrete-time multi-contagion ABS framework for networked populations. CSonNet implements the MCGDS model of Section 2 in the form of simulations. A simulation begins with reading in from file: (i) all graph instances G c ; (ii) all local functions F c over all nodes v i , 1 ≤ i ≤ n, and all contagions c j , 1 ≤ j ≤ n c ; (iii) the mapping M j, of contagion c j to graph G ; (iv) the initial state assignments I to all v i for each c j ; (v) the number n i of iterations (i.e., simulation instances) to run; (vi) the maximum number t max of time steps to run per iteration; and (vii) the number n wp of worker processes that perform the computations. The number n g of graphs and number n c of contagions are determined from G c and F c , respectively. In particular, the number and types of contagions are completely specified by the local functions f i,cj of Equation ( 1) in the sequence F c and K c . The local functions are not entered as equations in input files, but rather as models M st i,j , that implement these local functions for v i and c j . Section 3.2 below provides an example.
After reading all inputs, the main process of the simulation instantiates n wp worker processes that carry out iterations in parallel on multi-core hardware computing nodes. The main process provides to each worker process the appropriate graphs, local functions, mapping of contagion to graph, initial conditions for particular iterations, and other parameters that the worker process needs, and then starts the worker processes. An iteration consists of computing the dynamics based on the states of all nodes at time t = 0 and over all contagions, up through time t max . Specifically, an iteration consists of computations over the following nested loops: (i) over all time steps t ∈ [0, t max -1]; (ii) for each time, over all contagions c j , j ∈ [1, n c ]; and (iii) for each contagion, over all nodes v i , 1 ≤ i ≤ n. For each combination of (t, c j , v i ), Equation ( 1) is evaluated, generating v i 's next state s t+1 i,cj . If the next state is different from the current state of v i , then this next state is written to file, along with the corresponding iteration number, t, c j , and v i ; this is the simulation output.
State transition models M st i,j , which are specified in input files, represent the local functions of Section 2. Below are examples of two state transition model files; each specifies the contagion number, followed by a row of entities that constitute the elements of a rule; all subsequent lines contain particular rule names and parameter values for these rules. Hence, a model is composed of rules. Currently, there is a fixed but extensible set of rules. Moreover, the model files below are used to perform a simulation such as that in Section 2.2.
Contagion 1 (c 1 ) has one rule. It is a (deterministic) threshold model that applies to all nodes, and describes the transition from state 0 to state 1 when at least 3 neighbors of a node are in state 1 (the cause). The threshold change based on influence from the SIR model is zero, meaning that the SIR contagion does not affect the threshold-based contagion. The last value of 1 represents a minimum threshold for the nodes, in the event of a threshold decrease; the threshold cannot go below 1.
Contagion 2 (c 2 ) is an SIR model. There are two rules, one for the state transition S → I and one for the transition I → R. The transition S → I is governed by an edge probability of 0.006. If the infected (resp., susceptible) node is in state 1 for c 1 , then probability is reduced by the factor 0.5 (resp., 0.5). The Python list "[1,I]" indicates that states 1 from c 1 and state I from c 2 influence the transition S → I. Thus, c 1 influences c 2 . The second rule states that a node stays in the infected state for 10 time units before transitioning I → R.
The networks in Section 4.1 are used to perform strong scaling studies of simulation time in Section 4.2.
Networks used in performance analyses and the case studies of Section 5 are given in Table 1. Two are face-to-face human social contact networks and the third (Epinions) is an online social network.
Figure 2 shows strong scaling results for two-contagion simulations for each of the three networks. This is the total execution time over all worker processes (from the start of the first worker process to the end of the last executing process).
The city-based human contact networks (first two entries) were made with the procedures in [3]. Each network is the giant component of the network, since we run dynamics on these networks. Property computations were performed with [2]. The two contagions do not interact in Figure 2a; there is interaction between contagions in Figure 2b. The times are greater in the latter plot because all neighbors for both contagions must be iterated over to update the state of each contagion; these neighborhood iterations take place at each time step of each iteration for one simulation. Each data point, with ± one standard deviation error bars, represents the results of ten simulations for each set of conditions. The data in both plots indicate that CSonNet exhibits strong scaling for independent and interacting contagions. The simulation conditions are appreciably onerous. Each time data point is for a simulation with 100 iterations, each over 100 time steps. In practice (e.g., for epidemic simulations), less than 100 iterations are performed (often on the order of 30), and typically about 14-30 time steps (days) are simulated. Yet, computations by worker processes complete in under 300 seconds for n wp = 32 in Figure 2b. We are interested in determining the sizes of networks on which simulations can be run in under five minutes because typical studies involve large parametric studies with many simulations.
Two multi-contagion case studies are presented. The first case study is a Threshold-SIR system spreading on the Newport News network, and is motivated in part by [8,26]. The second is an SIR 1 -SIR 2 system spreading on the Danville network.
It is inspired by [16,22]. In the first case study, one contagion inhibits the other; in the second study, two contagions reinforce each other. The purpose of these case studies is to demonstrate multi-contagion capabilities of the code.
One contagion model (for c 1 ) is a Granovetter threshold model [12], with threshold θ = 3, and is used to model mask wearing during COVID. States 0 and 1 indicate that a person is not wearing a mask and is wearing a mask respectively. A node transitions from state 0 to state 1 if at least θ of its neighbors are in Consider a person v j in state I (i.e., infected with COVID) and another person who is v i in state S (i.e., is susceptible). Each of v i and v j may or may not be wearing a mask during an encounter. When a person wears a mask, the probability of transmission is reduced by a multiplicative factor attributed to each person. For each person v i wearing a mask, we let m i = 0.5; otherwise, the factor is m i = 1. Thus, the edge probability p e,i for transmission along the edge e between a susceptible person v i and an infected person v j is given by p e,i = w e,base,i • m i • m j . In our case study, w e,base,i = 0.006 and the infectious duration t inf = 10 days. COVID does not affect mask wearing in our model. The state transition model input files for this case study-for the interacting contagions simulation-were provided in Section 3.2. Figure 3 provides cumulative infection curves (for c 2 ) and cumulative numbers of nodes in state 1 (for c 1 ) for the Newport News network. Seed nodes for each contagion are chosen uniformly at random and separately for each contagion: 0.0062 fraction of nodes for c 1 and 0.0062 fraction of nodes for c 2 in each iteration. The particular seed nodes are different for each of the 100 iterations, but the same seeds are used across the two simulations, for comparison. Figure 3a provides baseline data for non-interacting contagions. Figure 3b shows results for the case of mask wearing (c 1 ) affecting the spread of COVID (c 2 ). The latter plot indicates that mask wearing reduces the spread of COVID [14].
In this SIR 1 -SIR 2 system, if contagion c 1 affects contagion c 2 , then this is realized by increasing the probability that a node v contracts c 2 , given that it has already contracted c 1 . More formally, consider a two-contagion system where the contagion models are SIR i and SIR j with i, j ∈ {1, 2} and i = j. If the simultaneously evolving contagions do not interact, then m k = 1 for contagion SIR k (k = 1, 2) by definition. However, if contagion SIR i affects SIR j , and a node v has already contracted contagion c i , then the edge weight (i.e., probability of infection) for c j , is w e,j = w e,base,j • m j . If 0 ≤ m j < 1, then contracting contagion c i reduces the probability of contracting c j . If m j > 1, then contracting contagion c i increases the probability of contracting c j . For SIR 1 , the base edge weight is w e,base,1 = 0.01, the infectious duration is t inf,1 = 10 days, and for interacting contagions m 1 = 8. For SIR 2 , w e,base,2 = 0.005, t inf,2 = 12 days, and for interacting contagions m 2 = 5. Since m 1 , m 2 > 1, the contagions, when interacting, reinforce each other. Figure 4 shows results from three simulations. Seed nodes for each contagion are chosen uniformly at random: 0.0077 fraction of nodes for c 1 and 0.0077 fraction of nodes for c 2 in each iteration. The particular seed nodes are different for each iteration, but the same seeds are used across the three simulations, for comparison. All plots are cumulative fractions of infected individuals as a function of time. Error bars at each time are ± one standard deviations over the 100 iterations per simulation. In Figure 4a, the two contagions do not interact. In Figure 4b, c 1 affects c 2 , but not vice versa; so the cumulative fraction of infections increases for c 2 while remaining unchanged for c 1 . In Figure 4c, the two contagions affect each other and the cumulative infected fractions increase for both contagions.