Competition, Collaboration, and Optimization in Multiple Interacting Spreading Processes

TL;DR

This paper introduces a novel dynamic message-passing framework for modeling and controlling two interacting spreading processes on networks, enabling optimal resource allocation for either promoting or limiting spread.

Contribution

It derives exact equations for two interacting processes on trees, extends them to general graphs, and develops an optimization framework for controlling spread under constraints.

Findings

Accurate modeling of interacting processes on networks.

Effective algorithms for maximizing or limiting spread.

Validated on synthetic and real-world networks.

Abstract

Competition and collaboration are at the heart of multi-agent probabilistic spreading processes. The battle on public opinion and competitive marketing campaigns are typical examples of the former, while the joint spread of multiple diseases such as HIV and tuberculosis demonstrates the latter. These spreads are influenced by the underlying network topology, the infection rates between network constituents, recovery rates and, equally importantly, the interactions between the spreading processes themselves. Here, for the first time we derive dynamic message-passing equations that provide an exact description of the dynamics of two interacting spreading processes on tree graphs, and develop systematic low-complexity models that predict the spread on general graphs. We also develop a theoretical framework for an optimal control of interacting spreading processes through an optimized resource allocation under budget constraints and within a finite time window. Derived algorithms can be used to maximize the desired spread in the presence of a rival competitive process, and to limit the spread through vaccination in the case of coupled infectious diseases. We demonstrate the efficacy of the framework and optimization method on both synthetic and real-world networks.

Figures

Tables

Table 1. Table 3: Comparing different deployment methods for a competitive scenario on various networks. On each network, we randomly choose 0.05 ​ N 0.05 𝑁 0.05N nodes as seeds of process B 𝐵 B and the same total budget for the competitive process A 𝐴 A , where N 𝑁 N is number of nodes in the network. Infection parameters are α A = 0.2 subscript 𝛼 𝐴 0.2 \alpha_{A}\!=\!0.2 and α B = 0.3 subscript 𝛼 𝐵 0.3 \alpha_{B}\!=\!0.3 . The initial parameter ϵ italic-ϵ \epsilon which forces the budget limits per node as in ( 213 ) is initially set to 0.1 but decays exponentially with the iteration steps. The optimization procedure is iterated 10 times and the best result is selected. The objective is to minimize the spreading of process B 𝐵 B ; the normalized total spreading ∑ i P B i ​ ( T ) / N subscript 𝑖 superscript subscript 𝑃 𝐵 𝑖 𝑇 𝑁 {\sum_{i}P_{B}^{i}(T)}/{N} at time T = 3 𝑇 3 T\!=\!3 is shown. In all methods budgets are available at time zero. Uniform Allocation assumes all the budgets are uniformly allocated to all nodes at time zero (except those infected with process B 𝐵 B ). In Free Spreading , no agents are allocated for process A 𝐴 A . In Blocking, competitive agents are not infectious ( p = 0 𝑝 0 p\!=\!0 ), all the budget is used to contain the spreading of B 𝐵 B . Best results are denoted by bold fonts.

Network	No. of nodes	DMP	Uniform	K-shell	HDA	Blocking	Free Spreading
Football [25]	115	0.5829	0.5674	0.8264	0.7731	0.7860	0.8264
Lesmis [47]	77	0.0834	0.2211	0.1962	0.0990	0.1016	0.3549
Karate [94]	35	0.4206	0.4902	0.5470	0.4212	0.4966	0.5472
Power [89]	4941	0.0500	0.0621	0.0641	0.0632	0.0501	0.0644
Polbooks [48]	105	0.2243	0.3733	0.3102	0.2779	0.4316	0.5744

Table 2. Table 4: DMP based resource deployment methods for the Football network under the same condition as in Fig.  6 . The objective is to minimize the spreading of B 𝐵 B , i.e. \textcolor black, minimize the fraction of B 𝐵 B nodes at time T = 5 𝑇 5 T\!=\!5 , ∑ i P B i ​ ( T ) / N subscript 𝑖 superscript subscript 𝑃 𝐵 𝑖 𝑇 𝑁 {\sum_{i}P_{B}^{i}(T)}/{N} .

Time	DMP optimal	DMP greedy
$t=1$	0.0767	0.0556
$t=2$	0.2157	0.1871
$t=3$	0.3109	0.4336
$t=4$	0.3211	0.5615
$t=5$	0.3211	0.5628

Table 3. Table 5: Comparing different resource-allocation approaches for collaborative processes on various networks. On each network, we randomly select 0.05 ​ N 0.05 𝑁 0.05N nodes as seeds of process B 𝐵 B and the same total budget for process A 𝐴 A to be optimally allocated, where N 𝑁 N is total number of nodes in the network. The infection parameters chosen are α A = 0.2 subscript 𝛼 𝐴 0.2 \alpha_{A}\!=\!0.2 , α B = 0.3 subscript 𝛼 𝐵 0.3 \alpha_{B}\!=\!0.3 , α A ​ B = 0.4 subscript 𝛼 𝐴 𝐵 0.4 \alpha_{AB}\!=\!0.4 and α B ​ A = 0.5 subscript 𝛼 𝐵 𝐴 0.5 \alpha_{BA}\!=\!0.5 . The initial parameter ϵ italic-ϵ \epsilon which forces the budget limits per node as in ( 213 ) is initially set to 0.1 but decays exponentially with the iteration steps. The optimization procedure is iterated 10 times and the best result is selected. The objective function in this case is to maximize the spreading of processes A 𝐴 A and B 𝐵 B , i.e., minimize the fraction of susceptible nodes ∑ i P S i ​ ( T ) / N subscript 𝑖 superscript subscript 𝑃 𝑆 𝑖 𝑇 𝑁 {\sum_{i}P_{S}^{i}(T)}/{N} at time T = 3 𝑇 3 T\!=\!3 . In DMP allocation-based optimization, we assume the seeding budget for B 𝐵 B is available at time zero. Uniform allocation assumes that the budget for A 𝐴 A is allocated uniformly to all free nodes initially. HDA deployment  [ 76 ] and K-shell  [ 46 ] seeding are used as before; in Free Spreading, no collaborative spreading budget A 𝐴 A is allocated. Lowest fraction of susceptible nodes is denoted by bold fonts.

Network	No. of nodes	DMP Allocation	Uniform Allocation	K-shell	HDA	Free Spreading
Football [25]	115	0.0536	0.0582	0.1736	0.1543	0.1735
Lesmis [47]	77	0.2222	0.3832	0.3151	0.3051	0.6451
Karate [94]	35	0.2771	0.3743	0.4527	0.2481	0.4529
Power [89]	4941	0.7652	0.7930	0.8434	0.9024	0.9355
Polbooks [48]	105	0.1524	0.2065	0.3347	0.2204	0.4255

Table 4. Table 6: Different vaccine-allocation policies in a collaborative process on difference benchmark networks. On each network, we randomly choose 0.05 ​ N 0.05 𝑁 0.05N of the nodes as seeds for process B 𝐵 B and 0.01 ​ N 0.01 𝑁 0.01N nodes as seeds for process A 𝐴 A ; the total vaccination budget is 0.05 ​ N 0.05 𝑁 0.05N , where N 𝑁 N is number of nodes in the network. Infection parameters are arbitrarily set to α A = 0.4 subscript 𝛼 𝐴 0.4 \alpha_{A}\!=\!0.4 , α B = 0.4 subscript 𝛼 𝐵 0.4 \alpha_{B}\!=\!0.4 , α A ​ B = 0.9 subscript 𝛼 𝐴 𝐵 0.9 \alpha_{AB}\!=\!0.9 and α B ​ A = 0.9 subscript 𝛼 𝐵 𝐴 0.9 \alpha_{BA}\!=\!0.9 . The parameter ϵ italic-ϵ \epsilon in ( 213 ) is set to 0.01 0.01 0.01 . The optimization procedure is iterated 10 times and the best result is selected. The objective is to minimize the spreading of processes A 𝐴 A and B 𝐵 B , i.e., maximize the fraction of susceptible nodes ∑ i P S i ​ ( T ) / N subscript 𝑖 superscript subscript 𝑃 𝑆 𝑖 𝑇 𝑁 {\sum_{i}P_{S}^{i}(T)}/{N} at the end of the process (time T = 3 𝑇 3 T\!=\!3 ). In DMP-based optimal allocation, we assume that all budgets are available at time zero. Uniform Allocation assumes all budgets to be uniformly allocated to all nodes at time zero. In Free Spreading, no vaccine is allocated. A detailed description of the vaccine allocation problem can be found in \textcolor blackAppendix I.

Network	Number of nodes	DMP Allocation	Uniformed Allocation	Free Spreading
Football [25]	115	0.1887	0.0607	0.0294
Lesmis [47]	77	0.9455	0.5519	0.4941
Karate [94]	35	0.4676	0.3491	0.3123
Power [89]	4941	0.8891	0.8650	0.8356
Polbooks [48]	105	0.5483	0.3071	0.2354

Equations

$$P_S^i(t)=\sum _{{\tau }_i^A>t}\sum _{{\tau }_i^B>t}{m}_{T_A,T_B}^i({\tau }_i^A,{\tau }_i^B),$$

$$P_A^i(t)=\sum _{{\tau }_i^A\le t}\sum _{{\tau }_i^B}{m}_{T_A,T_B}^i({\tau }_i^A,{\tau }_i^B),$$

$$P_B^i(t)=\sum _{{\tau }_i^B\le t}\sum _{{\tau }_i^A}{m}_{T_A,T_B}^i({\tau }_i^A,{\tau }_i^B),$$

$$P_{AB}^i(t)=\sum _{{\tau }_i^A\le t}\sum _{{\tau }_i^B\le t}{m}_{T_A,T_B}^i({\tau }_i^A,{\tau }_i^B).$$

$$S(i)+A(j)$$

$$S(i)+B(j)$$

$$v_A^i(t)=1-\prod _{j\in \partial i}(1-{\alpha }_{ji}^A\text{mathbbm}1[{\sigma }_j(t)=A]),$$

$$v_B^i(t)=1-\prod _{j\in \partial i}(1-{\alpha }_{ji}^B\text{mathbbm}1[{\sigma }_j(t)=B]),$$

$$Z_i=1-v_A^i(t)v_B^i(t).$$

$$q_{S\to A}^i(t)=\frac{v_A^i(t)(1-v_B^i(t))}{Z_i(t)},$$

$$q_{S\to B}^i(t)=\frac{v_B^i(t)(1-v_A^i(t))}{Z_i(t)},$$

$$q_{S\to S}^i(t)=\frac{(1-v_A^i(t)-v_B^i(t)+v_A^i(t)v_B^i(t))}{Z_i(t)},$$

$$S(i)+A(j)$$

$$S(i)+B(j)$$

$$A(i)+B(j)$$

$$S(i)+AB(j)\stackrel{{\alpha }_{ji}^A\times {\alpha }_{ji}^B}{\to }AB(i)+AB(j)$$

$$P_S^i(t)=P_S^i(0)\prod _{k\in \partial i\backslash j}{\theta }_{A,B}^{k\to i}(t,t),$$

$${\theta }_{A,B}^{k\to i}(t,t)=\sum _{{\tau }_k^A}\sum _{{\tau }_k^B}\prod _{t^{\prime }=0}^{t-1}(1-{\alpha }_{ki}^A\text{mathbbm}1[{\tau }_k^A\le t^{\prime }])$$

$$\times \prod _{t^{\prime \prime }=0}^{t-1}(1-{\alpha }_{ki}^B\text{mathbbm}1[{\tau }_k^B\le t^{\prime \prime }])m_{T_A,T_B}^{k\to i}({\tau }_k^A,{\tau }_k^B).$$

$$P_A^i(t)=\sum _{t^{\prime }\le t}{\mu }_A^i(t^{\prime }),$$

$$P_B^i(t)=\sum _{t^{\prime }\le t}{\mu }_B^i(t^{\prime }),$$

$${\mu }_A^i(t)=\sum _{{\tau }_i^B}{m}_{t,t}^i(t,{\tau }_i^B)=P_A^i(t)-P_A^i(t-1),$$

$${\mu }_B^i(t)=\sum _{{\tau }_i^A}{m}_{t,t}^i({\tau }_i^A,t)=P_B^i(t)-P_B^i(t-1).$$

$$P_{AB}^i(t)=P_A^i(t)+P_B^i(t)+P_S^i(t)-1.$$

$$m_{t,t}^i({\tau }_i^A,{\tau }_i^B)={\mu }_A^i({\tau }_i^A){\mu }_B^i({\tau }_i^B)$$

$$\text{if }{\alpha }_{ij}^{AB}={\alpha }_{ij}^A\text{ and }{\alpha }_{ij}^{BA}={\alpha }_{ij}^B.$$

$$\mathcal{L}=\underset{\text{objective}}{\underbrace{\mathcal{O}}}+\underset{\text{constraints}}{\underbrace{\mathcal{B}+\mathcal{P}+\mathcal{I}+\mathcal{D}}}.$$

$$\mathcal{O}=\sum _i(1-P_i^B(T)),$$

$$B_{\nu }=\sum _i{\nu }^i(0),$$

$$\mathcal{B}={\lambda }^{Bu}(B_{\nu }-\sum _i{\nu }^i(0)).$$

$$\mathcal{P}=ϵ\sum _i(\log (\overline{\nu }-{\nu }^i(0))+\log ({\nu }^i(0)-\underline{\nu }))$$

$$m_{T_A+1,T_B+1}^i({\tau }_i^A,{\tau }_i^B)=\sum _{\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i}}{W}^i({\tau }_i^A,{\tau }_i^B;\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i})\prod _{k\in \partial i}{m}_{T_A,T_B}^{k\to i}({\tau }_k^A,{\tau }_k^B\Vert {\tau }_i^A,{\tau }_i^B),$$

$$m_{T_A+1,T_B+1}^{i\to j}({\tau }_i^A,{\tau }_i^B\Vert {\tau }_j^A,{\tau }_j^B)=\sum _{\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i\backslash j}}{W}^i({\tau }_i^A,{\tau }_i^B;\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i})\prod _{k\in \partial i\backslash j}{m}_{T_A,T_B}^{k\to i}({\tau }_k^A,{\tau }_k^B\Vert {\tau }_i^A,{\tau }_i^B).$$

$$m_{T_A,T_B}^i({\tau }_i^A,{\tau }_i^B)=m_{{\tau }_i^A,{\tau }_i^B}^i({\tau }_i^A,{\tau }_i^B)\forall T_A>{\tau }_i^A,\text{and }{T}_B>{\tau }_i^B.$$

$$m_{T_A+1,T_B+1}^i({\tau }_i^A,{\tau }_i^B)=\sum _{\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i}}{W}^i({\tau }_i^A,{\tau }_i^B;\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i})\prod _{k\in \partial i}{m}_{T_A,T_B}^{k\to i}({\tau }_k^A,{\tau }_k^B\Vert \ast ,\ast ),$$

$$m_{T_A+1,T_B+1}^{i\to j}({\tau }_i^A,{\tau }_i^B\Vert {\tau }_j^A,{\tau }_j^B)=\sum _{\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i\backslash j}}{W}^i({\tau }_i^A,{\tau }_i^B;\{{\tau }_k^A,{\tau }_k^B{\}}_{k\in \partial i})\prod _{k\in \partial i\backslash j}{m}_{T_A,T_B}^{k\to i}({\tau }_k^A,{\tau }_k^B\Vert \ast ,\ast ).$$

Full text

I Introduction

Spreading processes have become increasingly more important in the fast moving modern world, where physical mobility is cheaper and easier than ever before, and where information is passed instantaneously on multilayered interwoven webs of contacts. Consequently, pandemics, whether physical or virtual in the form of computer viruses, Internet rumors or marketing campaigns spread very rapidly. \textcolorblackFor instance, the ongoing COVID-19 pandemic [92, 96] caused an unprecedented disruption in world’s functioning. Previously, \textcolorblackan outbreak of African swine fever caused by Asfivirus in China [97] in 2018 posed the risk of spreading globally. In 2017, a worldwide cyberattack by the WannaCry ransomware cryptoworm was estimated to have affected more then 200,000 computers across 150 countries, with total damages ranging from hundreds of millions to billions of dollars [66].

\textcolor

blackVery often, spreading processes do not diffuse on their own, but instead show a complex dynamics, characterized by a collaboration or competition between them. An example for a collaborative spreading process is the co-infection of HIV and tuberculosis, the latter being \textcolorblacka major factor that influences death rates of AIDS patients [81]. The risk of developing tuberculosis is estimated to be 16-27 times higher in people living with HIV than among uninfected individuals [10]. AIDS patients are more susceptible to tuberculosis, due to their weakened immune system, and tuberculosis can also activate the replication of HIV virus. Epidemiological studies have also shown that the co-infection also exists between HIV and many other diseases including malaria parasites [39], \textcolorblackherpes [69], fungal [49] and bacteria [6], and between Zika and Dengue viruses [23]. \textcolorblack\textcolorblackAn example of asymmetric collaborative spreading is given by Hepatitis D that can only transmit to people who already infected with Hepatitis B.

An example for a competitive process can be the spreading of rumors. Recently, the “Anti-Vaxx movement” in the US, has attracted the attention of the public, for instance through Twitter messages, leading to the belief of a growing number of parents that vaccination is a violation of human-rights, that vaccines cause autism, brain damage and do not benefit the health and safety of society. The idea has been spreading rapidly through social media. As a result, the measles virus, which was declared to be eliminated in 2000, is making a comeback [37]; consequently, the World Health Organization declared Vaccine hesitancy as one of the top 10 global threats [20] and social network platforms have been requested to block the spread of related information [85]. \textcolorblackA realistic competitive scenario is given by a situation where “valid information” and “unsubstantiated rumors” spread on the network simultaneously, and where individuals exposed to one tend to believe in its content and are less susceptible to the other.

\textcolor

blackTwo important problems that recurrently arise in the analysis of spreading processes are forecasting the dynamics, and the optimal use of resources to control the dynamics, for instance \textcolorblackto maximize or minimize the spread. Forecasting is based on a probabilistic modeling and inference of the system state, such as prediction of infection probability over time for given initial conditions and interaction type. Optimization is often refereed to resource allocation tasks such as the initial choice of best spreaders or to the best vaccination strategy to contain an outbreak of a disease. In this paper, we develop novel methods that address current gaps pertaining to both inference and optimization of interaction spreading processes.

There exists a large body of work on probabilistic modeling of spreading processes. Most commonly studied models include single spreading processes that follow the SIR (Susceptible, Infected, Recovered/Removed), SI (Susceptible, Infected), and SIS (Susceptible, Infected, \textcolorblackSusceptible) dynamics [73, 45], where variables can take up a small set of \textcolorblackstatuses such as $S$, $I$, $R$ and transition from one \textcolorblackstatus to another depending on their original \textcolorblackstatus and that of their neighbors. \textcolorblackExact prediction of the spread within these models is NP-hard [82], and therefore the dynamics has been approximately analyzed using a variety of mean-field methods, see [4, 7, 78, 74, 26, 27, 45] for a review. A mean-field method of the message-passing type that is particularly suited for approximating dynamics of continuous and discrete SIR-type models on sparse network has been introduced in [41, 84, 52], which in particular gives an exact prediction of the spread on tree graphs. When averaged over an ensemble of graphs, this method is equivalent to the Edge-based Compartmental Modeling (EBCM) method [61, 64, 83] derived using the cavity method type arguments and the correct choice of dynamic variables that allows one to close the system of equations. Yet another equivalent representation is given in terms of Dynamic Belief Propagation (DBP) equations on time trajectories that \textcolorblackwas presented in [3, 54, 30]. A framework introduced in [54] showed how starting from a DBP representation allows one to systematically derive closed-form Dynamic Message-Passing (DMP) equations for any models with unidirectional dynamics \textcolorblack(so that variable statuses cannot be revisited). In particular, the method of [54] not only recovers previously known DMP equations for simple SIR-type models [41, 64, 52, 57], but also allows one to derive DMP equations for more complex models with multiple neighbor-dependent transitions, where guessing correct dynamic variables becomes incomparably harder.

The analysis of multi-agent spreading is much more involved due the interaction between processes and its impact on the spreading dynamics. Numerical studies of multi-agent processes [11] have revealed the existence of new transitions, as a result of the cross-process interaction, and an aggressive spreading mode, which points to a percolation transition. These results highlight the risk of an unpredictable and violent outbreak in cooperative spreading scenarios. The most relevant studies to the current work focus on the analysis of multi-agent spreading in a competitive scenario on a specific network, using continuous equations similar to those of dynamic message-passing [62]; and on a two-stage infection process, which is a specific case for multi-agent spreading processes [70]. \textcolorblackAlso relevant to our work are studies of complex contagions [28, 63], characterized by the requirement for multiple transmissions before a network node changes \textcolorblackstatus. While this is not exactly the scenario of interacting processes that we analyze in this paper, its dynamics also depend on the infection-history similarly to the scenarios we examine. In this case, the interplay between topology and initial conditions may give rise to hybrid phase transitions when cascades are only possible for sufficiently prevalent initial infections. We can envisage similar phenomena in the scenarios studied here for some infection probabilities but have not observed it in the experiments carried out here as this has not been the focus of our study. Interestingly, recent work [33] shows a mapping between interacting multi-agent spreading processes and social reinforcement infections through multiple transmissions.

\textcolor

blackCompetitive [80, 62, 42, 19, 86, 87, 93] and collaborative [14, 32, 29, 5, 38, 18, 58, 77, 13, 12, 50] spreading processes have been studied in different contexts and in a variety of scenarios. The foci of many of these studies have been the fixed-point properties of the system, such as a phase diagram [42, 14, 19, 42, 87], describing regimes where one spreading process dominates the infection map or where both processes co-infect the system nodes and the type of transition between phases, epidemic thresholds and the infection cluster size [50, 13, 93, 80, 62, 32, 29, 5, 38, 18]. These analyses mostly do not require a full solution of the dynamics. Other studies focus on dynamical properties of the fraction of infected network nodes by one of the processes or both [87, 80] by investigating the corresponding differential equations to identify phenomena such as hysteresis [77] particularly in the SIS model scenario, the relation between topology and dynamics [32, 18] and the emergence of infected clusters [79] linked to temporal correlations. Additionally, most studies focus on analyzing networks of different degree distributions [18, 80, 62, 42, 14, 29, 38, 50], rather than specific network instances.

\textcolor

blackA recent attempt to extend \textcolorblackmessage passing equations to the case of cooperative epidemic spreading [65], which is most relevant to the current study, only focused on a particular case where transmission is independent \textcolorblackof the \textcolorblackstatus of the target node. \textcolorblackAdditionally, the work studies different degree distributions, rather than specific instances, and concentrates on the fixed point properties such as the phase diagram and infected cluster size, falling short of a complete description of the dynamics. Hence, unlike in the case of single dynamics, exact equations for \textcolorblackdescribing the complete dynamics of general interacting processes on tree graphs remain unknown. In this work, our first major contribution consists in a derivation of DMP equations for interacting \textcolorblackunidirectional spreading processes that are exact on tree graphs, including low-complexity message-passing equations for the case of collaborative interactions. Moreover, we study approximate schemes to these equations that result in simplified expressions that can be applied on general sparse networks.

Optimal resource deployment in various spreading settings has been mostly investigated in the case of a single spreading process. One of the most commonly studied \textcolorblackproblems is identifying the most influential spreaders, on which the deployment of resource at time zero would maximize the spread at a given end time. Most of these studies rely on the network’s topological properties and selection strategies are based on high-degree nodes [75], neighbors of randomly selected vertices [17], betweenness centrality [36], random-walk [35], graph-partitioning [16], and k-shell decomposition [46]. These approaches mostly ignore important dynamical aspects that impact on performance [8, 31]. A related approach termed network dismantling [67, 68, 9] aims at identifying the nodes which, if removed, lead to the fragmentation of the giant component and prevent the global percolation. The optimal deployment of immunization has been addressed using a belief propagation algorithm [1], based on cavity method techniques developed previously for deterministic threshold models [3, 30]. Several scenarios that incorporate the dynamical properties of the spreading process, such as the optimal seeding problem, where one allocates the set of initially infected nodes that maximizes the spread asymptotically, have been studied and analyzed [21, 43, 15]. The optimal seeding problem has been analyzed for the simple diffusion models of Independent Cascade (IC) and Linear Threshold, and the optimal seeding problem has been shown to be NP-hard [44] for both, i.e., there are no deterministic algorithmic solutions that grow polynomially with the system size. A different perspective is given by the study of scenarios with a finite time horizon as studied for the IC [22] and other spreading models [72]. Most relevant to the current study is the application of a recurrent optimization framework [56] to the DMP-based probabilistic modeling of spreading processes [54]. The framework also facilitates both open-loop resource allocation (a one-off preplanned assignment) and a closed-loop (dynamical resource deployment with feedback) under a limited remedial budget. We utilize a similar framework to investigate and optimize the dynamics of \textcolorblackmultiple spreading processes.

To the best of our knowledge, no analysis and optimization algorithms have been offered to address the general case of multi-process modeling and optimization, namely incorporating both detailed topologies and dynamical properties within a fixed time window for both inference and optimization. \textcolorblackSpecial cases, such as optimal seeding, have been addressed mostly via linearized fixed-point analysis [65] and for a simple dynamic that lends itself to single time-step optimization. Moreover, most optimization algorithms for single-agent processes \textcolorblackfollow the spread on a static network topology and cannot fully capture the intricate dynamics of \textcolorblackmultiple spreading processes\textcolorblack; they are therefore less effective for \textcolorblackthe optimization tasks \textcolorblackwe aim to solve. \textcolorblackAs a second major contribution, we build an optimization framework for both competitive and collaborative scenarios based on the derived DMP equations that enable an accurate probabilistic forecasting.

We demonstrate that the inference method that we construct in this paper provides an accurate dynamical description of both competitive and collaborative scenarios on both toy and large-scale problems; it is asymptotically exact on tree-like networks and provides a good approximation on networks with loops. We develop a related optimization algorithm for maximizing the spread within a given time window against a competing spread, as well as the containment of spreads in a collaborative spreading scenario through an optimized vaccination strategy that curbs one of the spreading processes. We demonstrate the efficacy of the suggested algorithm, offering excellent results with a scalable computational complexity.

The paper is organized as follows: In Section II, we derive exact and approximate DMP equations for general models of multiple interacting spreading processes. In Section III, we validate the efficacy of the probabilistic modeling by comparing the results with Monte Carlo simulations on synthetic and real networks. The optimization algorithm is introduced in Sec. IV for both competitive and cooperative scenarios and is tested on synthetic networks in Sec. V.1. In Section V.2 we apply the optimization algorithm to real world networks for demonstrating its usefulness to more realistic scenarios. These include both competitive and collaborative cases, and the optimal deployment of vaccines to contain an epidemic. A summary and outlook are provided in Sec. VI.

II Model and Dynamic Message-Passing Equations

\textcolor

blackSpreading models studied in this work are based on the \textcolorblackdiscrete-time SI process, where a couple of spreading agents are active in parallel and interact with one another. The implication is that the \textcolorblackstatus of a network vertex determines its susceptibility to be infected by either (or both) of the spreading processes. For instance, in mutually exclusive competitive processes that describes, for instance, the battle for public opinion, once a vertex has been infected by one process it cannot be infected by any other and retains its \textcolorblackstatus, which is also termed a “cross-immunity” [70]. In a collaborative spreading scenario that describes, for example, the spread of multiple diseases, being infected by one process increases the susceptibility of being infected by another according to some predefined conditional probability. Another variant of the model considered here is that of vaccination in the presence of collaboratively spreading diseases, where vaccination against one agent affects the spread of both processes. Although the framework for the various scenarios is similar, it does include some important differences and will therefore be developed separately. \textcolorblackIt is important to point out that the introduction of cooperative/competitive spreading processes cannot simply be reduced to a stochastic process with more \textcolorblackstatuses\textcolorblack, which would include co-infection statuses; the new interactions between \textcolorblackstatuses complicate the dynamics due to the dependence of the interaction probability on the \textcolorblackstatus of neighboring variables.

\textcolor

black In both competitive and collaborative scenarios, we assume that two SI-type processes are spreading in discrete time on a graph $G\!=\!(V,E)$ comprising the set of vertices $V$ and edges $E$ such that each node $i$ can generally be found in one of four \textcolorblackstatuses at any time \textcolorblackstep $t$: susceptible ($\sigma_{i}(t)\!=\!S$), infected by disease $A$, ($\sigma_{i}(t)\!=\!A$), infected by disease $B$, ($\sigma_{i}(t)\!=\!B$), or activated by both processes $A$ and $B$, ($\sigma_{i}(t)\!=\!AB$). In what follows, we define the spreading model in two scenarios, and derive the corresponding DMP equations.

II.1 From Dynamic Belief Propagation to Dynamic Message Passing

\textcolor

black DMP equations can be hard to guess beyond simplest models, but can be systematically derived starting from the general Dynamic Belief Propagation algorithm that approximates probabilities of time trajectories of individual nodes [54]. This is the approach that we will be undertaking here to obtain exact equations on tree graphs. For the two-processes dynamics considered in this work, the dynamics of a single node $i$ is fully described by a pair of activation times, $(\tau_{i}^{A},\tau_{i}^{B})$, where $\tau_{i}^{A}$ denotes the first time when node $i$ is found in the \textcolorblackstatus $A$, and similarly for $\tau_{i}^{B}$. For instance, $\tau_{i}^{A}=0$ means that node $i$ was initially in the active \textcolorblackstatus $A$, and we will denote by $\tau_{i}^{A}=*$ the situation where node $i$ did not get $A$-activated before some final observation time, in other words $*$ absorbs all the history that happens after the end of the observation window. For the convenience of presentation, in what follows we will consider two separate “observation windows”, for $A$ and $B$ processes.

\textcolor

black The starting point for our derivations are the general DBP equations [40, 2, 54] on the interaction graph, where the goal is to approximate the probability $m^{i}_{T_{A}+1,T_{B}+1}(\tau_{i}^{A},\tau_{i}^{B})$ that node $i$ exhibits a trajectory $(\tau_{i}^{A},\tau_{i}^{B})$ during the observation time window of length $T_{A}$ for process $A$ and $T_{B}$ for process $B$ \textcolorblack(we keep the flexibility of having two separate time windows although in most situations we use $T_{A}=T_{B}=t$). Exact equations that compute \textcolorblackthe probability $m^{i}_{T_{A}+1,T_{B}+1}(\tau_{i}^{A},\tau_{i}^{B})$ of a given time trajectory $(\tau_{i}^{A},\tau_{i}^{B})$ of node $i$ are explained in Appendix A. Due to the properties of Belief Propagation algorithm [60, 88], the fixed point solution of the DBP equations is guaranteed to be exact on tree graphs, and provides good estimates of marginal probabilities on loopy but sparse graphs [54].

\textcolor

black Given the computed value of the marginals $m_{T_{A},T_{B}}^{i}(\tau_{i}^{A},\tau_{i}^{B})$, one can straightforwardly define quantities of interest, such as probabilities for a given node $i$ to be found in a given \textcolorblackstatus:

[TABLE]

\textcolor

blackThe definition of $P_{A}^{i}(t)$ ($P_{B}^{i}(t)$) indicates that at time $t$ the vertex $i$ is $A$-activated ($B$-activated) irrespective of the other process. Therefore $P_{A}^{i}(t)$ ($P_{B}^{i}(t)$) represents the probability that at time $t$ vertex $i$ is in status $A$ ($B$) alone, or in status $AB$. Initial condition probabilities where the node is exclusively found is status $A$ ($B$) will be denoted as $P_{A^{*}}^{i}(0)$ ($P_{B^{*}}^{i}(0)$). In principle, one can solve the DBP equations to obtain $m^{i}_{T_{A}+1,T_{B}+1}(\tau_{i}^{A},\tau_{i}^{B})$, and then use the expressions (36)-(39) to obtain final aggregated expressions for the dynamic messages $P^{i}_{S}(t)$, $P^{i}_{A}(t)$, $P^{i}_{B}(t)$, and $P^{i}_{AB}(t)$. However, due to the generality of the DBP equations that are valid for any dynamics, this approach may not be the most efficient one: computing a single marginal $m^{i}_{t+1,t+1}(\tau_{i}^{A},\tau_{i}^{B})$ may require as many as $O(t^{2d})$ operations for a single marginal, where $d$ is the degree of the node $i$. On the other hand, for the concrete dynamics such as the processes considered here that often have a special structure, it is beneficial from a computational point of view to explore this structure in order to drastically reduce the complexity of computing the dynamic marginals $P^{i}_{S}(t)$, $P^{i}_{A}(t)$, $P^{i}_{B}(t)$, and $P^{i}_{AB}(t)$. Where this special structure is exploited to produce closed-form algebraic equations for computing $P^{i}_{S}(t)$, $P^{i}_{A}(t)$, $P^{i}_{B}(t)$, and $P^{i}_{\textcolor{black}{AB}}(t)$ with low algorithmic complexity, the resulting computational scheme will be referred to as Dynamic Message-Passing equations for the given process.

\textcolor

black Therefore, the procedure we adopt below for deriving DMP equations will be as follows: (i) specify DBP equations based off a given two-processes dynamics; (ii) where possible, exploit the structure in the dynamics to derive low-complexity closed-form DMP equations that iteratively compute the quantities of interest $P^{i}_{S}(t)$, $P^{i}_{A}(t)$, $P^{i}_{B}(t)$, and $P^{i}_{AB}(t)$ starting with the algebraic definitions (36)-(39), which will inherit exactness of prediction of tree graphs; and (iii) to further reduce the computational complexity or the algebraic form of the exact equations, derive approximate DMP equations that could be used as an algorithm for inference or optimization problems on general graphs.

II.2 Mutually Exclusive Competitive Processes

\textcolor

blackThe dynamics in mutually exclusive competitive processes can be made explicit by listing the allowed transitions and their respective probabilities at every discrete time step:

[TABLE]

In other words, the two infection processes $A$ and $B$ are mutually exclusive: any node can be infected by a neighboring nodes, assuming one of the two \textcolorblackstatuses, but once infected by one of the two processes, it cannot change its \textcolorblackstatus. Since the infection is \textcolorblackbased on a two-vertex interaction through an edge, the processes in (5) seem deceptively as two completely independent parallel processes; however, they clearly interact through the graph topology and the exclusivity of the adopted \textcolorblackstatuses. Indeed, we assume that at any given time step, the infection probabilities of process $A$ (denoted $\alpha^{A}_{ji}$) or process $B$ (denoted $\alpha^{B}_{ji}$) are \textcolorblacktreated as independent, but the probability of being infected by both $A$ and $B$ simultaneously is forced to be 0 \textcolorblack(thus creating probabilistic dependence between the two processes). For closing the \textcolorblackupdate rule we need to define what happens in the case where both processes jointly infect the vertex \textcolorblackin the same time step, resulting in an invalid \textcolorblackstatus. There are \textcolorblackmany possible ways to deal with this case that could be accommodated in both analysis and simulations, depending on the needs of a particular application. For the sake of simplicity, we consider the rule where the probabilities of transitioning to either of the \textcolorblackstatuses $A$ or $B$, or staying in the \textcolorblackstatus $S$ is proportionally renormalized in such a way that they sum to one. Alternatively, in simulation one could think of this procedure as resampling in the case where the joint infection occurs: If a joint infection \textcolorblackstatus by both processes $A$ and $B$ is sampled, it is rejected and resampling is carried out. \textcolorblackThis is done repeatedly until a valid \textcolorblackstatus without progressing the dynamics, such that no spurious probabilistic dependencies emerge. According to the vanilla dynamic rules defined above, the probability of transition to the \textcolorblackstatus $A$ from \textcolorblackstatus $\textcolor{black}{S}$ for a node $i$ is given by

[TABLE]

where $\partial i$ denotes the set of neighbors of node $i$, and $\mathbbm{1}$ is an indicator function. Similarly, define

[TABLE]

Then under the resampling procedure explained above, the final renormalized transition probabilities at time step $t$ read:

[TABLE]

where notation \textcolorblack$\sigma_{i}(t)=A/B/S$ refers to a node \textcolorblack$i$ at time $t$, being in one of the \textcolorblackstatuses $A$, $B$ or $S$.

\textcolor

black Given expressions (41)-(43), we can straightforwardly specify the respective DBP equations for the mutually exclusive competing dynamics of two processes that are given in Appendix B. In Table 1, we numerically verify that the DBP equations are exact on tree graphs. With a naive implementation, DBP marginals can be computed explicitly in time $\max(O(|E|T^{2(c-1)}),O(NT^{2c}))$, where $c$ is the maximum degree of the graph, and $T$ is the final observation time. It is important to note that transition probabilities defined as in (41)-(43) reflect the complexity of interaction between processes that results from the renormalization of probabilities. Indeed, these expressions depend on the particular realization of \textcolorblackstatuses for all neighbors, and hence lack any iterative structure at each time step. Due to this lack of structure, exact low-complexity DMP equations can not be straightforwardly derived for the chosen dynamics.

\textcolor

black Although DBP equations for the mutually exclusive competing processes are exact on trees (see Table 1), their polynomial but potentially high computational complexity makes them a less attractive choice in applications. To address this issue, we notice that real application problems are typically defined on sparse, but non-tree graphs\textcolorblack. Even if low-complexity DMP equations were available, they would still yield only approximate solutions on general loopy network instances. This observation motivates us to search for a tractable approximation of the message-passing equations on tree graphs: if the approximation is good, the resulting error on tree-like loopy networks may be similar compared to the application of exact DBP or DMP equations. In Appendix C, we implement this strategy and derive approximate DMP equations for the mutually exclusive competing scenario with a computational complexity $O(|E|T)$, where $|E|$ is the number of edges in the graph and $T$ is the observation window, which makes them scalable for very large sparse networks with millions of nodes. The approximation that we use is inspired by the fact that in the absence of renormalization the dynamics of each processes follows the dynamics of the usual SI-type process, and hence it is natural to try to perform the renormalization procedure at the level of dynamic marginals. In what follows, we present numerical tests that illustrate the accuracy of the employed approximation.

II.3 Collaborative Process

\textcolor

blackIn the two-processes collaborative scenario considered here, when a node is infected by one process, its susceptibility to the activation by another process increases (or decreases). Unlike in the mutually exclusive case, a node can be infected by both processes and the influence of the different processes is not necessarily symmetric. The dynamics in collaborative processes can be made explicit by listing the possible transitions and their respective probabilities at every discrete time step:

[TABLE]

\textcolor

black In this scenario, \textcolorblackstatus $AB$ can be regarded as a combination of \textcolorblackstatus $A$ and \textcolorblackstatus $B$. The non-triviality of the interaction between both processes comes from the fact that $\alpha^{AB}_{ij}$ and $\alpha^{BA}_{ji}$ are different from $\alpha^{A}_{ji}$ and $\alpha^{B}_{ji}$, respectively: when a node is infected by one process and becomes more (or less) vulnerable to another, and vice versa. Notice that unlike the mutually exclusive scenario where the \textcolorblackstatus $AB$ is forbidden, under the general collaborative scenario the process $S(i)\longrightarrow AB(i)$ is allowed, and in discrete time the rule

[TABLE]

follows from the transition rules above, simply as co-activation that happens at the same time.

\textcolor

black Following the scheme outlined above, we can start by forming a dynamic transition kernel that encapsulates the various transition rules and allows us to write the DBP equations for this spreading model. The resulting DBP equations are given in Appendix D. The special structure of the dynamic kernel is written as a sum of possible transition sequences factorized over the neighbors of a given node, making it possible to derive low-complexity and exact DMP equations for the collaborative model. The algebraic form of equations for the dynamic marginals is given a simple and intuitive meaning. The probability of finding node $i$ in the \textcolorblackstatus $S$ can be written as

[TABLE]

where $\theta_{A,B}^{k\rightarrow i}(t,t)$ is an aggregated dynamic message defined through the fundamental messages on time trajectories:

[TABLE]

From the definition of $\theta_{A,B}^{k\rightarrow i}(t,t)$, it is easy to “read off” its physical meaning: it corresponds to the probability that node $k$ did not send activation signals \textcolorblackeither $A$ \textcolorblackor $B$ before time $t$, while $i$ follows a fixed $(\tau_{i}^{A},\tau_{i}^{B})=(*,*)$ dynamics, i.e.\textcolorblack, it does not activate until time $t$ (\textcolorblacksee Appendices A and D for more details). This conveys the following meaning to the expression (17): $P_{S}^{i}(t)$ is given by the probability $i$ is in the \textcolorblackstatus $S$ at initial time, times the probability that \textcolorblacknone of its neighbors has activated it with any of the processes until time $t$ (which exactly factorizes over neighbors on a tree graph).

\textcolor

black In a similar way, the marginals corresponding to \textcolorblackstatuses $A$ and $B$ can be expressed as activations by the time $t$:

[TABLE]

where reduced marginals $\mu_{A}^{i}(t^{\prime})$ and $\mu_{B}^{i}(t^{\prime})$ are defined as follows:

[TABLE]

Finally, using the normalization of probabilities \textcolorblack(notice that by definition $P_{AB}^{i}(t)$ is contained in both $P_{A}^{i}(t)$ and $P_{B}^{i}(t)$), we finally get

[TABLE]

\textcolor

blackThe exact form of the DMP equations for collaborative processes along with detailed derivations are provided in the Appendix D. Due to exploitation of the structure properties, the DMP equations have a much lower computational complexity $O(|E|T^{2})$ for a final observation time $T$ compared to the DBP equations, while still providing exact predictions on tree graphs. We numerically verify this fact on a number of problem instances as shown in Table 2.

\textcolor

black In the search for simplified equations that are easier to use for inference and optimization purposes, as well as reduced computational complexity in time, we implement a similar strategy to the one used for the mutually exclusive scenario and derive approximate equations that could be used in lieu of exact DMP equations on general graph. To do so, we notice that the special case of transmission probabilities $\alpha_{ij}^{AB}=\alpha_{ij}^{A}$ and $\alpha_{ij}^{BA}=\alpha_{ij}^{B}$ for collaborative dynamics (12)-(15) corresponds to non-interacting spreading processes: activation by one process does not change the activation dynamics for the other. In this case, the DMP equations should simplify into the product of two independent Susceptible-Infected (SI) like processes:

[TABLE]

We use this observation to produce a simplified version of DMP equations, expanding the exact equations for interacting spreading processes around the non-interacting point. We keep certain first-order corrections in the update equations only, so that the resulting equations are similar to the SI-type equations. This approximate version is expected to be good as long as $\alpha_{ij}^{AB}$ and $\alpha_{ij}^{A}$, $\alpha_{ij}^{BA}$ and $\alpha_{ij}^{B}$ are similar. \textcolorblackIt would clearly break down if $\alpha_{ij}^{A}=0$ or $\alpha_{ij}^{B}=0$ while the other infection probabilities remain finite. Full derivation and expressions for the approximate DMP equations are provided in Appendix E.

III Inference using Approximate Dynamic Message-Passing Equations

\textcolor

black In Section II, we considered DBP and DMP equations that are exact on tree graphs, which follows from their derivation and supporting numerical checks. In this Section, our goal is to numerically establish the validity of approximate DMP equations introduced in the previous Section. These equations enjoy an improved computational complexity and a simpler algebraic form, and are expected to provide a good approximation in the regimes discussed above.

To validate the \textcolorblackapproximate DMP equations obtained for competitive/collaborative processes we test the accuracy of the inferred node values against numerical results obtained via Monte Carlo simulations. Testing is carried out on both synthetically generated networks and real instances.

III.1 Inference in Competitive Spreading

Validation is carried out on two synthetic networks generated by the package NetworkX, a tree network and a network with loops; and on a real world undirected benchmark Polbooks [48] network. The latter is provided as an example of a sparse network. With the given initial condition, we apply both Monte Carlo simulation and the DMP method to all models.

Exhaustive numerical experiments reveal that DMP-based inference provides an accurate marginal posterior probabilities for the variable \textcolorblackstatuses, with the exception of very high infection parameters values (very close to 1). Therefore, we do not examine cases with extreme infection parameter values in the examples provided. The location of seeds initializing the processes and target nodes to be observed are chosen randomly and have no significance.

In Fig. 1 we show two synthetically generated networks (a) a toy tree network of 10 nodes, (b) a network of 10 nodes with loops and (c) the Polbooks [48] network with 105 nodes. The choice of network has no significance, it is a standard benchmark network used in the literature, of books about US politics sold by Amazon [48]. The network comprises two sparsely connected hubs, where edges represent books (vertices) bought jointly by the same individuals. To validate the efficacy of DMP in modeling competitive scenarios we compare results obtained from running equations (A4)-(A9) against results obtained using Monte Carlo simulations. Simulations are carried out 10 times for gathering statistics, each round includes $10^{3}$ samples per node (about $10^{5}$ samplings in total, depending on the network size). The parameters used in the toy model tests are $\alpha_{A}\!=\!0.3$ and $\alpha_{B}\!=\!0.7$, and we observed the marginal posterior probabilities $P^{i=3}_{A}(t)$ in both Figs. 2(a) and (c), and $P^{i=3}_{S}(t)$ in Figs. 2(b) and (d) for the tree and loopy graphs, respectively. The seeds initializing the processes are placed on nodes 7 for process $A$ and on nodes 2 for process $B$. The choice of these particular nodes is arbitrary and has no significance. The results obtained show excellent agreement between theory and simulations. We also tested the accuracy of the method on the benchmark Polbooks network as shown in Fig. 2(e) and (f) for the parameters $\alpha_{A}\!=\!0.2$, $\alpha_{B}\!=\!0.2$. In this case, process $A$ starts from nodes 1 and 2, process $B$ from nodes 4 and 37 (again, both are arbitrary choices). The observed probabilities $P^{i=0}_{A}(t)$ in Fig. 2(e) and $P^{i=0}_{S}(t)$ in Fig. 2(e) show good agreement between theory and simulations.

III.2 Inference in Collaborative Spreading

Similar experiments were run for a collaborative process on the toy tree network, graphs with loops and the benchmark Football network. The latter is an undirected network of American football games between colleges during season Fall 2000 [25] which is quite uniformly connected.

The results are shown in Fig. 3 for the various cases.

The experimental results indicate that modeling based on \textcolorblackapproximate DMP equations is very accurate on tree-like networks, as expected for message-passing algorithms. It is less accurate on small loopy graphs at longer times, as expected, due to the small loops that violate the cavity method’s assumption (the specific 10-nodes network used includes 2 small loops). This effect is suppressed to some extent on larger networks where loops are typically longer as demonstrated on the real benchmark network examples.

\textcolor

blackAs discussed in Appendix E, the approximation quality is expected to degrade when $\alpha_{A}$ and $\alpha_{B}$ are very different from $\alpha_{AB}$ and $\alpha_{BA}$, respectively. Indeed, in this case the prediction of the dynamics by the approximate DMP equations can become inaccurate. To illustrate this point, let us consider an extreme example on a chain of three nodes and two links, with node 2 connected to nodes 1 and 3. Let us assume that $\alpha_{A}=1$, $\alpha_{B}=0$, $\alpha_{AB}=\alpha_{BA}=1$, i.e.\textcolorblack, that the infection $B$ can be transmitted only to nodes that are already infected with $A$. Interestingly, this scenario is relevant for several disease pairs such as Hepatitis D that can only be transmitted to individuals already infected with Hepatitis B. Given the initial condition $P_{A}^{1}(0)=1$, $P_{S}^{2}(0)=1$ and $P_{B}^{3}(0)=1$, node $1$ will become infected with disease B at time $t=3$: it first infects node $2$ with process $A$ at time $t=1$; this leads to infection of node $2$ by infection $B$ coming from node $3$ at time $t=2$; and finally node $2$ transmits infection $B$ to node $1$ at time $t=3$. This example represents an extreme case where the approximate DMP equations are not valid, and indeed preclude node $1$ from being infected by process $B$, which illustrates that they may not be exact even on tree graphs when the approximation criterion is not satisfied. However, it is easy to check that exact DMP equations provide an accurate answer in this case as well, as it should be. This counter example re-iterates the trade-off between the exactness and the computational complexity between exact and approximate DMP equations discussed in the previous section. However, as illustrated on numerical examples, in applications with \textcolorblackplausible parameters one could expect that the approximate DMP method provides a good description for cooperative spreading processes on sparse networks.

IV DMP-based Optimization Method

Competition and collaboration of spreading processes on graphs can be optimized through a judicious use of resource. We will demonstrate how managing a spreading process against an adversarial competing agent can be optimized within a given time frame; and how joint collaborative processes can be affected through best use of resource. The latter can take the form of spreading maximization while making use of the process interdependencies, or of containment through vaccination to impact on the spread of both processes.

We outline a general procedure for optimization in this section. Details for specific optimization problems we address here are given in \textcolorblackAppendix F for competitive processes and in \textcolorblackAppendix H for collaborative processes.

The core approach for optimization is based on a discretized variational method, whereby a functional over a time window (Lagrangian) is optimized through changes in control parameters throughout the time interval. The dynamics, resource constraints, initial conditions and other restrictions on the parameters used are enforced through the use of Lagrange multipliers. A similar method is used in optimal control.

We denote the components of the Lagrangian function used in a similar way to [56]:

[TABLE]

where in the ${\cal L}$ is the Lagrangian function, ${\cal O}$ is the objective to be optimized, ${\cal B}$ is the budget or constraints on the resource used, ${\cal I}$ represents the component that forces initial conditions, ${\cal P}$ are restrictions on the probabilities used and ${\cal D}$ represent the \textcolorblackdynamical constraints that here take form of approximate DMP equations. All of the terms ${\cal B}$, ${\cal I}$,${\cal P}$ and ${\cal D}$ are forced through the use of Lagrange multipliers.

For different problems, the constraints and objectives vary. We take the competitive process as an example. In this problem, we would like to find an optimal allocation of a limited number of spreaders (representing a budget, potentially time dependent) for process $A$, which minimizes the spreading of process $B$. These could represent a competition in a political or commercial setting. The objective function in this case is

[TABLE]

where $T$ is the end of the set time window. Our goal is to maximize this objective function, thus minimizing the spread of process $B$. The resources at our disposal for seeding nodes with process $A$ are represented by the budget constraint at time zero (although more elaborate budget constraint could be accommodated as in [56] and in some of the examples that follow)

[TABLE]

where $\nu^{i}(0)$ is the deployment of a fraction of the budget for process $A$ on node $i$. The budget constraint is forced through the Lagrange multiplier $\lambda^{Bu}$, such that

[TABLE]

The fraction/probability $\nu^{i}$ is kept within a certain range, determined by the upper and lower bounds $\overline{\nu}$ and $\underline{\nu}$, respectively. Also this is enforced using a Lagrange multiplier

[TABLE]

The term ${\mathcal{D}}$ is given by enforcing the approximate DMP equations for the competitive case using a set of Lagrange multipliers as detailed in \textcolorblackAppendix H for the competitive case. The remaining term ${\mathcal{I}}$ forces the initial conditions for the dynamics.

The extremization of the Lagrangian (27) is done as follows. Variation of $\mathcal{L}$ with respect to the dual variables (Lagrange multipliers) results in the DMP equations starting from the given initial conditions, while derivation with respect to the primal variables (control and dynamic parameters) results in a second set of equations, coupling the Lagrange multipliers and the primal variable values at different times. End conditions for the forward dynamics provide the initial conditions for the backward dynamics.

We solve the coupled systems of equations by forward-backward propagation, a widely used method in control and detailed in [56]. This method has a number of advantages compared to other localized optimization procedures. It is simple to implement, of modest computational complexity \textcolorblack$O(ET)$, where $E$ the number of edges in the graph and $T$ the time window) and does not require any adjustable parameters. The forward-backward optimization provides resource (budget) values to be placed at time zero (or of at any time within the time window if we so wish) in order to optimize the objective function, e.g., that of Eq. (28). \textcolorblackOne potential drawback of the method is the possible non-convergence of the dynamics to an optimal solution. This can be mitigated to some extent by solving the equations for the backwards dynamics using other available solvers and by storing the best solutions found over time, \textcolorblackor approaching a fixed point via gradient descent. In general, as the functions used become more nonlinear it will become more difficult to obtain optimal solutions, although we have not experienced significant problems in the cases studied here.

V Numerical Study of the Optimization Algorithm

To validate and demonstrate the efficacy of the optimization method we carry out experiments on both synthetic and realistic networks. Before embarking on a large-scale application we study the performance of the derived method on a tree-like network of 30 nodes.

V.1 Validation of the Optimization Algorithm

To validate the DMP optimization algorithm on a problem that could be exhaustively studied and intuitively presented, we restrict the study to a small exemplar tree-like synthetic model. Moreover, we select a small number of nodes (3) on which resource could be deployed. The objective is to maximize the spread of both agents or to minimize the spreading of one of them (the disease control scenario) in both competitive and collaborative processes.

The 30-node network used for carrying out these experiments is presented in Fig. 4. Comparison between the DMP-based optimization algorithm and the exhaustive search is implemented in the following way: We consider a scenario where the entire budget is available at time $t\!=\!0$. The optimization problem minimizes the spreading of process $B$ in a competitive processes through judicious budget allocation of the seeds for process $A$, i.e.\textcolorblack, we aim at minimizing $\sum_{i}P_{i}^{B}(T)$, where $T$ is the end time of the process.

In the experiments, resources for process $A$ were deployed on three nodes: [math], $19$ and $28$ (determined by the choices for nodes [math] and $19$ due to the total budget constraint). The fixed seed for process $B$ is node $3$; these choices are arbitrary and insignificant. The objective function landscape has been explored by sampling for different parameter values as denoted by the green points in Fig. 5(a). The collaborative scenario with specific infection parameters is plotted in Fig. 5(b). A different set of controllable nodes is presented in Fig. 5(c) and Fig. 5(d) for the competitive and collaborative cases, respectively. The maximum value obtained for the objective function is contrasted with the results obtained using the DMP-based optimization procedure, marked by the red dot, in all cases. It is clear that there is good agreement between the DMP-optimal values obtained and the optima discovered through sampling\textcolorblack, although in some of the cases they are not identical. The corresponding objective function values for Fig. 5 are denoted on the various subfigures.

V.2 Results on Real Networks

V.2.1 Competitive processes

We study the performance of our optimization algorithm on networks assuming that all resources are available at time zero. Other optimization strategies where resources’ availability is time-dependent could also be considered [56] as shown later on. One of the problems in comparing the performance of our method to other approaches is the limited number of dedicated competing optimization methods for the multiple agents scenario, especially in cases where resource availability and its deployment are spread over the whole time window. We therefore compare the DMP-optimized spreading process with known heuristics for the single agent scenario such as a uniform allocation of resources over all nodes (except seed nodes induced with process $B$), the high degree deployment strategy [76] HDA and K-shell decomposition [46]. Free spreading refers to an uncontrolled spread of process $B$. Blocking refers to the allocation of a competing agent that is non-infectious ($p\!=\!0$). In all cases, resources are used to contain the spreading of process $B$.

Multi-agent spreading processes exist in different settings, from social network, energy grid and road network to the graph of random interaction between individuals. There is no specific network type that is more relevant to the scenarios we examine. We therefore chose a set of sparsely connected benchmark networks of different characteristics to test the efficacy of our algorithm compared to competing approaches. The description of the different networks and their specific properties appear in the corresponding references. We aim to show that the suggested methods work well on all of the sparse networks examined and we therefore expect it to work well on most other sparse networks.

In these experiments we allocated seeds of process $B$ on $0.05N$ nodes at time zero, where $N$ is the total number of nodes in the network. The same amount of resource was allocated to the controlled competing process $A$. Infection parameters are $\alpha_{A}\!=\!0.2$ and $\alpha_{B}\!=\!0.3$. This choice of parameters is arbitrary and the performance for other parameter choices provides qualitatively similar results. The free parameter that forces the upper/lower limits of the resource variables (213) is set to $\epsilon\!=\!0.1$ initially and decays exponentially with iteration steps. Experimental results show that decaying $\epsilon$ in this manner leads to a improved performance, arguably since it allows one to obtain solutions which are closer to the limit values. The optimization procedure is iterated 10 times and the best result is selected. The objective is to minimize the spreading of process $B$ and the normalized total spreading ${\sum_{i}P_{B}^{i}(T)}/{N}$ at time $T\!=\!3$ is shown in Table 3. The short time window used is due to the small diameter of the networks.

In competitive processes we observe that the containment of process $B$ is carried out effectively by optimal deployment of spreading agents of process $A$, compared to a static blocking, HDA deployment, K-shell and uniform seeding. This scenario corresponds to marketing, rumor spreading/fake news and opinion setting. We also evaluate the improvement obtained for a given budget in blocking the spread, which allows one to allocate the appropriate budget for containment (e.g., addressing the spread of fake news or antivaxxing rumors by releasing verifiable information). Clearly, key factors in determining the spread are the actual infection parameters associated with the various processes, which can be obtained through data analysis. The Football network is highly connected and we speculate that this is the reason for the success of uniformly allocating the spreading agents.

When the budget allocation is done dynamically at different times, several optimization procedures can be used. \textcolorblackConventional stochastic optimal control [24] is based on \textcolorblackplanning ahead for the entire time window, taking into account future uncertainties. An optimal solution in discrete time can be stated through a solution of the Bellman equation that in our setting would result in an algorithm with a high computational complexity due to the curse of dimensionality [59]. To mitigate this issue, we adopt a procedure that reminds updates in the closed loop control, which updates the resource allocation at every step based on the information on the realization of the stochastic dynamics. This approach does not strictly guarantee the solution optimality, but quantitatively takes into account realization of uncertainties while keeping computational complexity under control. We carry out the optimization for the end time $T$ at each step, on a shrinking time window, incorporating the newly available information. Considering a case where one process spreads randomly while for the other a unit budget is available per time step to be deployed optimally. Given the information at hand for each time step, one may then consider the following dynamic allocation strategies: (a) DMP greedy deploys the unit to optimize the objective function for the next time step; while (b) DMP optimal optimizes the objective function for the end-time $T$, while incorporating information available at each time step. Note that the latter represents a \textcolorblackclosed-loop like optimization in the sense that updated information on realization of dynamics at each time step is incorporated to produce a better resource allocation plan.

To demonstrate the efficacy of our method and the differences between the DMP-greedy and DMP-optimal resource deployment, we use the Football network [25] and infect node 1 at time 0 (blue point left of the center, total budget for $B$ is 1); we then deploy optimally a budget of 1 for process $A$ at each time step $t\!=\!1,2,\ldots,5$. The infection parameters used are $\alpha_{A}\!=\!\alpha_{B}\!=\!0.7$. The results are shown in Fig. 6 where the heat bar represents the dominating process per node through the value $P_{A}^{i}(t)\!-\!P_{B}^{i}(t)$. Red and blue represent dominating processes $A$ and $B$, respectively. It is clear that DMP-optimal (b) is much more effective than DMP-greedy (a) in restricting the spread of process $B$ by maximizing the spread of process $A$. Numerical comparison between the two methods for the same network and conditions are presented in Table 4. It is clear that while DMP-greedy is successful at earlier time steps, DMP-optimal minimizes the spread of process $B$ at $T\!=\!5$. A second example on a more densely connected network is given in \textcolorblackAppendix G.

To demonstrate the improvement in reducing the spread of agent $B$ given a budget $b_{B}$ against an optimally deployed budget $b_{A}$ of process $A$, we plotted the ratio of infected $A$ and $B$ \textcolorblackstatuses at time $T\!=\!5$ against the budget ratio $b_{A}/b_{B}$. This has been done for the Lesmis network [47] (77 nodes, representing the co-appearance of characters in the novel Les Miserables) with infection probabilities $\alpha_{A}\!=\!\alpha_{B}\!=\!0.5$ and time window $T\!=\!5$. The results shown in Fig. 7(a) have been averaged over 5 instances for both uniform and optimal DMP-based deployment. It is clear from the figure that the optimal DMP-based deployment provides much better results than uniform deployment, which exhibits a linear increase. The saturation for high ratios is limited by the size of the graph.

We observed several instances in which the ratio of infected probabilities exhibits a fast transition at specific points. Figure 7(b) shows a fast transition towards a process-$A$ dominated network as the ratio between the budgets allocated exceeds a certain value $b_{A}/b_{B}\!\approx\!0.28$. The $y$-axis represents the ratio of process probabilities $\frac{\sum_{i\!=\!1}^{N}P_{A}^{i}(T)}{\sum_{k\!=\!1}^{N}P_{B}^{k}(T)}$. These results clearly depend on the topology, budgets and infection rates used. The example given here is based on the Lesmis network with infection probabilities $\alpha_{A}\!=\!0.5$ and $\alpha_{B}\!=\!0.7$, fixed budget $b_{B}\!=\!1$, time window $T\!=\!5$ and optimal DMP-based deployment.

To evaluate the interplay between infection parameter values and budget allocated to each of the processes, we studied a competitive case on the Lesmis network, with a varying ratio between budgets ($b_{B}\!=\!1$, $b_{A}\!=\!0.5,\ldots,4.5$) and infection probabilities ($\alpha_{B}\!=\!0.7$, $\alpha_{A}\!=\!0.1,\ldots,0.7$). The points where the processes end up with equal probabilities at $T\!=\!5$ are plotted in Fig. 7(c) for the DMP-optimized deployment of resource $A$ (green line) and for uniform deployment (red line). The results are averaged over randomly chosen 5 initial position for the seed of $B$. From the figure it is clear that in this case DMP-optimized deployment can effectively mitigate significantly inferior infection rates or budget ratios (area above the green curve); only when both ratios are very low the $B$ process dominates the network after $T$ steps (yellow area). Uniform deployment results in much inferior performance (area above the red line).

V.2.2 Collaborative processes

Results obtained for collaborative scenarios, shown in Table 5, exhibit a similar behavior, showing that collaborative processes can be optimized to spread quickly due to the mutually-supportive role played by the two processes. In this case, some nodes ($0.05N$) are infected by process $B$ and we allocate a given budget of process $A$ such that the joint spreading will be maximized and the number of non-infected nodes (in \textcolorblackstatus $S$) minimized. This could represent, for instance, the spread of opinions on the basis of political affiliation. Also in this case the DMP-based optimization algorithm works well, with the exception of the Football network where uniform spreading seems to be successful, arguably due to the same reasons as in the competitive case.

The final set of experiments we carried out relates to the vaccination policy in collaborative spreading processes. In this case, two collaborative processes $A$ and $B$ spread throughout the system and the task is to minimize the spread through a vaccination process to one of the two, say $B$. The role of vaccination in this case is in reducing the infection probability by the level of vaccination deployed. For instance, we adopt a simplistic model whereby if the vaccine deployed at a given node is $b$, its probability of being infected will reduce to $\alpha_{B}\!-\!b$. Similarly, the infection parameter $\alpha_{BA}$ will also decrease by the same amount to $\alpha_{BA}\!-\!b$. Other infection models can be easily accommodated. The optimization task is in the deployment of a given vaccination budget such that the fraction of non-infected sites ${\sum_{i}P_{S}^{i}(T)}/{N}$ at time $T$, is maximized.

The experimental results demonstrate that the DMP-based optimization shows excellent performance.

Generally, for the different network topologies and size, DMP-based optimization appears to be more effective than other methods. The main cases where it does not offer the best result are in complex networks with bounded connectivity fluctuation, where it has been shown that uniformly applied immunization strategies are highly effective [76].

To demonstrate the efficacy of DMP-optimal resource deployment for vaccination in the case of a collaborative spreading process we use the main road network of England [34]. \textcolorblackLivestock epidemics often spreads through the transport of infected livestock on the road network (as was the case in the 2001 Foot and Mouth epidemic in the UK). We consider the spread of highly infectious coupled processes through the road network starting from the areas of Greater London (node 10, process $A$) and Leeds (node 11, process $B$); the budget for both $A$ and $B$ is one unit and infection parameters $\alpha_{A}\!=\!\alpha_{B}\!=\!0.2$, $\alpha_{BA}\!=\!\alpha_{AB}\!=\!0.99$. The process is runs for 20 steps. The vaccination budget is one unit per time step and is effective against process $B$ only, while being ineffective for process $A$. The results shown in Fig. 8 demonstrate the efficacy of the DMP-optimal vaccination strategy aimed at minimizing $\sum_{i}{(1\!-\!P_{S}^{i}(t))}$ (b) in contrast to the free spreading of both infections (a). Blue/red represent uninfected/infected \textcolorblackstatuses, respectively; more specifically, the node color represents $1\!-\!P_{S}^{i}(t)$, where red and blue correspond to 1 and 0, respectively. As we can see in Fig. 8, the infection spreading around London remains the same, with or without the deployment of vaccine as it is mainly infected by process $A$, and hence the vaccination has no effect, while the spread emanating from the Leeds area is effectively blocked.

VI Summary and future work

Competition and collaboration between spreading processes are prevalent on social and information networks, and on interaction-networks between humans or livestock, to name but a few. To better understand the expected spread and dynamics of diseases, marketing material, opinions and information it is essential to infer and forecast the spreading dynamics. Moreover, the judicious use of limited resource will help contain the spread of epidemics, win political and marketing campaigns, and better inform the public in the battle against unsubstantiated or misleading information.

Most current studies of multi-spreading processes such as [11] heavily rely on computer simulations and heuristics with varying results from one instance to another. Obtaining reliable estimates therefore comes at a high computational cost in order to obtain meaningful statistics. We advocate the use of a principled probabilistic technique for the analysis and employ the same framework to develop and offer an optimization algorithm for best use of the resource available.

Both competitive and collaborative spreading processes are investigated here using the probabilistic DMP framework and provide an accurate description of the spreading dynamics, as validated on both synthetic and real networks. \textcolorblackFor the first time, we derive exact message-passing equations for two interacting \textcolorblackunidirectional dynamic processes, and study approximate DMP equations that benefit from a simpler form and a lower algorithmic complexity; this represents the first major novel contribution of this work. This probabilistic description facilitates the study of multi-agent spreading processes on general sparse networks at a modest computational cost, opening up new possibilities for obtaining insights into the characteristics of such scenarios.

\textcolor

blackAs a second major contribution in this work, we employ a scalable optimization framework that incorporates the DMP dynamics to deploy limited resource for an objective function to be optimized at the end of the time window. This is a difficult hard-computational problem since the resource is deployed at earlier times, the two processes dynamically interact throughout the time window and the aim is to optimize the objective function at its end.

One of the scenarios we considered is that of resource deployment in a mutually exclusive competitive case where the deployment of an agent’s resources aims to maximize its spread at the end of the time window with respect to the spread of an adversarial process. This can be seen as a competition between adversarial agents, each of which aims to contain the spread of the other, making it highly relevant in the battle for public opinion. The collaborative scenarios we optimized include constructive interaction between the spreading processes, such that infection by one process facilitates the infection by another. The aim in this case is to exploit the exposure optimally in order to maximize the spread of both processes. Finally, we also examined optimal vaccination strategies in the case of collaborative spreading processes, where vaccination is effective with respect to one process only but helps to reduce the spread of both due to the interaction between them.

The optimization algorithm has been tested in several different scenarios and on a range of small-scale and real networks, showing excellent performance with respect to the existing alternatives of high-degree allocation, K-shell based approaches, free spreading and uniform budget allocation. We demonstrated that the optimized deployment makes a huge difference in the use of limited resources and allows for a balanced competition even with significantly inferior resource availability or lower infection rates. We see a significant potential in the use of such principled algorithms to make the most of limited resource.

Moreover, the suggested framework is highly adaptive and can accommodate: targeted spreading [56], where only specific vertices are available and at specific times (as in the case of critical vote and time-sensitive campaigns); temporal deployment of resources, either of the process of interest or of the adversarial process, where resources are available at different times within the time window (e.g., due to limited production or shipment restrictions); determining vaccination policies and for more complex collaborative scenarios. \textcolorblackWe anticipate that the developed methods will also be useful for development of scalable [91] algorithms for learning [55, 51] of spreading models in the context of multiple interacting spreading processes.

The application of our method \textcolorblackfor suggesting vaccination policy to contain the spread of collaborative health conditions, for process containment (e.g., the spread of fake news, anti-vaxxing messages), and for effective information dissemination (e.g., marketing or opinion-setting material), are promising and timely applications that should be explored.

Acknowledgements.

\textcolorblackWe would like to thank Ho Fai Po for compiling the UK road network data. H.S. and D.S acknowledge support by the Leverhulme trust (RPG-2018-092) and the EPSRC Programme Grant TRANSNET (EP/R035342/1). A.Y.L. acknowledges support from the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20200121ER.

Appendix A Dynamic Belief Propagation

\textcolor

black In the two-process dynamics, the dynamics of a single node $i$ is described by a pair of activation times, $(\tau_{i}^{A},\tau_{i}^{B})$, where $\tau_{i}^{A}$ denotes the first time when node $i$ is found in the \textcolorblackstatus $A$, and similarly for $\tau_{i}^{B}$. For instance, $\tau_{i}^{A}=0$ means that node $i$ was initially in the active \textcolorblackstatus $A$, and we will denote by $\tau_{i}^{A}=*$ the situation where node $i$ did not get $A$-activated before some final observation time, i.e.\textcolorblack, in some sense $*$ absorbs all the history that happens after the end of the observation window. For the convenience of presentation, in what follows we will consider two separate “observation windows”, for $A$ and $B$ processes.

\textcolor

black The starting point for our derivations are the general Dynamic Belief Propagation (DBP) equations [40, 2, 54] on the interaction graph, where the goal is to approximate the probability $m^{i}_{T_{A}+1,T_{B}+1}(\tau_{i}^{A},\tau_{i}^{B})$ that node $i$ has a trajectory $(\tau_{i}^{A},\tau_{i}^{B})$ during the observation time window of length $T_{A}$ for process $A$ and $T_{B}$ for process $B$:

[TABLE]

where $\partial i$ denotes neighbors of $i$ on the graph, and the sum runs over all possible values of $\tau_{k}^{A}$ and $\tau_{k}^{B}$, i.e.\textcolorblack, $\{0,1,\ldots,t_{A/B},*\}$. The transition kernel $W^{i}$ is defined through the dynamical rules for a given dynamics; we will provide an explicit expression for the kernel below. The probabilities $m_{T_{A},T_{B}}^{i\to j}(\tau_{i}^{A},\tau_{i}^{B}\|\tau_{j}^{A},\tau_{j}^{B})$ have a sense of the probability that node $i$ has a trajectory of length $T_{A}$ for process $A$ and $T_{B}$ for process $B$ (in what follows referred to as $(T_{A},T_{B})$-trajectory) equal to $(\tau_{i}^{A},\tau_{i}^{B})$ subject to a fixed trajectory $(\tau_{j}^{A},\tau_{j}^{B})$ of node $j$, and satisfy the following consistency relations:

[TABLE]

The fixed point solution of equations (32)-(33) is guaranteed to be exact on tree graphs, and provides good estimates of the marginal probabilities on loopy but sparse graphs [60, 54].

\textcolor

black Notice that $m^{i\to j}_{T_{A},T_{B}}(\tau_{i}^{A},\tau_{i}^{B}\|\tau_{j}^{A},\tau_{j}^{B})$ is not a conditional probability (e.g.\textcolorblack, Bayes rule does not apply). Let us state the following fact that will be crucial for the derivations below. (This fact is due to the properties of the DBP equations and the dynamics specified by $W^{i}$.)

Fact 1** (Fundamental property of messages)**

For $\tau_{i}^{A}\neq*$ and $\tau_{i}^{B}\neq*$,

[TABLE]

A similar statement holds for messages.

\textcolor

black\textcolorblackThis means that future events cannot affect current or past events. Due to Fact 1, the DBP equations can be considerably simplified. Indeed, Eqs. (32) and (33) basically state that the right hand sides do not depend on the final times in the observation windows, i.e.\textcolorblack, on times $(T_{A}+1,T_{B}+1)$. In particular, this means that for $\tau_{i}^{A}=0$ or $\tau_{i}^{B}=0$, the right hand side does not depend on $\tau_{i}^{A}$ or $\tau_{i}^{B}$, respectively. For other cases, choosing $T_{A}=\tau_{i}^{A}-1$ and $T_{B}=\tau_{i}^{B}-1$, we see that since the time window of messages are $(\tau_{i}^{A}-1,\tau_{i}^{B}-1)$, they do not depend on the precise value of $\tau_{i}^{A}$ and $\tau_{i}^{B}$. Because of that, the following fact holds.

\textcolor

black

Fact 2** (Simplified DBP equations)**

For the dynamics considered, the following DBP equations hold:

[TABLE]

[TABLE]

To simplify the notations, in what follows we will use the notation $m^{k\to i}_{T_{A},T_{B}}(\tau_{k}^{A},\tau_{k}^{B})\equiv m^{k\to i}_{T_{A},T_{B}}(\tau_{k}^{A},\tau_{k}^{B}\|*,*)$. Given the computed value of the marginals $m_{T}^{i}(\tau_{i}^{A},\tau_{i}^{B})$, one can straightforwardly define quantities of interest, such as probabilities for a given node $i$ to be found in a given \textcolorblackstatus:

[TABLE]

In the last expression, we used Fact 1, and in the previous expressions, we used the definition of $*$ and hence the equalities of the type

[TABLE]

\textcolor

black Hence, in principle, to each node $i$ one can associate a $(t+1)\times(t+1)$ matrix of probabilities $m_{t,t}^{i}(\tau_{i}^{A},\tau_{i}^{B})$, and to each edge $(ij)$ associate a \textcolorblackfourth order tensor of probabilities $m_{t,t}^{i\to j}(\tau_{i}^{A},\tau_{i}^{B}\|\tau_{j}^{A},\tau_{j}^{B})$; solve the fixed point equations (35), e.g., through iteration using some initialization for the conditional probabilities; use (34) to compute the quantities $m_{t,t}^{i}(\tau_{i}^{A},\tau_{i}^{B})$; and finally, use (36)-(39) to assemble the probabilities of finding node $i$ in one of the \textcolorblackstatuses $S$, $A$, $B$, or $AB$. However, the naive implementation of this scheme would require as many as $O(t^{2(d-1)})$ operations for an estimation of a single message in (35) and $O(t^{2d})$ operations for a single marginal (34), which still produces a polynomial-time algorithm for bounded-degree graphs, but can quickly become prohibitive for large $d$.

\textcolor

black On the other hand, the argument above applies to general dynamics parametrized by the flipping times $\tau_{i}^{A},\tau_{i}^{B}$ for node $i$, while dynamics of interest (such as processes considered here) often has a special structure, and it is beneficial from the computational point of view to explore this structure in order to drastically reduce the complexity of the computation of the marginal probabilities (36)-(39). Therefore, while the problem is conceptually solvable in polynomial time at the level of the basic DBP equations (34)-(35), below we take care of deriving a computational scheme that would allow for an efficient computation of marginal probabilities (36)-(39). We refer to this scheme as to the Dynamic Message-Passing (DMP) equations.

Appendix B Exact DBP Equations for Mutually Exclusive Competitive Processes

\textcolor

black As explained in the Main Text, the transition rules for competitive dynamics at time step $t$ are defined as follows:

[TABLE]

where

[TABLE]

and $\sigma_{i}(t)$ denotes the \textcolorblackstatus of node $i$ at time $t$. For two subsets $\Theta_{A}\subset\partial i$ and $\Theta_{B}\subset\partial i$, \textcolorblackthe transmission probabilities introduced \textcolorblackin (41)-(43) can be explicitly written as follows:

[TABLE]

\textcolor

black Given the usual parametrization of marginals and messages $m^{i}(\tau_{i}^{A},\tau_{i}^{B})$, only the following messages can be non-zero for this dynamics: $m^{i}(\tau_{i}^{A},*)$, $m^{i}(*,\tau_{i}^{B})$, and $m^{i}(*,*)$. The dynamic messages are then defined as follows:

[TABLE]

with the initial conditions $m^{i}_{t,t}(0,*)=P^{i}_{A}(0)$ and $m^{i}_{t,t}(*,0)=P^{i}_{B}(0)$.

\textcolor

black As discussed in Appendix A, the complexity of a general DBP message evaluation requires $O(t^{2d})$ steps, where $d$ is the node degree. Under the mutually exclusive scenario, some of the combinations of flipping times are forbidden since node can not be simultaneously activated under both $A$ and $B$ processes. To reflect this structure, we rewrite DBP equations for $t>0$ in a computationally sub-optimal, but compact way:

[TABLE]

DBP messages in the expressions above are computed in a similar way, where subsets $\Theta_{A}$ and $\Theta_{B}$ are chosen from $\partial i\backslash j$:

[TABLE]

In Section II of the Main Text, we show numerically that these equations are exact on tree graphs.

Appendix C Approximate DMP Equations for Mutually Exclusive Competitive Processes

\textcolor

black In the search of low-complexity equations that could be used on general graphs, in this Section we derive approximate DMP equations for mutually exclusive competitive processes. As outlined in the Main Text, our approximation scheme is based on the idea that in the absence of renormalization, the dynamics of each processes follows the dynamics of a simple SI process. Hence, we derive equations that perform \textcolorblackthe renormalization procedure at the level of dynamic marginals.

\textcolor

black Let us start by reminding the equations for the SI model. We will closely follow the notations used in [53, 54]. For a single process $A$, let $\theta_{A}^{i\rightarrow j}(t)$ denote the probability that no infection message $A$ has been passed from node $i$ to node $j$ until time $t$, and $\phi_{A}^{i\rightarrow j}(t)$ denote the probability that no infection message $A$ has been passed until time $t\!-\!1$, while node $i$ is infected with $A$ at time $t$. Exact DMP equations for a single-process SI model read [54]:

[TABLE]

\textcolor

black For the two competitive and mutually exclusive processes $A$ and $B$, we approximate the total spreading as a product of two independent spreading processes that follow equations (59)-(64). Our goal will be to derive approximate dynamic marginals, that we denote by $\widehat{P}_{S}^{i}(t)$, $\widehat{P}_{A}^{i}(t)$, and $\widehat{P}_{B}^{i}(t)$. For simplicity, let us assume that $\alpha^{A}_{ij}<1$, and hence $\theta_{A}^{i\to j}$ type messages are non-zero (treatment of the case of deterministic spreading is straightforward and results in several additional equations). Under the chosen approximation, the non-normalized probability of staying in the susceptible \textcolorblackstatus can be written as:

[TABLE]

\textcolor

black Under this approximation, the transition to \textcolorblackstatus $A$ happens when the target node is susceptible, and the $B$ infection is not transmitted:

[TABLE]

We then compute the renormalized messages for the three \textcolorblackstatuses $\widehat{P}_{S/A/B}^{i\rightarrow j}(t)$:

[TABLE]

\textcolor

black The dynamic marginals are calculated iteratively in a similar fashion:

[TABLE]

\textcolor

blackWe numerically study the performance of these approximate DMP equations in the Main Text.

Appendix D Exact DBP and DMP Equations for Collaborative Processes

D.1 DBP Equations for Collaborative Processes

\textcolor

black Let us start by reminding the transition rules for the collaborative spreading, as discussed in the Main Text:

[TABLE]

\textcolor

black We start by specifying the transition kernel $W^{i}(\tau_{i}^{A},\tau_{i}^{B};\{\tau_{k}^{A},\tau_{k}^{B}\}_{k\in\partial i})$ based on the dynamic rules (69)-(72):

[TABLE]

In the expression above, it is assumed that when $\tau_{i}^{A}=*$ or $\tau_{i}^{B}=*$, the corresponding transition does not happen during the given observation window; or, alternatively, happens eventually, beyond the observation horizon. $P^{i}_{S}(0)$, $P^{i}_{A^{*}}(0)$, $P^{i}_{B^{*}}(0)$, and $P^{i}_{AB}(0)$ denote probabilities that at initial time, node $i$ is initialized in the \textcolorblackstatuses $S$, $A$ exclusively, $B$ exclusively, or $AB$, respectively.

\textcolor

black We can now provide explicit forms of the expressions that correspond to the transitions inside curly brackets\textcolorblack; specifically, in the following we refer to the numbered items in (73)-(78). In the case $\tau_{i}^{A}=*$ and $\tau_{i}^{B}=*$ we have:

[TABLE]

[TABLE]

\textcolor

black If either $\tau_{i}^{A}=*$ or $\tau_{i}^{B}=*$, the activation above is replaced by an absence of activation. For instance, in the case \textcolorblackof $\tau_{i}^{B}=*$, the last term in \textcolorblackEq. (79) would read:

[TABLE]

for the $(t+1,t+1)$-trajectory of node $i$.

\textcolor

black The DBP equations for collaborative processes are given by (34)-(35) using the $W^{i}(\tau_{i}^{A},\tau_{i}^{B};\{\tau_{k}^{A},\tau_{k}^{B}\}_{k\in\partial i})$ defined above. In what follows, we explain how \textcolorblackto derive the equivalent low-complexity DMP equations for the collaborative processes from the DBP equations (34)-(35), state their full iterative form, and finally expand the exact DMP equations around the non-interacting point, which yields approximate DMP equations that have a more compact form and a lower algorithmic complexity.

D.2 Definitions of Dynamic Messages

\textcolor

black Let us introduce the mathematical definitions for the dynamic messages as the functions of \textcolorblackfundamental probabilities on time trajectories, $m_{T_{A},T_{B}}^{i\to j}(\tau_{i}^{A},\tau_{i}^{B})\equiv m_{T_{A},T_{B}}^{i\to j}(\tau_{i}^{A},\tau_{i}^{B}\|*,*)$, that will be used in the computations scheme below. First, it is useful to consider reduced marginals:

[TABLE]

They have the physical meaning of probabilities of individual activation with one of the processes, when we do not care about the activation by the other process. In a similar way, we also define reduced messages:

[TABLE]

\textcolor

black Similarly to dynamic marginals (36)-(39), one can define the aggregated dynamic messages. In principle, we can defined them for a general fixed dynamics $(\tau_{j}^{A},\tau_{j}^{B})$ of the “cavity” node $j$, but in what follows we will only use the messages for the fixed dynamics $(\tau_{j}^{A},\tau_{j}^{B})=(*,*)$.

[TABLE]

These quantities have the same physical interpretation, except that they are defined on the “cavity” graph where node $j$ is following a fixed $(\tau_{j}^{A},\tau_{j}^{B})$-dynamics.

\textcolor

black Let us also define the following dynamic messages that are weighted combinations of messages:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These messages have the meaning of the probabilities that node $k$ did not send $A$ and $B$ activation signals to node $i$ on the “cavity” graph where node $i$ is following a fixed $(*,*)$-dynamics.

\textcolor

black Next, we define the following quantities:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These messages have the meaning of the probabilities that node $k$ did not send $A$ and $B$ activation signals to node $i$ on the “cavity” graph where node $i$ is following a fixed $(*,*)$-dynamics, but node $k$ is in an active \textcolorblackstatus itself.

\textcolor

black Finally, we define auxiliary expressions:

[TABLE]

[TABLE]

D.3 Update Equations for Dynamic Messages

\textcolor

black We start with the definition

[TABLE]

From the DBP equations on time trajectories (34), we get

[TABLE]

where

[TABLE]

Using the definition (94)

[TABLE]

For simplicity, let us denote $\theta_{A,B}^{k\rightarrow i}(t+1,t+1)\equiv\theta_{A,B}^{k\rightarrow i}(t+1,t+1\|*,*)$. With this simplified notation, we finally write

[TABLE]

We will use the same simplified notation for all quantities of the type $\#^{k\to i}(...)\equiv\#^{k\to i}(...\|*,*)$. Using the previously defined quantities $\phi_{A}$, $\phi_{B}$ and $\phi_{A,B}$, through (99), (100) and (101), we derive the following update equations:

[TABLE]

\textcolor

black Using the equality $\mathbbm{1}[\tau_{k}^{A}\leq t_{A}]=\mathbbm{1}[\tau_{k}^{A}\leq t_{A}-1]+\mathbbm{1}[\tau_{k}^{A}=t_{A}]$, we further get

[TABLE]

where $\theta_{A}$ and $\theta_{B}$ have been defined above. Finally, using

[TABLE]

we get in the end:

[TABLE]

[TABLE]

Additionally,

[TABLE]

where $m_{t,t}^{k\rightarrow i}(\tau_{k}^{A},\tau_{k}^{B})\equiv m_{t,t}^{k\rightarrow i}(\tau_{k}^{A},\tau_{k}^{B}\|*,*)$. Finally, we need to compute iteratively $\theta_{A}^{k\rightarrow i}(t,t)$ and $\theta_{B}^{k\rightarrow i}(t,t)$:

[TABLE]

where by definition $\theta^{k\to i}_{A}(0,t)=\mu^{k\to i}_{B}(t)=P^{k\to i}_{B}(t)-P^{k\to i}_{B}(t-1)$, and similarly $\theta^{k\to i}_{B}(t,0)=\mu^{k\to i}_{A}(t)=P^{k\to i}_{A}(t)-P^{k\to i}_{A}(t-1)$. Initial and border conditions of that type can be directly obtained from the definitions of the dynamic messages and their connections to messages on time trajectories.

\textcolor

black From the scheme above, we still need to compute the update equations for the marginals and the messages on time trajectories of the type $m_{t,t}^{i}(\tau_{i}^{A},\tau_{i}^{B})$ and $m_{t,t}^{i\rightarrow j}(\tau_{i}^{A},\tau_{i}^{B})$. The derivation of the equations for these messages is given below. Once the full or reduced marginals on time trajectories are computed, the dynamic marginals of interest are easy to compute according to the following equations:

[TABLE]

[TABLE]

[TABLE]

D.4 Computation of $m_{t,t}^{i}(\tau_{i}^{A},\tau_{i}^{B})$ and $m_{t,t}^{i\rightarrow j}(\tau_{i}^{A},\tau_{i}^{B})$

\textcolor

black It is convenient to break down the computation depending on a particular combination of $(\tau_{i}^{A},\tau_{i}^{B})$, including cases of times [math] and $*$. We will provide a detailed derivation of one of the most general cases, all the others being can be treated similarly. Let us consider the case of finite-time $(\tau_{i}^{A},\tau_{i}^{B})$. Several cases of the type depicted in Fig. 9 need to be considered. For the sake of illustration, we will consider the subcase $*\neq\tau_{i}^{B}>\tau_{i}^{A}>0$.

\textcolor

black Using Fact 35, we the simplified DBP equations read

[TABLE]

and

[TABLE]

Using the definitions for the dynamics in the case considered, i.e.\textcolorblack, term (1) in $W^{i}(\tau_{i}^{A},\tau_{i}^{B};\{\tau_{k}^{A},\tau_{k}^{B}\}_{k\in\partial i})$, as well as for the dynamic messages (95), we obtain

[TABLE]

and, in the same way,

[TABLE]

Through the definition, we immediately get

[TABLE]

where the former quantity is defined through (102). Working with this definition, we get the following update equation for $\phi_{A,B,AB}^{k\to i}(t_{A},t_{B},t_{AB};\tau_{i}^{A})$:

[TABLE]

The term $\prod_{t^{\prime\prime}=0}^{t_{B}-1}(1-\alpha_{ki}^{B}\mathbbm{1}[\tau_{k}^{B}\leq t^{\prime\prime}])$ is not active for the times considered here and is equal to one because $\tau_{k}^{B}=t_{AB}$ and $t_{AB}\geq t_{B}$, and hence we get the final update equation

[TABLE]

where all $m_{t_{AB},t_{AB}}^{k\to i}(t^{\prime},t^{\prime\prime})$ are known from previous step computations.

D.5 Exact DMP Equations for Collaborative Processes

\textcolor

black In this Section, we summarize the derivation procedure outlined above, and state the full form of the DMP equations for collaborative processes.

\textcolorblack Initial conditions from the problem definition:

[TABLE]

\textcolor

black Initialization of dynamic messages:

[TABLE]

\textcolor

black Border conditions for dynamic messages:

[TABLE]

\textcolor

black Expression for the dynamic marginals:

[TABLE]

\textcolor

black Update equations:

[TABLE]

\textcolor

black

[TABLE]

\textcolor

black The update expressions for the messages $m_{t,t}^{i\to j}(\tau_{i}^{A},\tau_{i}^{B})$ are similar to the expressions for the marginals $m_{t,t}^{i}(\tau_{i}^{A},\tau_{i}^{B})$ above, except $\partial i$ is replaced with $\partial i\backslash j$:

\textcolor

black

[TABLE]

\textcolor

black The dynamic messages entering the expressions above are updated as follows.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

\textcolor

black Order of the updates:

Only once, at time zero: initialize messages 2. 2.

If initialization is not given for a message, it means that its initial value is fixed through a border condition (reduction to a different message), and will need to be initialized in this way for each update 3. 3.

Update $\phi$’s 4. 4.

Update $\theta$’s 5. 5.

Update $m$ messages and marginals 6. 6.

Update $\mu$ messages and marginals 7. 7.

Compute dynamic marginals $P_{S}^{i}(t),P_{A}^{i}(t),P_{B}^{i}(t),P_{AB}^{i}(t)$.

Appendix E Approximate DMP Equations for Collaborative Processes

\textcolor

black From the dynamic rules for the interacting spreading processes (69)-(72), we notice that the special case of transmission probabilities $\alpha_{ij}^{AB}=\alpha_{ij}^{A}$ and $\alpha_{ij}^{BA}=\alpha_{ij}^{B}$ corresponds to non-interacting spreading processes: activation by one process does not change the activation dynamics for the other. In this case, the DMP equations should simplify into the product of two independent SI-like processes:

[TABLE]

We use this observation to produce a simplified version of DMP equations, expanding exact equations for interacting spreading processes around the non-interacting point. Several expansions are possible, including the most straightforward one where we expand to a certain order in $\alpha_{ij}^{A}-\alpha_{ij}^{AB}$ and $\alpha_{ij}^{B}-\alpha_{ij}^{BA}$, which will result in bulky expressions. Instead, here we produce a simplified set of DMP equations by keeping certain first-order corrections in the update equations only, so that the resulting equations are similar to the SI-type equations. Using the approximation

[TABLE]

and the general expression

[TABLE]

we compute the corresponding contribution of each term. For simplicity, in what follows we assume that $\alpha_{ij}<1$ and hence $\theta^{k\to i}$ type messages are non-zero (treatment of the case of deterministic spreading is straightforward and results in several additional equations).

\textcolor

black Simplified approximate computation of $m_{t,t-1}^{i}(t,*)$: We start with the following approximations in the dynamic messages:

[TABLE]

[TABLE]

The update equation

[TABLE]

becomes under this approximation

[TABLE]

or simply

[TABLE]

Using these simplified expressions, we get after some algebra

[TABLE]

\textcolor

black Simplified approximate computation of $\sum_{t^{\prime}\leq t-1}m_{t,t-1}^{i}(t,t^{\prime})$: Similarly, we start with the approximations in the dynamic messages:

[TABLE]

The approximate version of the update equations

[TABLE]

then reads

[TABLE]

or simply

[TABLE]

Using these simplified expressions, we get after some algebra

[TABLE]

\textcolor

black Final form of the approximate DMP equations: The final update equations that we will need,

[TABLE]

with

[TABLE]

gives

[TABLE]

\textcolor

black Combining all of the computation above, we finally obtain the following approximate form of the DMP equations that should provide a good approximation around the non-interacting point:

[TABLE]

[TABLE]

[TABLE]

\textcolor

black In the Main Text, we numerically establish the approximation power of these equations. In what follows, we use these approximate DMP equations for performing optimization due to their simpler form and lower computational complexity.

Appendix F Optimization of Mutually Exclusive Competitive Processes

For competitive processes we study two optimization problems, the multi-agent seeding problem and disease containment, the only difference between the two is the objective function. For disease containment the objective function to be maximized is $\mathcal{O}\!=\!\sum_{i}\left(1-P_{i}^{B}(T)\right)$ and for the multi-agent seeding problem $\mathcal{O}\!=\!\sum_{i}(1-P_{i}^{S}(T))~{}.$ The budget constraint at time zero is enforced by the Lagrange multiplier $\lambda^{Bu}$ is:

[TABLE]

where permitted $\underline{\nu}<\nu^{i}<\overline{\nu}$ value restrictions are enforced by the term

[TABLE]

In this case, the restrictions used are $\overline{\nu}\!=\!1$ and $\underline{\nu}\!=\!0$.

Initial conditions are forced through a set of Lagrange multipliers $\lambda$ for the various parameter values

[TABLE]

The DMP equations for the dynamics are forced through a set of Lagrange multipliers

[TABLE]

Derivatives with respect to the dynamics parameters give rise to the optimization (dynamical) equations for the Lagrange multipliers. The case described here is that of the containment problem (minimizing spread of an adversarial process).

[TABLE]

[TABLE]

[TABLE]

The equations for ${\partial\mathcal{L}}/{\partial P_{B}^{i\rightarrow j}(t)}$, ${\partial\mathcal{L}}/{\partial P_{B}^{i}(t)}$, ${\partial\mathcal{L}}/{\partial\widehat{P}_{B}^{i\rightarrow j}(t)}$, ${\partial\mathcal{L}}/{\partial\widehat{P}_{B}^{i}(t)}$, ${\partial\mathcal{L}}/{\partial\theta_{B}^{i\rightarrow j}(t)}$, ${\partial\mathcal{L}}/{\partial\phi_{B}^{i\rightarrow j}(t)}$ are similar to their $A$ process counterparts with an exchange of variables $A\leftrightarrow B$.

For simplicity one can write equation (218) as:

[TABLE]

The same method as in [56] is used here to solve the quadratic equation (219), writing $\nu^{i}(0)$ as a function of $\psi_{i}$:

[TABLE]

In this scenario, the positive square root solution satisfies restriction (213). Given the budget restriction one can obtain $\nu^{i}(0)$ numerically through

[TABLE]

One can use the following update procedure to obtain optimal solution iteratively:

1. Step 1:

Carry out the forward iteration using given initial conditions, Eqs. (65)-(68) 2. Step 2:

Compute all Lagrangian multipliers at $t\!=\!T$. 3. Step 3:

Compute $\lambda_{ij}^{S}(t\!-\!1),\lambda_{i}^{S}(t\!-\!1),\lambda_{i}^{A}(t\!-\!1),\lambda_{i}^{B}(t\!-\!1),\lambda_{ij}^{\phi_{A}}(t\!-\!1),\lambda_{ij}^{\phi_{B}}(t\!-\!1)$ using obtained Lagrangian multipliers at time $t$ from Eqs.(216)-(217). 4. Step 4:

Compute $\lambda_{ij}^{A}(t\!-\!1),\lambda_{ij}^{B}(t\!-\!1),\lambda_{i}^{\widehat{A}}(t\!-\!1),\lambda_{i}^{\widehat{B}}(t\!-\!1),\lambda_{i}^{\widehat{S}}(t\!-\!1)$ using obtained Lagrangian multipliers and Eqs.(216)-(217). 5. Step 5:

Compute $\lambda_{ij}^{\widehat{A}}(t\!-\!1),\lambda_{ij}^{\widehat{B}}(t\!-\!1),\lambda_{ij}^{\widehat{S}}(t\!-\!1)$ using obtained Lagrangian multipliers and Eqs.(216)-(217). 6. Step 6:

Compute $\lambda_{ij}^{\theta_{A}}(t\!-\!1),\lambda_{ij}^{\theta_{B}}(t\!-\!1)$ using obtained Lagrangian multipliers. 7. Step 7:

Repeat Steps 1-5 back in time until all Lagrangian multipliers have been obtained for the range $t\!=\!0,\ldots,T$. 8. Step 8:

Solve Eqs. (220) and (221) numerically and update $\nu^{i}$. 9. Step 9:

Repeat Step 1 to Step 8 until convergence. 10. Step 10:

Calculate the posterior marginals on the basis of the messages using Eqs. (C).

A similar approach is used in collaborative processes using the corresponding equations for the dynamics and the Lagrange multipliers’ dynamics.

Appendix G Optimization of Competitive Processes - Example

We also demonstrate the efficacy of DMP-optimal algorithm on the denser network of the 1994 world metal trade network [71] .The network consists of 80 countries; we will ignore their respective trade volumes and use the dense topology for the current example. The infection probabilities used for the two processes are $\alpha_{A}\!=\!\alpha_{B}\!=\!0.5$; the budget for $B$ is 1, allocated at node 1 (Argentina) at $t\!=\!0$, and a budget of 1 per time step is assigned using the DMP-optimal strategy for process $A$ within the time window of $T\!=\!3$. The results are shown in Fig. 11, where the subfigures represent the containment of process $B$ at different times. The heat bar represents the dominating process per node through the value $P_{A}^{i}(t)\!-\!P_{B}^{i}(t)$. Red/blue represent dominating processes $A/B$, respectively.

Appendix H Optimization of Collaborative Processes

Optimization in collaborative processes is highly relevant in many health-related problems such as the co-epidemic of HIV and tuberculosis, and in the marketing of correlated products or opinions. In the co-epidemic case, identifying the most catastrophic health-risk scenario is an important aspect of evaluating the risk of a co-epidemic event; while in the marketing area, it aims at the best joint advertising campaign.

We consider two main scenarios in our model: (a) The multi-agent seeding problem, i.e., computing the optimal allocation of one agent which minimizes the number of susceptible vertices at a certain time $T$ for a given spread of the second agent. (b) Optimal allocation of vaccines against one of the processes for maximizing the impact on the spread of both processes. Assume we do not have any vaccines for HIV but vaccines for tuberculosis are available (such as Bacillus Calmette-Guérin [90]); the optimization problem would be articulated as: given a certain budget of vaccines for tuberculosis, what is the best vaccination policy to minimize the spreading of both?

The Lagrangian used is similar to the one of the competitive case but with different dynamics and initial conditions. The dynamics equations are enforced through a set of Lagrange multipliers $\lambda$ as before

[TABLE]

The initial conditions in this case are of the form:

[TABLE]

Appendix I Optimization of Vaccine Allocation

The multi-process seeding problem in the collaborative case has a similar structure to that of the containment task in competitive processes. The vaccine allocation problem is of a slightly different nature. There are many ways in which it can be formulated, the model illustrated below is the one we have chosen to use.

If a node receives a unit of vaccine before being exposed to infected neighbors, it will be immune to that process (denoted as $B$ in our model), which means the infection message from its neighbor will decrease to 0. However, sometime a node receives part of a unit of vaccine; specifically, if one allocates a certain amount of vaccine, say $b$ to a node, the probability of it being infected will decrease to $\alpha_{B}\!-\!b$ (bounded by 0 from below), where $\alpha_{B}$ is the initial infection parameter. Meanwhile, the parameter $\alpha_{BA}$ will also decrease to $\alpha_{BA}\!-\!b$. Therefore the budget for vaccines can be expressed by the budget for decreasing $\alpha_{B}$ and $\alpha_{BA}$ in Eqs. (205)-(211), which now take the following form:

[TABLE]

The budget restriction in this case is:

[TABLE]

where $Bu$ is the total vaccination budget. Another restriction for $b(i)$ allocated on a node (not exceeding the infection probability) is expressed similarly to before as:

[TABLE]

Finally, the derivative of the part of the Lagrangian that enforces the dynamics $\mathcal{D}$ is differentiated with respect to $b(i)$, required to complete the set of equations provides

[TABLE]

The same procedure (following after Eq.(219)) can be implemented for the optimization in this case.