INTRODUCTION

During the COVID-19 pandemic, humanity faced not only a problem of a spread of the pandemic itself, but also an unexpected phenomenon of formation of negative attitude towards vaccination in the society. Since this latter process does significantly affect development and resulting consequences of the pandemic, it is very important to build qualitative understanding of factors influencing formation of public opinion in the particular context of epidemic dynamics.

To develop a quantitative model describing interrelated evolution of epidemics and epidemics-related attitude it is natural to use a two-layered network with interconnected strata containing networks on which there takes place formation of epidemics and epidemics-related attitude/opinion respectively. The particular model we build combines an Ising-type mechanism of formation of binary yes/no-type opinion operating on the social contact network coupled to the SIRS (Susceptible-Infected-Recovered-Susceptible) epidemic mechanism (see [1, 2] and references therein) operating on the network of physical contacts in real space. The model does thus describe two interrelated spreading phenomena taking place in the two-layer network in which two networks are glued at nodes representing agents.

Behavioural motives behind epidemics were previously analysed in game-theoretic setting in the literature on vaccination decision formation (for detailed review see [3] and references therein). However, to the best of authors’ knowledge, these models mostly analyse the vaccination-related free-rider type problems. The COVID 19 pandemic showed that at the early stages of viral spreading vaccination decisions are to a large extent based not on a proper factual analysis of contagion probability with or without vaccination, but mostly on subjective analysis of opinion of acquaintances, information from media, etc., as factual data on vaccination efficiency, virus spreading and its danger are not yet available or seem to be not enough reliable. One of possible ways of describing public opinion formation in the binary choice setting is to use an Ising-type model [411]. The present paper introduces an Ising-type model describing vaccination-related decisions.

The main results of the paper are the following.

First, we uncover a nontrivial dependence between the noise amplitude influencing the vaccination opinion formation and the fraction of recovered agents (i.e., people who caught the virus and later recovered from it). It turns out that an increase of the noise amplitude triggers two oppositely directed effects. On the one hand, public opinion formation on vaccination takes more time and, therefore, epidemics spreading is faster. On the other hand, the probability of reaching an equilibrium where almost all individuals are vaccinated and, therefore, epidemics terminates, is higher. A superposition of these oppositely directed trends leads, in principle, to an existence of some “optimal” noise amplitude. The position of this optimum is, however, strongly dependent on the values of model parameters.

Second, we uncover a new effect related to the external field influence. In the model under consideration the external field represents an intensity of the official information delivery. We show that if initial beliefs on vaccination and agents’ attitude to official information sources are highly correlated and clusterized, growing intensity of the official information delivery intensity can lead to a decrease of probability of reaching a preferable equilibrium with respect to vaccination attitude and, consequently, to potential growth of epidemics—related risks.

NETWORK

The bilayer network under consideration comprises the network representing physical contacts and that of social contacts representing existing communication channels. Both constituent networks have common nodes representing agents. The first layer is a scene of the SIRS epidemic spreading process. The second (social communication) layer describes communication interactions influencing vaccination decisions.

Let us start with the description of the first layer. It is well known that a topology of a contact network is one of the key determinants of the properties of e-pidemic spreading process. Therefore, to construct a  realistic description of epidemic dynamics it is n-atural to use a real contact network as a basic building block for constructing the contact layer of our bilayer graph. The particular real contact network we chose is the Contact high-school dataset from https://www.cs.cornell.edu/arb/data/contact-high-school/ [12, 13]. This data contains information about interactions recorded by sensors worn by high school students. The original dataset contains information about simplexes, i.e., sets of nodes representing groups of students interacting at some particular moment so that it describes a contact hypergraph. The data is recorded with a time resolution of 20 s. In our computation we use a weighted clique projection of this hypergraph where weights on a link connecting two nodes are equal to number of simplexes in which they appear together. From this weighted clique projection we take only links with weights higher than 20 (i.e. with total time of interaction exceeding 7 min during one week). Using these links we constructed an unweighted graph \({{G}_{0}} = ({{V}_{0}},{{E}_{0}}),\) where \({{V}_{0}}\) is the set of nodes represented students and \({{E}_{0}}\) is the set of edges modelling real-space contacts between students. The number of students participating in the research described in [12] was equal to 327. The network \({{G}_{0}}\) contains only connected nodes (this is ensured by preserving only the links with weights higher than 20). For the constructed network \({{G}_{0}}\) we have \({\text{|}}{{V}_{0}}{\text{|}} = 320\) and \({\text{|}}{{E}_{0}}{\text{|}} = 1174\).

This network \({{G}_{0}}\) is, however, too small for building a realistic model of epidemic spreading. To construct a larger network with realistic properties one can use two approaches. The first one is to estimate key properties of the small graph and reproduce them in the random graph of acceptable size. The second approach is to serialize \({{G}_{0}}\) and add some random links to obtain a connected graph. Although the first method looks more reliable and natural for the modelling purposes, in the particular case of our graph \({{G}_{0}}\) it is hardly applicable, as this graph contains only 320 nodes. In such a small graph it is difficult to identify such network mesoscopic structures as communities of connected nodes. It is well known that real complex networks are characterized by the existence of tightly connected communities of nodes which are interconnected through the so-called weak ties, bridges or local bridges (see [14] Section 3 and [2] Section 14). This property is known to have strong impact on epidemic dynamics [1517] so that a synthetic network should posess realistic community structure. One of the methods of graph randomisation preserving community structure and degree distribution is proposed in [17]. In this paper authors compare characteristics of epidemic spread in the standard configuration model and in two models preserving community structure of real networks. In the first model internal links of real communities are taken as given (experimental) and randomization is performed only for inter-community links, whereas in the second model intra-community links are also randomized. The results of this analysis show that evolution of epidemics in the classical configuration model (without community structure) is significantly different from that in the two cases of random graphs with community structure. At the same time, the dynamics in the models with community structure is approximately the same. Considering these results, in our study we suppose that the contact network \({{G}^{{{\text{cnt}}}}}\) comprises a set of communities each having the \({{G}_{0}}\) topology. As the epidemic dynamics including internal randomisation and without it is approximately the same we choose the simpler option without internal randomisation. However, in the future it is of interest to provide a robustness check using the model with intra-communities randomisation. Formally, for the modelling purposes we construct a synthetic contact network inheriting the main features of the original network \({{G}_{0}}\) and having a natural community structure as follows:

(1) we take \({{n}^{{{\text{copies}}}}}\) of \({{G}_{0}}\)\(\{ G_{0}^{{(1)}}, \ldots ,G_{0}^{{({{n}^{{{\text{copies}}}}})}}\} \);

(2) we combine \(\{ G_{0}^{{(1)}}, \ldots ,G_{0}^{{({{n}^{{{\text{copies}}}}})}}\} \) into a single disconnected graph \(G_{0}^{{{\text{cnt}}}}\) with \({{n}^{{{\text{copies}}}}}\) similar disconnected components isomorphic to \({{G}_{0}}\). In this graph we therefore have \({{n}^{{{\text{copies}}}}}{\text{|}}{{V}_{0}}{\text{|}}\) nodes and \({{n}^{{{\text{copies}}}}}{\text{|}}{{E}_{0}}{\text{|}}\) edges;

(3) we randomly add 0.5% of new edges (\(0.005{{n}^{{{\text{copies}}}}}{\text{|}}{{E}_{0}}{\text{|}}\)) into \(G_{0}^{{{\text{cnt}}}}\) obtaining in such a way connected contact network \({{G}^{{{\text{cnt}}}}}\).

In this study we use \({{n}^{{{\text{copies}}}}} = 4\). Therefore, we work with contact network consisting of the four tightly connected communities and few random bridges between them.

We now turn to the description of the layer containing the social communication network. It is common knowledge that Web-based social networks drastically simplify communication between people. Thanks to different modern communication methods people can now share their opinions being even in different continents. This argument shows that social communication network includes much more links than that of physical contacts. In particular, in the network of social contacts we expect to have all the links that exist in the physical contact layer, as in the web-based social networks (such as Facebook, Twitter, etc.) people typically add their “physical familiars” into the friends/followees/etc. set. However, in the web-social networks people can easily follow many other individuals. It is important that a distribution of these additional (with respect to the physical contact layer) between the nodes is not uniform. It is well known that such social-media are characterized by an existence of highly influential nodes (hubs) having enormously large degrees as compared to those of most of other nodes. The simplest way to reproduce this property is to use the preferential attachment mechanism [18]. Formally, we construct the social communication layer network \({{G}^{{{\text{scl}}}}}\) from \({{G}^{{{\text{cnt}}}}}\) as follows:

(1) we introduce the parameter \({{n}^{{{\text{add}}}}}\) representing the number of additional links for each node;

(2) for each node we construct \({{n}^{{{\text{add}}}}}\) of additional random links through preferential attachment with other nodes.

In the study we use \({{n}^{{{\text{add}}}}} = 10\).

DYNAMICS

The proposed model describes an interrelated dynamics of infection and vaccination and attitude towards it developing on top of the two networks, those of real space physical contacts and social interaction correspondingly. Each node/agent i participates in both of these dynamical processes. Its state at each given time t is characterised by the set \((\Xi _{t}^{i},v_{t}^{i},b_{t}^{i})\), where \(\Xi \) corresponds to one of the states S (susceptible), I (infected) or R (recovered (released)), \(v = \pm 1\) corresponds to a node being vaccinated/not vaccinated and \(b = \pm 1\) is an indicator of positive/negative attitude towards vaccination so that \(v_{t}^{i} = + 1\) if the node i is vaccinated at time t and \(v_{t}^{i} = - 1\) if it is not while \(b_{t}^{i} = + 1\) if the agent i stands for vaccination and \(b_{t}^{i} = - 1\) otherwise. The connection between the dynamic processes appearing in the two strata is thorugh the probability to catch the virus: once an agent gets the vaccine (\(v_{t}^{i} = + 1\)) the probability to become infected drops. Therefore, the opinion dynamics directly influences the epidemics spread.

Let us first describe the SIRS epidemic spread driven by transitions \(S \to I \to R \to S\) over the physical contact network layer. We start from putting at time \({{t}^{{{\text{start}}}}}\) one random node into the state I (infected). At each following discrete time step \(t \in [{{t}^{{{\text{start}}}}},T]\) we proceed as follows.Footnote 1

1. For each susceptible node \(i\) we calculate the probability \(p_{{i,t}}^{{\text{I}}}\) of becoming infected. This probability is determined by the two factors: (1) the number of infected node’s neighbours \(n_{{i,t}}^{{\text{I}}}\) and (2) the individual’s status with respect to vaccination. We have

$$p_{{i,t}}^{{\text{I}}} = 1 - {{(1 - p)}^{{n_{{i,t}}^{{\text{I}}}}}},$$
(1)

where p is the probability p of being infected due to interaction with an infected neighbouring node defined by

$$p = {{p}^{{nv}}}{{\delta }_{{v_{t}^{i}, - 1}}} + {{p}^{v}}{{\delta }_{{v_{t}^{i},1}}},\quad {{p}^{{nv}}} > {{p}^{v}}.$$
(2)

The values of \({{p}^{v}}\) and \({{p}^{{nv}}}\) are model parameters.

2. Each infected node stays in the state I for \({{\tau }^{{\text{I}}}}\) time periods. After that the node changes its state to R (recovered (released)).

3. A node stays in the state R for \({{\tau }^{R}}\) time periods. After that the node changes its state to S (susceptible).

Simultaneously with the SIRS epidemic in the social network layer there takes place a dynamical process of vaccination opinion formation and vaccination as such. At each time point \(t \in [0,T]\) the state \({\mathbf{s}}_{t}^{i}\) of the node \(i\)’s is \({\mathbf{s}}_{t}^{i} = (b_{t}^{i},v_{t}^{i})\). We suppose that:

• vaccination is irreversible; i.e., if the agent i gets the vaccine at time \(t = {{\tau }_{i}}\), then \(v_{t}^{i} = - 1,\;\forall t < {{\tau }_{i}}\) and \(v_{t}^{i} = + 1,\;\forall t \geqslant {{\tau }_{i}}\);

• once an agent is vaccinated, he always supportive of vaccination later on; i.e., \(b_{t}^{i} = + 1,\;\forall t \geqslant {{\tau }_{i}}\). Therefore, starting from the vaccination time point agents act as “positive zealots” [1922]; i.e., inflexible agents who can’t change their opinions during the game.

We also assume that an agent decides to get vaccinated if he has been supportive of vaccination during some preceding time period (memory size) θ; i.e.,

$$t = {{\tau }_{i}} \Leftrightarrow b_{{t - \theta + 1}}^{i} = b_{{t - \theta + 2}}^{i} = \ldots = b_{t}^{i} = + 1.$$
(3)

The process of formation of opinion on vaccination is driven by transitions \(b_{t}^{i} \to - b_{t}^{i}\) described in the framework of noisy binary choice game characterised by the choice utility of the form

$${{U}_{i}}(b_{t}^{i},t) = {{\chi }_{i}}Hb_{t}^{i} + J\sum\limits_{j \ne i} a_{{ij}}^{{{\text{scl}}}}b_{t}^{i}b_{{t - 1}}^{j} + {{\varepsilon }_{{b_{t}^{i}}}},$$
(4)

where \(H > 0\) describes universal external influence, for example, appeals for vaccination in media, \({{\chi }_{i}} = \pm 1\) parametrises the agent i’s attitude towards this influence, J is the parameter determining the magnitude of influence of the agent i’s neighbours’ beliefs, \(a_{{ij}}^{{{\text{scl}}}}\) is the social communication graph adjacency matrix element and \({{\varepsilon }_{{b_{t}^{i}}}}\) is idiosyncratic noise. The game described by (4) is a random field Ising game, see, e.g., [8].

In what follows we assume that \({{\varepsilon }_{{b_{t}^{i}}}}\) are taken from the Gumbel distribution so that

$${{F}_{ < }}(x) = {\text{Prob}}({{\varepsilon }_{{b_{t}^{i}}}} < x) = {{e}^{{ - {{e}^{{ - \frac{x}{\beta }}}}}}},$$

where the parameter β sets the noise scale.

At each time step \(t \in [0,T]\) we have the following sequence of computations:

1. For each non-vaccinated agent the values of two random variables \({{\varepsilon }_{ + }}\) and \({{\varepsilon }_{ - }}\) are generated independently.

2. For each non-vaccinated agent i we calculate \({{U}_{i}}(b_{t}^{i} = + 1,t)\) and \({{U}_{i}}(b_{t}^{i} = - 1,t)\) using (4) and realisations of the noise variables \({{\varepsilon }_{ + }}\) and \({{\varepsilon }_{ - }}\).

3. Each non-vaccinated agent at time \(t\) decides if he/she supports vaccination (chooses one of \(b_{t}^{i} = \pm 1\)). In particular, the decision is made according to the following expression

$$b_{t}^{i} = \left\{ \begin{gathered} + 1,\quad {\text{if}}\quad {{U}_{i}}( + 1,t) \geqslant {{U}_{i}}( - 1,t), \hfill \\ - 1,\quad {\text{otherwise}}{\text{.}} \hfill \\ \end{gathered} \right.$$
(5)

4. Each non-vaccinated agent \(i\) for whom \(b_{{t - \theta + 1}}^{i} = b_{{t - \theta + 2}}^{i} = \ldots = b_{t}^{i} = + 1\) takes the vaccine; i.e., \(v_{t}^{i} = + 1\).

In the next section we describe the simulations results using the following parameters specifications: \(J = 1\), \({{\tau }^{{\text{I}}}} = 7\), \({{\tau }^{R}} = 90\), \({{p}^{{nv}}} = 0.15\), \({{p}^{v}} = 0.01\), \(T = 100\), \({{t}^{{{\text{start}}}}} = 15\). Other parameters values (\(H,\;{{\chi }_{i}},\;\theta ,\;\beta \)) varied across different scenarios as described in the Results section. Qualitative characteristics of the effects described in the next section are independent from the values of fixed parameters, this values, however, determine quantitative amplitudes and timelines of the considered processes.

RESULTS

In this section we describe the results of numerical simulations of the model. Two assumptions on initial agents’ beliefs \(\{ b_{{t = 0}}^{i}\} \) are studied:

• initial agents’ beliefs are chosen uniformly at random between ±1;

• initial agents’ beliefs are clustered: in two out of four agents communities the common initial belief equals +1, in the other two it equals –1.

We also assume that for all agents \(b_{{t = 0}}^{i} = {{\chi }_{i}}\) for all i so that initial attitudes towards vaccination are assumed to be fully correlated to the attitudes towards generic external information on it.

We first analyse the case of non-clustered uniform distribution of initial beliefs. In Fig. 1 we show the evolution of the fraction of vaccinated population for different levels of noise and two magnitudes of external influence: \(H = 0\) and \(H = 2\). We see that resulting trajectories belong to two distinct classes: (1) “good” ones where evolution results in vaccination of the whole population and (2) “bad” ones where it results in vaccination of only a small fraction of population. Time evolution of mean fractions within these two groups is also presented.

Fig. 1.
figure 1

(Color online) Dynamics of the fraction of vaccinated population for different noise levels and H = (blue lines) 0 and (red lines) 2 for the uniform initial belief distribution. The memory size is \(\theta = 5\). Light solid lines represent trajectories of 100 different simulation algorithm calls. Bright solid and dashed lines show the evolution of the mean fraction over the so called “good” and the “bad” simulations, respectively.

We see that: (1) for stronger noise the number of trajectories reaching equilibrium state with the whole population being vaccinated grows; (2) as noise level increases, the rate at which totally vaccinated equilibrium is reached decreases; (3) the evolution does not depend on the value of H. To make the first conclusion more evident in Fig. 2 we show the dependence of the fraction of “good” simulations on the noise amplitude and memory size.

Fig. 2.
figure 2

(Color online) Dependence between the fraction of simulations with “good” equilibrium and noise amplitude β for different memory sizes θ and uniform initial belief distribution. The fraction is calculated over 100 simulations for each combination of \((\beta ,\theta )\).

The described effect of the growth of the yield of “good” equilibria with growing noise level \(\beta \) results from the asymmetry in zealots formation. In our model only positive zealots may occur. Therefore, when increasing the noise we increase random fluctuations with respect to low noise phase. These fluctuations lead to symmetrical more frequent opinion changes from +1 to –1 and vise versa. However, random fluctuations from –1 to +1 sometimes lead to opinion freezing—vaccination effect, while fluctuations from +1 to –1 do not.

The opinion formation layer influences the SIRS epidemic dynamics through significant decrease in probability to catch the virus. In Fig. 3 we plot the fractions of recovered people in different simulations when the initial agents’ beliefs are uniformly distributed. We see that in the high noise phase the number of simulations with small recovered population fraction is the highest. We plot only the case of \(H = 0\) as for the uniform initial beliefs distribution results for \(H = 0\) and \(H = 2\) are approximately the same.

Fig. 3.
figure 3

(Color online) Fractions of recovered population in 100 simulations in case of the uniform initial belief distribution, \(\theta = 5\), \(H = 0\) and different noise amplitudes β. Mean over simulations infected fractions dynamics is shown in the incut.

We now turn to the case of clustered initial beliefs distribution. In Fig. 4 the vaccinated population fraction for this case is presented. The conclusions about noise amplitude effects remain the same as in the case of uniform initial beliefs distribution. The distinguishing feature of this case is in the effects related to the external influence. We see that the influence of H does strongly depend on the noise amplitude. In the low noise phase an increase of H diminishes the “good” equilibrium yield, while in the case of moderate noise, oppositely, positive field values lead to an increase in this yield. In the high noise phase the results with \(H = 0\) and \(H = 2\) are approximately the same.

Fig. 4.
figure 4

(Color online) Dynamics of fraction of the vaccinated population in case of different noises and different field values in case of the clustered initial beliefs distribution. Memory size \(\theta = 5\). The blue and red lines correspond to H = 0 and 2, respectively. Light solid lines represent trajectories of 100 different simulation algorithm calls. Solid bright lines show the position of the mean fraction over the, so called, “good” simulations and dashed lines—over the “bad” ones.

Figure 5 shows the properties of epidemic spread corresponding to the scenarios described in Fig. 4. The significant effect of the introduction of external influence H for the medium noise case (\(\beta = 5\)) is shown. We see that the fraction of released population drops significantly with increasing \(H\). Considering the results presented in the Fig. 4 for the case \(\beta = 0.1\), one may expect an increase of recovered fraction when we turn on H in the low noise phase. However, this does not happen and the fraction of recovered population is approximately the same as for the case of \(H = 0\) (not presented in the Fig. 5). This happens because differences between \(H = 0\) and \(H = 2\) for \(\beta = 0.1\) in the Fig. 4 arise only around \(t = 40\). At this time point the epidemic is already fading. Note, however, that in case of more complex epidemic spread description the effect will evidently arise.

Fig. 5.
figure 5

(Color online) Fractions of recovered population in 100 simulations in case of the clustered initial beliefs, \(\theta = 5\) and different noise amplitudes β. For noise \(\beta = 5\) both results for \(H = 0\) and \(H = 2\) are shown. For \(\beta = 0.1,\;30\) only results for \(H = 0\) are provided as the results for \(H = 2\) are approximately the same. Mean over simulations infected fractions dynamics is shown in the incut.

CONCLUSIONS

The model of coupled dynamics of vaccination opinion formation and SIRS epidemic in the two-layer network is presented. The model describes non-trivial dependence between the properties of vaccination dynamics and (1) noise amplitude, (2) initial opinion clustering in the society and (3) magnitude of external influence. The simulation results show in particular that if initial social opinion about vaccination is clustered and the level of noise is low, the presence of generic external influence may worsen the situation. The proposed model might in future be developed in several directions. In particular, it is promising to investigate the influence of the “hubs” opinion and to introduce the feedback of epidemic spread properties onto vaccination opinion formation process.