1 Introduction

While typically non-life-threatening for healthy individuals, seasonal influenza is responsible for tens of thousands of deaths (Centers for Disease Control and Prevention 2010a) and tens of billions of dollars of lost earnings (Molinari et al. 2007) each year in the United States alone. Given the capacity of universal vaccination to mitigate these consequences and to protect especially vulnerable populations (i.e. children, pregnant women and the elderly), the United States Centers for Disease Control and Prevention recommends that everyone 6 months of age or older get a flu vaccine every year (Centers for Disease Control and Prevention 2010b). Despite the CDC recommendation, each year, roughly half of the American population still does not get a flu shot (Centers for Disease Control and Prevention 2016).

The paper Galvani et al. (2007) reports that this empirically observed pattern is in accord with predictions of a model that conceptualizes vaccination decisions as strategies in a multi-player game. In these so-called vaccination games, rational individuals are assumed to reach a Nash equilibrium at which each of them follows a strategy that minimizes each individual’s expected cost given what the other individuals in the population do. At Nash equilibrium, the vaccination coverage will not be sufficient for herd immunity, and the overall cost to the population will not be minimized. This misalignment between what is optimal for a society (vaccination levels at the herd immunity threshold) and what is optimal for individuals (freeloading on the rest of the population’s vaccination) is often referred to as the vaccination dilemma (Fu et al. 2011a). The vaccination dilemma was first described in Fine and Clarkson (1986), and then later independently in Geoffard and Philipson (1997), albeit not phrased in game-theoretic terminology. The first papers that cast this observation in terms of Nash equilibria in vaccination games were Bauch et al. (2003) and Bauch and Earn (2004). Since these seminal papers appeared, the study of vaccination games has mushroomed. A recent survey article Wang et al. (2016) includes a total of 777 citations, including one book-length treatment Manfredi and d’Onofrio (2013).

While these general results on vaccination games explain how individual vaccination decisions are likely to lead to lower vaccination coverages than the societally optimal herd immunity threshold, they are based on several assumptions that will not be satisfied in a real population. People may not have an accurate perception of the costs of vaccination (such as likelihood of side effects) and infection, and they don’t always behave rationally. On the one hand, these factors can exacerbate the dilemma inherent in voluntary vaccination. On the other hand, this creates opportunities for designing public policy that would eliminate or alleviate the vaccination dilemma by offering appropriate incentives or effective dissemination of useful information. Therefore most of the current research on vaccination games focuses on understanding the processes by which people arrive at their decisions to vaccinate or remain unvaccinated, and how these processes will influence the resulting vaccination coverage.

In this paper we focus on the role of imitation. Imitation is prevalent in much of everyday decision-making, in particular when the environment is complex or largely unknown. It can be a very successful procedure for finding advantageous strategies in social games (Rendell et al. 2010). Social scientists and psychologists have long recognized the importance of imitation, and it has recently moved into the focus of economists (Apesteguia et al. 2007). An important concern in the study of imitation is that inertia, resistance to change, may be present in nearly all decision-making processes (Forsell and Aström 2012).

The influential paper Fu et al. (2011a) presented a model that incorporates imitation into the decision-making process for vaccinations against flu-like infections. In this model it is assumed that during a vaccination campaign that precedes the actual outbreak of the disease, each focal individual i independently makes a decision based on comparing his or her own costs C(i) in the preceding season with the cost C(j) of one other randomly chosen individual j. Then i either follows the same strategy as in the previous season, or switches to j’s strategy of the previous season. The probability of switching is given by a so-called Fermi function

$$\begin{aligned} p_{switch}(i\rightarrow j) = \frac{1}{1 + e^{-\beta (C(i) - C(j))}}. \end{aligned}$$
(1)

Thus the strategy of a much better performing player is readily adopted, whereas it is unlikely, but not impossible, to adopt the strategies of worse performing players. The parameter \(\beta \) incorporates the uncertainties in the strategy adoption, originating in either the variation of costs or in mistakes in the decision-making process. In the limit \(\beta \rightarrow 0^+\) player i is unable to retrieve any information from player j and switches to the strategy of j by tossing a fair coin (Hauert and Szabó 2005).

Fermi functions take their origin from the Fermi-Dirac distribution of statistical physics; see, for example, Reif (2008) or Schroeder (1999). The first use of Fermi functions for determining the probabilities of switching to another strategy is usually attributed in the literature to Blume (1993). Fermi functions are examples of what are referred to as “smoothed imitation” (Szabó and Fáth 2007). Other examples of smoothed imitation would be functions of the form

$$\begin{aligned} p_{switch}(i\rightarrow j) = \mu + \frac{\nu }{\alpha + e^{- \beta (C(i) - C(j)) + \gamma }} \end{aligned}$$
(2)

for some suitable choice of parameters \(\alpha , \mu \), and \(\nu \) that gives probabilities.

There is evidence that switching probabilities as in (2) may be more realistic. In fact, Traulsen et al. (2010) reports the results of behavioral experiments on imitation in prisoner’s dilemma games. Analysis of observed behaviors and curve fitting lead to probabilities of switching to the other strategy of the form (2) with \(\alpha = 1, \nu = 1- \mu \) and \(\mu = 0.28 \pm 0.07, \beta = 0.67 \pm 0.28, \gamma = -\,0.11 \pm 0.23\) for the switch from cooperation to defection. For the switch from defection to cooperation they found \(\mu = 0.25 \pm 0.01, \beta = 0.99 \pm 0.23, \gamma = 0.79 \pm 0.14\).

While in the curve fitting performed by Traulsen et al. (2010) the parameter \(\alpha \) was always fixed at \(\alpha =1\), as a matter of mathematical convenience, we will here investigate models where \(\alpha \) can vary, but \(\mu , \nu \), and \(\gamma \) are fixed. Lemma 3 below then shows that each of our models has a counterpart in the general form of (2) with \(\alpha = 1\) that exhibits the same qualitative features of the dynamics as our models, in particular, admits equilibria that are arbitrarily close to the societal optimum. See also the discussion surrounding (2) below.

While in ODE-based models such as Bauch (2005) the likelihood of imitation is often assumed to be proportional to the difference in costs, Fermi functions are used in almost all published discrete-time models of vaccination games with imitation. Sometimes the cost of the focal player is compared with an average cost of several other players (Fukuda et al. 2014; Ichinose and Kurisaku 2017; Iwamura and Tanimoto 2018; Li et al. 2017), but the probability of switching still follows the pattern of (1). In view of the results of Traulsen et al. (2010), the question naturally arises whether the particular form of smoothing functions given by (1) significantly influences the predictions of models based on it; a question that has received surprisingly little attention in the literature thus far. A notable exception is Zhang et al. (2011), where two distinct parameter settings in (2) that depend on the current strategy of the focal player were assumed. For this type of setup, it is intuitively clear that the equilibrium may shift towards the strategy with the more favorable parameters for being imitated.

In this work, we will demonstrate that the choice of the form of the smoothing function itself can significantly alter the model’s predictions even if the parameters in (2) do not depend on the current strategy of the focal player. Most notably, we show that a suitable choice of the smoothing function alone can drive the system to equilibrium vaccination coverages that are arbitrarily close to the societal optimum of herd immunity.

There are four parameters in (2) in addition to the \(\beta \) of (1), but for mathematical convenience, we will focus on models with switching probabilities of the form

$$\begin{aligned} p_{switch}(i\rightarrow j) = \frac{1}{\alpha + e^{-\beta (C(i) - C(j))}} \end{aligned}$$
(3)

that has only one additional parameter \(\alpha \ge 1\). As will be shown in Section 4, this modified switching probability models a situation where individuals only rarely compare their costs with others, but are eager to switch when they do. Moreover, for every model that uses (2) with \(\gamma < 0\) there exists a corresponding model with switching probabilities of the form (3) that for the same cost and disease transmission parameters predicts the same equilibria and intervals where the vaccination coverage decreases or increases. The additional parameters of (2) may influence the stability of the equilibria and how fast vaccination coverages approach them, but they will not influence the positions of the equilibria; see Lemma 3. For a more detailed discussion of generalized Fermi functions, in particular of the role of the parameter \(\alpha \), see Sect. 4.

We focus here on the case where the cost of vaccination is low relative to the cost of the disease, which is the realistic one for flu vaccinations. We also assume uniform mixing of the population to eliminate all aspects of the structure of contact networks that may confound the effects of using the modified switching probability (3) in place of (1). Under these assumptions, for \(\alpha = 1\) it is reported in Fu et al. (2011a) that the equilibrium coverage is always even lower than the Nash equilibrium, which is already lower than the societal optimum at herd immunity. In stark contrast, we found that for all sufficiently large values of \(\beta \), as long as we also choose \(\alpha \) large enough relative to \(\beta \), the equilibrium vaccination coverage predicted by our model can be arbitrarily close to the societally optimal value of herd immunity.

2 Our model

Our model assumes an infection with seasonal outbreaks that follow dynamics of type SIR in each season with no carry-over of immunity to the next season. These assumptions would, in particular, apply to influenza in regions or countries where outbreaks are seasonal. It is a deterministic difference equation model that predicts the time evolution of vaccination coverage from season to season. We designed the model in such a way that it is as similar as possible to the version of the model of Fu et al. (2011a) with uniform mixing. Thus, in particular, the vaccine is assumed to be 100% effective, each individual is assumed to base the vaccination decision only on comparison of previous costs with one other individual, the stochasticity inherent in individual vaccination decisions enter the model only via mean values, costs of vaccination and infection are conceptualized as mean costs, and in addition to the network-based versions the models of Fu et al. (2011a) we ignore the patterns by which people make contact. This exact replication of the assumptions of  the non-network model in Fu et al. (2011a) allows us to cleanly disentangle the effects of our new parameter \(\alpha \) from all other potentially confounding aspects.

Each time step n represents a year or flu season, and \(V_n\) represents the proportion of individuals in the population who decide to get vaccinated in season n. Decisions on whether or not to get vaccinated may depend on individual experience in the previous flu season, the current strategy of a host, and on imitation of one randomly chosen other host. These decisions are assumed to be made by all hosts independently and simultaneously prior to any flu outbreak. They collectively determine the vaccination coverage \(V_n\) in season number n. After individuals make their vaccination decisions, the probability \(x_n = x(V_n)\) of infection of any unvaccinated individual in flu season number n is then calculated based on a standard SIR model. Our approach of combining a discrete vaccination dynamical system with a continuous model of the inter-year disease dynamics is similar to that of Vardavas et al. (2007), which was one of the first game-theoretic models of influenza vaccination.

The above description is written in the language of an agent-based stochastic process, but we did not actually implement and study the model for finite populations. Instead, our version assumes a very large population and is a deterministic compartment-level model, based on expected proportions.

We let \(c_v > 0\) denote the mean cost of vaccination, and \(c_i > 0\) denote the mean cost of infection. Costs are treated as fixed positive numbers here that represent average costs. Throughout this paper, we assume that the relative cost of vaccination defined as \(c = \frac{c_v}{c_i}\) satisfies the inequality \(c < 0.5\). This seems realistic for infections like seasonal influenza (Fu et al. 2011b).

The model is initialized by randomly assigning strategies for the first season. Before the next flu season \(n+1\), each player i updates his or her strategy as follows:

  • Player i picks a randomly chosen other player j.

  • Player i compares his or her own actual cost C(i) in the current season to the actual cost C(j) of player j in the current season.

  • Player i switches to player j’s strategy with probability

    $$\begin{aligned} p_{switch}(i\rightarrow j) = \frac{1}{\alpha + e^{-\beta (C(i) - C(j))}} \end{aligned}$$
    (4)

    and retains the current strategy with probability \(1 - p_{switch}(i\rightarrow j)\).

In the second stage of each season, after all individuals in the population have made their vaccination decisions, there will be a flu outbreak. It is assumed to develop according to a standard ODE-based SIR model with basic reproductive ratio \(\mathcal {R}_0> 1\) that is kept fixed over all seasons. For seasonal influenza in the continental U.S., typical values for \(R_0\) range from 1.5 to 3 (Fu et al. 2011b). Our simulations reported here were done for \(R_0 = 2.5\) that exemplifies this range, and for \(R_0 = 1.3\) and 5 to test robustness of our findings.

More precisely, the outbreak in a given season is assumed to occur according to a model of the formFootnote 1

$$\begin{aligned} \begin{aligned} s'&= -\xi si, \\ i'&= \xi si - \omega i, \\ r'&= \omega i, \end{aligned} \end{aligned}$$
(5)

where s, i and r denote the proportions of the population that are susceptible, infectious and recovered, respectively. The basic reproductive ratio of the model is \(\mathcal {R}_0=\xi /\omega > 1\). The limit \(s_\infty \) of the fraction s of susceptible individuals as \(t \rightarrow \infty \) depends only on \(\mathcal {R}_0\) and satisfies [see, for example Diekmann et al. (2013)]:

$$\begin{aligned} \ln s_\infty - \ln s_0 = \mathcal {R}_0\left( s_\infty - s_0\right) . \end{aligned}$$
(6)

Note that \(t = \infty \) here should be interpreted as the end of the given season.

In flu season n, the initial proportion of susceptible individuals is \(s_0 = 1 - V_n\) and the probability of each susceptible individual being infected,

$$\begin{aligned} x(V_n) = \frac{s_0 - s_\infty }{s_0}, \end{aligned}$$

is then calculated by solving (6). Note that this derivation is valid under the assumption of a 100% effective vaccine where the probability of being infected is 0 for a vaccinated individual.

Recall that we model the dynamics of vaccination coverage from season to season, where \(V_n\) is the expected proportion of hosts who decide to vaccinate in season n. Now let us derive a formula for \(V_{n+1}\) as a function of \(V_n\) from the underlying assumptions of our model about individual decision-making. Consider a very large population of size N. Note that in the first stage of season \(n+1\), it is expected that a proportion of \(1-V_n\) among approximately \(V_n N\) focal individuals who did vaccinate in season n will compare their cost \(c_v\) with the cost born by individuals who did not vaccinate. For individuals who did not vaccinate, the probability of infection in season n was \(x_n = x(V_n)\). Thus approximately \(x_nV_n(1-V_n)N\) of those comparisons will be with individuals that experienced infection with a cost of \(c_i\), and approximately \((1-x_n)V_n(1-V_n)N\) will have escaped infection and born no cost whatsoever. As a result of these comparisons, some players decide to switch to the other strategy, in accordance with the switching probabilities given by (4). Therefore, the number of individuals that got vaccinated in season n but do not in season \(n+1\) is given by

$$\begin{aligned} V_n N (1-V_n) \left( x_n \frac{1}{\alpha + e^{-\beta (c_v-c_i)}} + (1-x_n) \frac{1}{\alpha + e^{-\beta (c_v-0)}} \right) . \end{aligned}$$

An analogous calculation shows that the number of individuals that did not vaccinate in season n but do in season \(n+1\) is given by

$$\begin{aligned} (1-V_n)NV_n \left( x_n \frac{1}{\alpha + e^{-\beta (c_i-c_v)}} + (1-x_n) \frac{1}{\alpha + e^{-\beta (0-c_v)}} \right) . \end{aligned}$$

Dividing by N, we obtain the expected change \(\varDelta (n)\) in the proportion of vaccinators

$$\begin{aligned} \begin{aligned} \varDelta (n)&= V_{n+1} - V_n \\&= (1-V_n)V_n \left( \frac{x_n}{\alpha +e^{-\beta (c_i-c_v)}} + \frac{1-x_n}{\alpha +e^{\beta c_v}} - \frac{x_n}{\alpha +e^{-\beta (c_v-c_i)}} - \frac{1-x_n}{\alpha +e^{-\beta c_v}}\right) , \end{aligned} \end{aligned}$$
(7)

Thus our difference equation model for the evolution of flu vaccination is given by the updating function \(J: [0,1] \rightarrow [0,1]\) that maps \(V_n\) to \(V_{n+1}\) and is defined as follows:

$$\begin{aligned} \begin{aligned} J(V_n)&= V_{n+1} = V_n + \varDelta (n)\\&= V_n + (1-V_n)V_n \\&\quad \times \left( \frac{x_n}{\alpha +e^{-\beta (c_i-c_v)}} + \frac{1-x_n}{\alpha +e^{\beta c_v}} - \frac{x_n}{\alpha +e^{-\beta (c_v-c_i)}} - \frac{1-x_n}{\alpha +e^{-\beta c_v}}\right) . \end{aligned} \end{aligned}$$
(8)

Here \(x_n\) is the fraction calculated according to the epidemic model (5), so that

$$\begin{aligned} x_n = x(V_n) = \frac{s_0 - s_\infty }{s_0} = \frac{ (1-V_n) - s_\infty }{1-V_n} \end{aligned}$$

and \(s_\infty \) is the solution to

$$\begin{aligned} \ln s_\infty - \ln (1-V_n) = \mathcal {R}_0\left( s_\infty - (1-V_n)\right) . \end{aligned}$$

3 Predictions of the model

3.1 Preliminary observations: costs, Nash equilibria, and the societal optimum

The expected costs \(c_n^V\) and \(c_n^U\) for vaccinators and nonvaccinators in season n will be:

$$\begin{aligned} c_n^V = c_v \quad \text{ and } \quad c_n^U = c_i x_n. \end{aligned}$$

Let

$$\begin{aligned} V_{hit} = 1 - \frac{1}{R_0} \end{aligned}$$

denote the vaccination coverage at the herd immunity threshold. Then the probability x of an unvaccinated individual getting infected is a strictly decreasing function on the interval \([0, V_{hit}]\) such that \(x(V_{hit}) = 0\). Thus for vaccination coverage \(V_n = V_{hit}\) we would have \(c_n^V = c_v > c_n^U = 0\) and the rational choice would be not to vaccinate. Conversely, when \(c_i x(0) > c_v \) and \(V_n = 0\), then \(c_n^V < c_n^U\), and the rational choice would be to vaccinate, which implies the existence of a unique Nash equilibrium with vaccination coverage \(V_{Nash} \in (0,V_{hit})\).

Now let us consider a societally optimal vaccination coverage \(V_{opt}\). This is supposed to minimize the following function that represents the average cost to the entire population:

$$\begin{aligned} PC(V) = c_v V + c_i (1-V) x(V), \end{aligned}$$

where for \(V \in [0, V_{hit}]\) the function \(x = x(V)\) is determined by the SIR model of the outbreak.

An important observation is that when \(c_i \ge c_v\), then PC(V) is strictly decreasing on the interval \([0, V_{hit}]\) so that, in particular, \(V_{hit}\) is the societally optimal vaccination coverage. Proofs of this observation can be found, for example, in Fu et al. (2011b) and Xin et al. (2018).

3.2 Preliminary observations: equilibria in our model

Let \(J: [0,1] \rightarrow [0,1]\) denote the updating function that maps the vaccination coverage \(V_n\) in season n to the vaccination coverage \(V_{n+1}\) in the next season. We were primarily interested in finding equilibria \(V^*\) where \(J(V^*) = V^*\) and determining their stability properties. Note that the meaning of “equilibrium” is different here than in the phrase “Nash equilibrium.” In our model players are assumed to base their decisions on imitation, while Nash equilibria result from a different model that assumes rational choice.

Notice that the updating function J that maps \(V_n\) to \(V_{n+1}\) and is defined in (8) maps [0, 1] continuously into [0, 1]. If \(V_n = 1\), then all individuals vaccinate, and if \(V_n = 0\), then nobody vaccinates. In either case, no player can switch strategies, since there is nobody in the population to imitate who would follow another strategy. So \(J(0) = 0\), \(J(1) = 1\), and \(V^{**} := 0\), \(V^{***} := 1\) are always equilibria.

However, we would be more interested in interior equilibria \(V^* \in (0,1)\) for our model. Note that in our conceptualization of imitation there is always a positive probability that a player will stick with the current strategy rather than imitate another one, even if that other strategy gave a vastly lower cost. Thus \(J(V) = 0\) only if \(V=0\) and \(J(V) = 1\) only if \(V=1\), which means that the interior (0, 1) is invariant under J. If both \(V^{**}\) and \(V^{***}\) are repelling, then the updating function J maps some interval \([\varepsilon , 1-\varepsilon ]\) into itself, and at least one such interior equilibrium \(V^*\) is guaranteed to exist by Brouwer’s fixed point theorem.

If \(V^*\) exists, then it must be in the interval \((0, V_{hit})\). To see this, consider \(V_n \in [V_{hit}, 1)\). For very large population sizes N, we will observe approximately \(V_n(1-V_n)N\) comparisons that may induce a player to switch from vaccinating to not vaccinating, and the same number of comparisons that may induce a player to switch from not vaccinating to vaccinating. In all these comparisons, players who vaccinated will have born a positive cost, while players who did not vaccinate will not have born any cost. Thus for any choice of the parameters \(\alpha , \beta \) we will observe more switches from vaccinating to not vaccinating than vice versa, and it follows that \(V_{n+1} < V_n\).

In particular, it follows that the equilibrium \(V^{***} = 1\) is always repelling. Lemma 2 of the next subsection shows that an interior equilibrium exists in our model only if \(V^{**} = 0\) is repelling and gives precise conditions on the parameters for when this is the case. Lemma 1 of the next subsection shows that if an interior equilibrium \(V^*\) exists, it must be unique.

3.3 Our main theorem

The following theorem is the main result of this paper. It shows that for all sufficiently large \(\beta \) and all choices of \(\alpha \) that are sufficiently large relative to \(\beta \) the interior equilibrium \(V^*\) will be arbitrarily close to the societal optimum \(V_{hit}\).

Theorem 1

Fix any \(0 \le V^- < V_{hit}\) and assume the parameters \(\alpha , \beta \) of our model satisfy the inequalities

$$\begin{aligned}&1-e^{-2\beta (c_i-c_v)} - \frac{2(1-x(V^-))}{x(V^-)}e^{-\beta (c_i-2c_v)}> 0. \\&\quad \alpha > \text{ max }\{1,\ e^{\beta (c_i-c_v)} + e^{-\beta (c_i-c_v)}-2e^{\beta c_v}-2e^{-\beta c_v}\}. \end{aligned}$$

Then

  1. (a)

    If \(0< V_n < V^-\), then \(V_{n+1} > V_n\).

  2. (b)

    If an interior equilibrium \(V^*\) exists, then \(V^* \ge V^-\).

The inequalities in Theorem 1 only provide sufficient conditions that serve the purpose of giving a rigorously derived qualitative result rather than providing good estimates of the thresholds. In order to keep the complexity of the calculations within reasonable limits, we used considerable simplifications, so the inequalities are far from sharp. For example, with \(R_0 = 2.5\), \(c_i = 12\), \(c_v = 1\) and \(V^- = 0.5\), the first inequality implies that \(\beta > \beta _0\), where \(\beta _0\) is approximately 0.1281. In the second inequality, the \(\alpha \) threshold value increases exponentially fast as \(\beta \) increases. Under the same parameter setting, when \(\beta = 0.1281\), the second inequality gives \(\alpha > \alpha _0\) where \(\alpha _0\) is 1, when \(\beta = 0.2\), \(\alpha _0\) is approximately 5.0555, and when \(\beta = 1\), \(\alpha _0\) is approximately 59,868.

Proof of Theorem 1

Note that it suffices to prove part (a); part (b) is then an immediate consequence.

Let \(V^-, \alpha , \beta \) be as in the assumption, and let \(V_n < V^-\).

We will show that the sign of \(\varDelta (n) = V_{n+1}-V_n\) is positive. From (7) we get:

$$\begin{aligned} \begin{aligned} \varDelta (n)&= V_{n+1} - V_n \\&= (1-V_n)V_n\left( \frac{x_n}{\alpha +e^{-\beta (c_i-c_v)}} + \frac{1-x_n}{\alpha +e^{\beta c_v}} - \frac{x_n}{\alpha +e^{-\beta (c_v-c_i)}} - \frac{1-x_n}{\alpha +e^{-\beta c_v}}\right) \\&= (1-V_n)V_n\left( \frac{x_n(\alpha +e^{-\beta (c_v-c_i)}) - x_n(\alpha +e^{-\beta (c_i-c_v)})}{(\alpha +e^{-\beta (c_i-c_v)})(\alpha +e^{-\beta (c_v-c_i)})}\right) \\&\quad + (1-V_n)V_n\left( \frac{(1-x_n)(\alpha +e^{-\beta c_v})-(1-x_n)(\alpha +e^{\beta c_v})}{(\alpha +e^{\beta c_v})(\alpha + e^{-\beta c_v})}\right) \\&= (1-V_n)V_n\left( \frac{x_n(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)})}{(\alpha +e^{-\beta (c_i-c_v)})(\alpha +e^{\beta (c_i-c_v)})} + \frac{(1-x_n)(e^{-\beta c_v}-e^{\beta c_v})}{(\alpha + e^{\beta c_v})(\alpha +e^{-\beta c_v})}\right) \\&= \frac{(1-V_n)V_n}{(\alpha +e^{\beta c_v})(\alpha + e^{-\beta c_v})}\left( quot(\alpha ,\beta )x_n(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) \right. \\&\quad \left. +\, (1-x_n)(e^{-\beta c_v} - e^{\beta c_v})\right) , \end{aligned} \end{aligned}$$

where

$$\begin{aligned} quot(\alpha , \beta ) = \frac{(\alpha +e^{\beta c_v})(\alpha +e^{-\beta c_v})}{(\alpha +e^{\beta (c_i-c_v)})(\alpha + e^{-\beta (c_i-c_v)})}. \end{aligned}$$

Let

$$\begin{aligned} f(\alpha , \beta ) = quot(\alpha ,\beta )x_n(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) + (1-x_n)(e^{-\beta c_v} - e^{\beta c_v}). \end{aligned}$$

Now to show that \(\varDelta (n) > 0\), it suffices to show that \(f(\alpha , \beta ) > 0\).

For \(\beta , V_n\) as in the assumptions we have \(x_n > x(V^-)\), and hence

$$\begin{aligned} \begin{aligned} 1-e^{-2\beta (c_i-c_v)} - \frac{2(1-x(V^-))}{x(V^-)}e^{-\beta (c_i-2c_v)}&> 0,\\ 1-e^{-2\beta (c_i-c_v)} - \frac{2(1-x_n)}{x_n}e^{-\beta (c_i-2c_v)}&> 0,\\ 1-e^{-2\beta (c_i-c_v)} + \frac{2(1-x_n)}{x_n}e^{-\beta c_i} - \frac{2(1-x_n)}{x_n}e^{-\beta (c_i-2c_v)}&> 0. \end{aligned} \end{aligned}$$

Then

$$\begin{aligned} \begin{aligned}&\frac{1}{2}x_n(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) + (1-x_n)(e^{-\beta c_v}- e^{\beta c_v})\\&\quad = \frac{x_n}{2}e^{\beta (c_i-c_v)}\left( 1-e^{-2\beta (c_i-c_v)} + \frac{2(1-x_n)}{x_n}e^{-\beta c_i} - \frac{2(1-x_n)}{x_n}e^{-\beta (c_i-2c_v)}\right) > 0. \end{aligned} \end{aligned}$$

Moreover, the following inequalities are all equivalent:

$$\begin{aligned} \begin{aligned}&quot(\alpha , \beta ) = \frac{(\alpha +e^{\beta c_v})(\alpha +e^{-\beta c_v})}{(\alpha +e^{\beta (c_i-c_v)})(\alpha + e^{-\beta (c_i-c_v)})}> \frac{1}{2},\\&2\alpha ^2 + (2e^{\beta c_v} + 2e^{-\beta c_v})\alpha + 2> \alpha ^2 + (e^{\beta (c_i-c_v)}+e^{-\beta (c_i-c_v)})\alpha + 1,\\&\alpha ^2 + (2e^{\beta c_v} + 2e^{-\beta c_v}- e^{\beta (c_i-c_v)} - e^{-\beta (c_i-c_v)} )\alpha + 1 > 0. \end{aligned} \end{aligned}$$

Thus for \(\alpha > \text{ max }\{1,\ e^{\beta (c_i-c_v)} + e^{-\beta (c_i-c_v)}-2e^{\beta c_v}-2e^{-\beta c_v}\}\) we have

$$\begin{aligned} quot(\alpha , \beta ) > \frac{1}{2}. \end{aligned}$$

Then

$$\begin{aligned} \begin{aligned} f(\alpha , \beta )&= quot(\alpha ,\beta )x_n(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) + (1-x_n)(e^{-\beta c_v} - e^{\beta c_v})\\&> \frac{1}{2}x_n(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) + (1-x_n)(e^{-\beta c_v} - e^{\beta c_v})\\&> 0. \end{aligned} \end{aligned}$$

\(\square \)

Lemma 1

There can be at most one equilibrium \(V^*\) in the interval (0, 1).

Proof

In this argument we will treat x as a variable that is a function of V. This function is strictly decreasing, and hence invertible, on the interval \([0, V_{hit}]\). Consider the functions:

$$\begin{aligned} \begin{aligned} g(x)&:= quot(\alpha ,\beta )(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)})x + (1-x)(e^{-\beta c_v} - e^{\beta c_v}),\\ G(V)&:= quot(\alpha ,\beta )(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)})x(V) + (1-x(V))(e^{-\beta c_v} - e^{\beta c_v}), \end{aligned} \end{aligned}$$
(9)

where \(quot(\alpha , \beta )\) is as in the proof of Theorem 1. It follows from our calculations of \(\varDelta (n)\) in that proof that at an interior equilibrium \(V^*\) with \(x^* := x(V^*)\) we must have \(g(x^*) = G(V^*) = 0\).

Note that g(x) is a linear function with slope

$$\begin{aligned} m = quot(\alpha ,\beta )(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) + e^{\beta c_v} - e^{-\beta c_v} > 0. \end{aligned}$$
(10)

Thus the system can have at most one interior equilibrium. \(\square \)

Lemma 2

The following conditions are equivalent:

  1. (a)

    An interior equilibrium \(V^* \in (0, V_{hit})\) exists.

  2. (b)

    The equilibrium \(V^{**} = 0\) is repelling.

  3. (c)

    One of the following equivalent inequalities holds:

    $$\begin{aligned} \begin{aligned}&quot(\alpha ,\beta )x(0)(e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}) + (1-x(0))(e^{-\beta c_v} - e^{\beta c_v})> 0,\\&x(0)\left( quot(\alpha ,\beta )\left( e^{\beta (c_i-c_v)} -e^{-\beta (c_i-c_v)} \right) + e^{\beta c_v} - e^{-\beta c_v} \right)> e^{\beta c_v} - e^{-\beta c_v},\\&x(0) > \frac{1}{quot(\alpha ,\beta )\left( \frac{e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}}{e^{\beta c_v}-e^{-\beta c_v}}\right) +1}. \end{aligned} \end{aligned}$$

Proof

Let g(x), G(V) be defined as in (9). Then the conditions in part (c) are simply saying that \(g(x(0)) = G(0) > 0\).

Here \(x(0) = x(V^{**}) \) is the predicted final size of an outbreak with no vaccination whatsoever. Note that the function g(x) is linear in x and increasing by (10), while x(V) is nonincreasing, and strictly decreasing on \([0, V_{hit}]\). Thus G(V) will be strictly decreasing on the interval \([0, V_{hit}]\). Since \(x(V_{hit}) = 0\) and thus

$$\begin{aligned} G(V_{hit}) = e^{-\beta c_v} - e^{\beta c_v} < 0, \end{aligned}$$

it follows from the IVT that \(V^*\) exists if, and only if, \(G(0) > 0\).

Here and in the next subsection we will work with the following formula for the updating function J as defined in (8):

$$\begin{aligned} \begin{aligned} H(V)&:= \frac{(1-V)V}{(\alpha +e^{\beta c_v})(\alpha + e^{-\beta c_v})},\\ J(V)&= V + H(V)G(V). \end{aligned} \end{aligned}$$
(11)

Thus when \(G(0) > 0\), then \(J(\varepsilon ) > \varepsilon \) for sufficiently small \(\varepsilon > 0\) and \(V^{**} = 0\) will be repelling. On the other hand, when \(G(0) \le 0\), then we must have \(G(V) < 0\) for all \(V \in (0,1]\), the equilibrium \(V^{**}\) will be locally asymptotically stable and globally attracting on [0, 1), while \(V^*\) does not exist. \(\square \)

Lemma 2 predicts that the vaccination coverage in a population that starts with \(V(0) \in (0,1)\) will evolve towards the equilibrium \(V^{**} = 0\) where nobody vaccinates if, and only if,

$$\begin{aligned} x(0) \le \frac{1}{quot(\alpha ,\beta )\left( \frac{e^{\beta (c_i-c_v)}-e^{-\beta (c_i-c_v)}}{e^{\beta c_v}-e^{-\beta c_v}}\right) +1}. \end{aligned}$$
(12)

When the inequality in (12) is reversed, the vaccination coverage may evolve towards an interior equilibrium \(V^*\) even when \(V^{**} = 0\) is the Nash equilibrium. For details see Xin et al. (2018).

3.4 Simulated evolution of vaccination coverage

Figure 1 shows equilibrium vaccination coverage, \(V^*\), as a function of the parameters determining the shape of the switching probabilities, \(\alpha \) and \(\beta \). Vaccine coverage at equilibrium is calculated via numerical simulation of the model in Eq. (8) with stopping criteria of \(|V_{n+1} -V_n| < 10^{-5}\). In our simulations, the cost parameter values are chosen according to those estimated in Galvani et al. (2007), which analyzed responses to surveys of health opinions, behaviors and outcomes from the Health Promotion at Work longitudinal study of university employees. These parameter values were also used in Fu et al. (2011a). For a detailed derivation of the parameter values from the findings of Galvani et al. (2007), see Fu et al. (2011b).

The locations of apparent interior equilibria \(V^*\) that were found in these simulations closely match those predicted by the theoretical analysis; see Figure 10 of Xin et al. (2018). To demonstrate the result of Theorem 1, contours are included for the vaccination coverage achieved at the Nash equilibrium.

The simulation results in Fig. 1 qualitatively agree with the behavior predicted by Theorem 1. Specifically, we see that open-minded imitation leads to vaccination coverage higher than that of the Nash equilibrium provided that \(\beta \) is not excessively small and that \(\alpha \) is sufficiently large relative to \(\beta \). Moreover, numerical simulation in the finite window of \(\alpha , \beta \) produces maximum vaccination coverage levels that are quite close to the societal optima (e.g. maximum vaccine coverages of \(V=0.597, 0.798\) compared with \(V_{hit} = 0.6, 0.8\) in Fig. 1b, c, respectively). While it is true that the vaccination coverages obtained with open-minded imitation do not outperform those of the Nash equilibrium by large amounts, we note that this is a consequence of the fact that Nash equilibrium coverages are close to societal optimal because the costs of vaccination and infection differ significantly [i.e. \(c_v=1\) and \(c_i=12\) as in Fu et al. (2011a, b)].

Fig. 1
figure 1

Equilibrium vaccine coverage, \(V^*\), as a function of \(\alpha \) and \(\beta \) via numerical simulation of Eq. (8) with \(c_v = 1\), \(c_i = 12\) and stopping criteria of \(|V_{n+1} -V_n| < 10^{-5}\). Labeled contours illustrate where imitation matches the vaccination coverage of the Nash equilibrium, \(V_{Nash}\), and the societal optimal (i.e. herd immunity threshold), \(V_{hit}=1-1/\mathcal {R}_0\), when applicable. All simulations start at initial vaccination coverage of \(V_0=0.5\). Typical estimates of \(R_0\) for seasonal influenza range from around 1.5 to 3 (Fu et al. 2011b)

Full size image

The lower right corner of Fig. 1a (i.e. simulation results for large \(\alpha \) and small \(\beta \)) shows the situation where imitation appears to produce vaccination coverage that is even higher than the societal optimum. However, recall from Section 1 that the focal player is unable to retrieve any information from comparing to another player as \(\beta \rightarrow 0^+\) and note that \(p_{switch}\) is a decreasing function of \(\alpha \). Therefore, the high levels of vaccination coverage are not a result of the imitation dynamics but rather a consequence of the fact that the values of \(\alpha \) and \(\beta \) make the probability of switching strategies so low that vaccine coverage stays at or near the initial vaccination coverage of \(V_0 = 0.5\) for a long time.

3.5 Stability of the interior equilibrium

When \(\alpha \) is sufficiently large, the interior equilibrium will be locally stable and will be monotonically approached. However, for \(\alpha \) close to 1 trajectories may approach \(V^*\) via damped oscillations, or \(V^*\) may even become unstable, with trajectories exhibiting sustained oscillations instead of approaching \(V^*\).

Let us assume that the interior equilibrium \(V^*\) exists and work with the function J defined in (11). A sufficient condition for stability of \(V^*\) is given by

$$\begin{aligned} -1< \frac{dJ}{dV}(V^*) < 1. \end{aligned}$$

Similarly, a sufficient condition for monotone approach to \(V^*\) is given by

$$\begin{aligned} 0< \frac{dJ}{dV}(V^*) < 1. \end{aligned}$$

By differentiating J with respect to V we find that for all \(V \in (0,1)\):

$$\begin{aligned} \frac{dJ}{dV}(V) = 1 + \frac{dH}{dV}(V)G(V) + H(V)\frac{dG}{dV}(V). \end{aligned}$$

At the interior equilibrium \(V^*\) we have \(G(V^*) = 0\), so that:

$$\begin{aligned} \begin{aligned} \frac{dJ}{dV}(V^*)&= 1 + \frac{dH}{dV}(V^*)G(V^*) + H(V^*)\frac{dG}{dV}(V^*) \\&= 1 + \frac{dH}{dV}(V^*)(0) + H(V^*)\frac{dG}{dV}(V^*) = 1 + H(V^*)\frac{dG}{dV}(V^*)\\&= 1 + H(V^*) m \frac{dx}{dV}(V^*) \\&< 1. \end{aligned} \end{aligned}$$
(13)

Since \(H(V) > 0\) for \(V \in (0,1)\), the last inequality in (13) follows from (10) and the fact that x(V) is a strictly decreasing function on \([0, V_{hit}]\); no special assumptions on \(\alpha , \beta \) needed so far.

Now consider the term \(H(V^*) m \frac{dx}{dV}(V^*)\) of the third line of (13). This term is always negative. For local stability of \(V^*\) we need

$$\begin{aligned} -2 \le H(V^*) m \frac{dx}{dV}(V^*), \end{aligned}$$
(14)

and for monotone approach we need that

$$\begin{aligned} -1 \le H(V^*) m \frac{dx}{dV}(V^*). \end{aligned}$$
(15)

It remains to investigate bounds on \(\frac{dx}{dV}(V)\). This involves fairly straightforward but somewhat tedious calculations; for details see Subsection 4.2 of Xin et al. (2018). These calculations show that in our model \(\alpha \ge 2\) is a sufficient condition for local stability of \(V^*\) and \(\alpha \ge 4\) is a sufficient condition for monotone approach to \(V^*\).

If we extend our model to allow for switching probabilities of the form (16), then these conditions are no longer sufficient for large enough values of the parameter \(\nu \), but for any given choice of the parameters \(\nu \) and \(\gamma \), local stability of \(V^*\) and even monotone approach to \(V^*\) will still be guaranteed for all sufficiently large \(\alpha \).

The conditions \(\alpha \ge 2\) for local stability of \(V^*\) and \(\alpha \ge 4\) for monotone approach are only sufficient in our model, but not necessary. However, some conditions on the parameters are needed even for stability. We numerically explored the value of \(H(V^*)m\frac{dx}{dV}(V^*)\) and found regions of the parameter space where (15) and even (14) fail. In order to better visualize these regions, we color-coded them by distinguishing values in certain relevant intervals. More specifically, we defined a function \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) \) in the following way:

  • \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) = -4\) if \(H(V^*)m\frac{dx}{dV}(V^*)< -2\),

  • \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) = -1.5\) if \(-2 \le H(V^*)m\frac{dx}{dV}(V^*) < -1.05\),

  • \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) = -0.7\) if \(-1.05 \le H(V^*)m\frac{dx}{dV}(V^*) < -1\),

  • \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) = -0.5\) if \(-1 \le H(V^*)m\frac{dx}{dV}(V^*) < 0\).

In Figs. 2 and 3 below, we plot the resulting color-coded partitions of the parameter space, together with some sample trajectories of \(V_n\) in regions of interest. Additional details about these observations are given in Xin et al. (2018). In particular, numerical explorations of the phase transitions for stability of \(V^*\) and monotone approach to this equilibrium are reported in Subsubsection 5.1.2 of Xin et al. (2018).

Fig. 2
figure 2

Numerical exploration of the stability of the interior equilibrium when \(\mathcal {R}_0=12\). Given that stability is determined by the quantity defined by \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) \) in (14) and (15), parameter sets in the blue/green region guarantee stability (as illustrated by simulations with \(\alpha = \beta =1\) and \(\alpha =1, \beta =5\)) and that parameter sets in the yellow region guarantee monotone convergence as well (see simulation with \(\alpha =2, \beta =1\)) (colour figure online)

Full size image
Fig. 3
figure 3

Numerical exploration of the stability of the interior equilibrium when \(\mathcal {R}_0=13\). Given that stability is determined by quantity defined by \(Int\left( H(V^*)m\frac{dx}{dV}(V^*)\right) \) in (14) and (15), parameter sets in the blue/green region (see \(\alpha =1, \beta =2\) simulation) guarantee stability, parameter sets in the dark blue region (see \(\alpha = 1, \beta =3\) simulation) do not guarantee stability and parameter sets in the yellow region (see \(\alpha = 2, \beta =2\) simulation) guarantee monotone convergence (colour figure online)

Full size image

4 Discussion: generalized Fermi functions and the interpretation of \(\alpha \) as a degree of open-mindedness

In the introduction we mentioned three functional forms for the probability of focal player i switching to the strategy of player j. Here we will discuss in more detail the relationship between these forms and the roles of their parameters. For convenience, let us repeat their formulas here. A fairly general form of smoothed imitation is given by

$$\begin{aligned} q_{switch}(i\rightarrow j) = \mu + \frac{\nu }{\alpha + e^{- \beta (C(i) - C(j)) + \gamma }}, \end{aligned}$$
(16)

where the parameters satisfy the inequalities \(\alpha , \nu > 0\), \(\mu , \beta \ge 0\).

In this paper, we focus on the case \(\mu = \gamma = 0\) and \(\nu = 1\) with \(\alpha \ge 1\), so that:

$$\begin{aligned} p_{switch}(i\rightarrow j) = \frac{1}{\alpha + e^{-\beta (C(i) - C(j))}}. \end{aligned}$$
(17)

For \(\alpha = 1\), we recover the classical Fermi functions

$$\begin{aligned} p_{switch}(i\rightarrow j) = \frac{1}{1 + e^{-\beta (C(i) - C(j))}}. \end{aligned}$$
(18)

The functional forms (16) and (17) are equivalent in the sense that they exhibit the same directions of change and therefore produce the same equilibria with rescaled \(\alpha \) and the same \(\beta \). More precisely, we make the following observation:

Lemma 3

Consider two models with the same parameters \(\beta , \mathcal {R}_0\) and with parameters \(\alpha '\) and  \(\alpha \), respectively, that satisfy \(\alpha ' = e^\gamma \alpha \) and \(\alpha \ge 1\). Assume that these models are constructed as the one described here, but in the first model the switching probabilities \(q_{switch}(i\rightarrow j)\) are given by (16) for \(\alpha '\), while in the second model the switching probabilities \(p_{switch}(i\rightarrow j)\) are given by (17) for \(\alpha \). Then for any given vaccination coverage \(V_n\), the sign of \(V_{n+1} - V_n\) in the second model will be the same as the sign of \(V_{n+1} - V_n\) in the first model.

Proof

Let \(V_n\) be the vaccination coverage for season n, and let \(\varDelta (n)\) denote the change of vaccination coverage from season n to season \(n + 1\). Let us consider the second model first, as it is the simpler one. In this model the switching probabilities are given by (17), so that

$$\begin{aligned} \begin{aligned} \varDelta _{model 2}(n)&= V_{n+1} - V_n = (1-V_n)V_n sp, \ \text{ where }\\ sp&= \frac{x_n}{\alpha +e^{-\beta (c_i-c_v)}} + \frac{1-x_n}{\alpha +e^{\beta c_v}} - \frac{x_n}{\alpha +e^{-\beta (c_v-c_i)}} - \frac{1-x_n}{\alpha +e^{-\beta c_v}}. \end{aligned} \end{aligned}$$

On the other hand, in the first model where the switching probabilities are given by (16) with \(\alpha \) replaced by \(\alpha ' = e^\gamma \alpha \),

$$\begin{aligned} \begin{aligned} \varDelta _{model 1}(n)&= V_{n+1} - V_n = (1-V_n)V_nsq,\ \ \text{ where } \\ sq&= x_n\left( \mu + \frac{\nu }{\alpha ' + e^{-\beta (c_i-c_v) +\gamma }}\right) + (1-x_n)\left( \mu +\frac{\nu }{\alpha ' + e^{\beta c_v + \gamma }}\right) \\&\quad - x_n\left( \mu +\frac{\nu }{\alpha '+e^{\beta (c_i-c_v)+\gamma }} \right) - (1-x_n)\left( \mu +\frac{\nu }{\alpha '+e^{-\beta c_v +\gamma }}\right) \\&= \nu \left( \frac{1}{e^{\gamma }}\right) \, sp. \end{aligned} \end{aligned}$$

Since \(\nu > 0\), the sign of \(V_{n+1}-V_n\) in the two models is always the same. \(\square \)

The parameter \(\beta \) plays similar roles in our generalized Fermi functions (17) as in the classical version (18). As the upper panel of Fig. 4 shows, when \(\beta \) increases, the function \(p_{switch}(i\rightarrow j)\) becomes more like a binary switch that reacts to the sign of the cost difference \(C(j) - C(i)\).

The parameter \(\alpha \) has two effects. The first is that high values of \(\alpha \) will make it less likely for the focal player to switch to the other strategy. See the middle panel of Fig. 4 for an illustration. In particular, high values of \(\alpha \) will have a stabilizing effect on the interior equilibrium of our model; see Sect. 3.5. Less frequent imitation could also be achieved by choosing low values of \(\nu \), and in view of Lemma 3 we can see that the frequency of imitation alone does not change the location of the interior equilibrium. The main result of this paper, that imitation with our generalized Fermi functions (17) for sufficiently large \(\alpha \) will give equilibria that are arbitrarily close to \(V_{hit}\), can be explained only by the second effect of \(\alpha \), which is a more subtle one. To illustrate this second effect, let us think of a two-step process for making the decision to imitate. In the first step, the focal player would make a decision on whether to consider imitating another player (with probability \(p_{imitate} = \alpha ^{-1}\)) or simply do the same as in the last season (with probability \(1 - \alpha ^{-1}\)). In the first case, the focal player i would then compare costs C(i) and C(j) for one randomly chosen player j and switch with conditional probability

$$\begin{aligned} p_{switch \, | \, imitate}(i\rightarrow j) = \frac{\alpha }{\alpha + e^{-\beta (C(i) - C(j))}}. \end{aligned}$$
(19)

This two-step-procedure is equivalent to the one-step decision given by (17). The lower panel of Fig. 4 shows how the function given by (19) depends on \(\alpha \). We can see that for large \(\alpha \), the focal player is very likely to switch to the other strategy once a decision to consider imitating has been made, even if that other strategy might be slightly worse than the focal player’s current strategy.

It is interesting to note that in the limit of \(\alpha \rightarrow \infty \), the absolute switching probability given by (17) approaches 0 while the conditional switching probability 19 will approach 1 for each fixed setting of the cost parameters. This shows that Theorem 1 is due to a balance between these two effects of increasing \(\alpha \), as neither can produce the result by itself.

We believe that high values of \(\alpha \) can be thought of as representing open-mindedness, understood as a willingness to experiment with new strategies unless there is strong evidence that they are not working well. While open-mindedness can be defined as willingness to consider ideas and opinions that are new or different to your own (Cambridge Academic Content Dictionary 2017), our working definition is more closely tied to actions than to attitudes. Note that our parameter \(\alpha \) has a straightforward translation into probabilities of adopting new courses of action. Probabilities of taking certain courses of action might be more objectively measurable in an experimental setting than attitudes. However, it remains an interesting open question beyond the expertise of the authors how to design actual experimental protocols for estimating these probabilities.

Fig. 4
figure 4

Dependence of the switching probabilities on the parameters and on the cost difference. Upper left panel: \(p_{switch}\) for \(\alpha = 1.2\). Upper right panel: \(p_{switch}\) for \(\beta = 0.5\). Lower panel: \(p_{switch\, | \, imitate}\) for \(\beta = 0.5\)

Full size image

Let us also mention that the same effect of open-mindedness can be achieved by considering large negative values for \(\gamma \) in (16). To see this, note that

$$\begin{aligned} \begin{aligned} q_{switch}(i\rightarrow j)&= \mu + \frac{\nu }{\alpha + e^{- \beta (C(i) - C(j)) + \gamma }} \\&= \mu + \frac{\nu e^{-\gamma }}{\alpha e^{-\gamma } + e^{- \beta (C(i) - C(j))}} \end{aligned} \end{aligned}$$

and recall that by Lemma 3 the parameters \(\mu \) and \(\nu e^{-\gamma }\) do not influence the locations of the equilibria in our models.

In Zhang et al. (2011), functional forms for the switching probabilities as in (16) for \(\mu = 0\) and \(\nu = \alpha = 1\) were considered. The authors suggested that positive values of \(\gamma \) represent inertia and negative values of \(\gamma \) represent “eagerness to switch.” In the context of (16), the relation between \(\gamma \) and inertia is not straightforward, as low values of both \(\mu \) and \(\nu \) also give small overall switching probabilities. However, the above calculations do show a direct correspondence between high values of \(\alpha \) and negative values of \(\gamma \), with both of them having a plausible interpretation in terms of open-mindedness.

5 Summary and directions of future work

The problem of designing effective public policy for inducing people to vaccinate against vaccine-preventable diseases is of great societal urgency (Manfredi and d’Onofrio 2013). To solve this problem, we need to understand how people really make vaccination decisions and how certain factors that enter this process influence the outcome.

The work presented here focuses on the role of imitation of successful others, which has been well-documented to be an important component of human decision-making (Apesteguia et al. 2007). Our model builds on the version of the model in Fu et al. (2011a) that assumed uniform mixing in the population. The only difference is including an additional parameter \(\alpha \) in the functional form for the probability of switching to another strategy. This parameter can be loosely interpreted as a degree of open-mindedness and has similar effects as some of the parameters in functional forms of these probabilities that were empirically derived in Traulsen et al. (2010).

Our results demonstrate that for sufficiently high values of \(\alpha \) the predicted equilibrium coverage will be arbitrarily close to the societally optimal value \(V_{hit}\) that gives herd immunity. They were independently confirmed in three different ways: analytically for large regions of the parameter space by way of proving Theorem 1, by numerical explorations of the equilibria predicted in the proof of this theorem [see Subsection 5.2 of Xin et al. (2018)] and by studying the equilibria that are being approached in simulated evolution of the vaccination coverages (see Sect. 3.4). These findings stand in stark contrast to the results reported in Fu et al. (2011a) for the standard Fermi functions (i.e. \(\alpha =1\)) in which imitation leads to vaccination coverage even below the Nash equilibrium when the cost of vaccination is small relative to the cost of infection and uniform mixing is assumed. As we deliberately kept our model as basic as possible and excluded all other factors, such as community structure, incentives, or misperceptions, these radical differences in the predictions can be due only to choosing a high value of \(\alpha \). We conclude that “open-minded imitation,” based on switching probabilities of the form (3) with suitable choices of \(\alpha \), provides a possible avenue for attenuating the vaccination dilemma.

The findings presented here open up a number of avenues for future research, in at least four directions. First, from the purely mathematical point of view it may be of some interest to derive tighter bounds in Theorem 1 and also to find complete characterizations of the local and global asymptotic stability of the interior equilibrium \(V^*\) that refine the partial results of Sect. 3.5.

Our model, like all models of this kind and in particular, our precursor Fu et al. (2011a), makes a number of simplifying assumptions about disease transmission and human decision-making. While many of these assumptions presumably don’t significantly alter the predictions, at least at the qualitative level, we have demonstrated here that the particular form of the smoothing function for the switching probabilities can matter a lot. Thus the second direction would be to examine the effect of the parameter \(\alpha \) in extensions of our model that incorporate a number of additional aspects of voluntary vaccination dynamics that have been deliberately set aside here. We have already done some preliminary work on extending our model to the important and biologically more realistic case when the vaccine has an efficacy of less than 100% and found a similar pattern as the one reported here in that equilibrium coverage can get arbitrarily close to the societal optimum. However, the dependence of the optimal choice of \(\alpha \) on \(\beta \) is more complicated than in Theorem 1 and still needs to be worked out in more detail. Other aspects that might be incorporated into more detailed versions of our model are restrictions of disease transmission and/or imitation to edges of contact networks, as has already been studied for the case of classical Fermi functions in Fu et al. (2011a) and a number of related papers; see Wang et al. (2015) for a review. Similarly, one can study the effects of a mixture of rational decision-making and imitation (d’Onofrio et al. 2012; Ndeffo Mbah et al. 2012), incentives (Li et al. 2017; Zhang et al. 2013, 2014), misperceptions of costs (Bauch and Bhattacharyya 2012; Coelho and Codeço 2009; Cojocaru et al. 2007; Reluga et al. 2006; Voinson et al. 2015), altruism (Shim et al. 2012; Szolnoki et al. 2012; Zhang 2013), peer pressure (Ichinose and Kurisaku 2017; Wu and Zhang 2013), presence of individuals who remain committed to vaccinating or not vaccinating without ever imitating others (Fukuda and Tanimoto 214; Liu et al. 2012), the effects of other available control measures or treatment options (Andrews and Bauch 2015; Chen 2006; Ida and Tanimoto 2018; Jijón et al. 2017; Li et al. 2017; Wang et al. 2014), or variability of \(\mathcal {R}_0\) from season to season. Going further than varying \(\mathcal {R}_0\) from season to season, we could even model the seasonality of influenza continuously throughout the year rather than our current on/off model for flu season (Buonomo et al. 2018). One can also study the effects of our parameter \(\alpha \) when imitation is based on comparison with the average cost of a larger sample of other individuals (Fukuda et al. 2014; Ichinose and Kurisaku 2017; Iwamura and Tanimoto 2018; Li et al. 2017), or when weighted averages of costs over a number of previous seasons are being compared (Zhang et al. 2012).

The third direction would be to study applications of our functional form (3) to domains other than vaccination games. The paper Zhang et al. (2011) considers prisoner’s dilemma games in a finite population with switching probabilities of the form \(\frac{1}{1 + e^{- \beta (C(i) - C(j)) - \tau _s}}\), where \(\tau _s\) depends on the strategy of the focal player. Since \(-\tau _s\) is the same as \(\gamma \) in our notation, in view of our earlier observations we can consider these as rescaled versions of (3). The authors of Zhang et al. (2011) obtained a number of analytical results, but their focus is on fixation probabilities, time to fixation, and stochastic stability of the equilibria, which is different from ours. The literature on applications of evolutionary games with the structure of a prisoner’s dilemma is vast, and vaccination games are only a small part of it. So it seems likely that our version of (3) or its counterpart in Zhang et al. (2011) could find many applications outside of vaccination games. Let us also remark that there may be some applications to evolutionary computation (Eiben and Smith 2015) as well. In this field the balance between exploration and exploitation is of paramount importance, and our interpretation of high \(\alpha \) as open-mindedness bears some resemblance to shifting this balance towards the former for high values of this parameter.

Finally, it would be important to get a better understanding of how our parameter \(\alpha \) relates to actual decision-making by real people. If the benefits of open-minded imitation in the sense of high values of this parameter can also be demonstrated in more detailed and realistic models than the one studied here, then it will be of interest to explore how public policy could enhance more “open-minded” decision-making. This third direction of suggested follow-up work would require a multidisciplinary effort that goes far beyond the realm of mathematical modeling.