The Effects of Imitation Dynamics on Vaccination Behaviours in SIR-Network Model

We present a series of SIR-network models, extended with a game-theoretic treatment of imitation dynamics which result from regular population mobility across residential and work areas and the ensuing interactions. Each considered SIR-network model captures a class of vaccination behaviours influenced by epidemic characteristics, interaction topology, and imitation dynamics. Our focus is the resultant vaccination coverage, produced under voluntary vaccination schemes, in response to these varying factors. Using the next generation matrix method, we analytically derive and compare expressions for the basic reproduction number R0 for the proposed SIR-network models. Furthermore, we simulate the epidemic dynamics over time for the considered models, and show that if individuals are sufficiently responsive towards the changes in the disease prevalence, then the more expansive travelling patterns encourage convergence to the endemic, mixed equilibria. On the contrary, if individuals are insensitive to changes in the disease prevalence, we find that they tend to remain unvaccinated. Our results concur with earlier studies in showing that residents from highly connected residential areas are more likely to get vaccinated. We also show that the existence of the individuals committed to receiving vaccination reduces R0 and delays the disease prevalence, and thus is essential to containing epidemics.


Introduction
Vaccination has long been established as a powerful tool in managing and controlling infectious diseases by providing protection to susceptible individuals [1]. With a sufficiently high vaccination coverage, the probability of the remaining unvaccinated individuals getting infected reduces significantly. However, such systematic programs by necessity may limit the freedom of choice of individuals. When vaccination programs are made voluntary, the vaccination uptake declines as a result of individuals choosing not to vaccinate, as seen in Britain in 2003 when the vaccination program for Measles-Mumps-Rubella (MMR) was made voluntary [2]. Parents feared possible complications from vaccination [2,3] and hoped to exploit the 'herd immunity' by assuming other parents would choose to vaccinate their children. Such hopes did not materialize precisely because other parents also thought similarly.
Under a voluntary vaccination policy an individual's decision depends on several factors: the social influence from one's social network, the risk perception of vaccination, and the risk perception of infection, in terms of both likelihood and impact [4]. This decision-making is often modelled using game theory [3,[5][6][7][8][9][10][11], by allowing individuals to compare the cost of vaccination and the potential cost of non-vaccination (in terms of the likelihood and impact of infection), and adopting imitation dynamics in modelling the influence of social interactions. However, many of the earlier studies [12][13][14] were based on the key assumption of a well-mixed homogeneous population where each individual is assumed to have an equal chance of making contact with any other individual in the population. This population assumption is rather unrealistic as large populations are often diverse with varying levels of interactions. To address this, more recent studies [7][8][9][10][11]15] model populations as complex networks where each individual, represented by a node, has a finite set of contacts, represented by links [16].
It has been shown that there is a critical cost threshold in the vaccination cost above which the likelihood of vaccination drops steeply [6]. In addition, highly connected individuals were shown to be more likely to choose to be vaccinated as they perceive themselves to be at a greater risk of being infected due to high exposure within the community. In modelling large-scale epidemics, the population size (i.e., number of nodes) can easily reach millions of individuals, interacting in a complex way. The challenge, therefore, is to extend epidemic modelling not only with the individual vaccination decision-making, but also capture diverse interaction patterns encoded within a network. Such an integration has not yet been formalized, motivating our study. Furthermore, once an integrated model is developed, a specific challenge is to consider how the vaccination imitation dynamics developing across a network affects the basic reproduction number, R 0 . This question forms our second main objective.
One relevant approach partially addressing this objective is offered by the multi-suburb (or multi-city) SIR-network model [17][18][19][20][21] where each node represents a neighbourhood (or city) with a certain number of residents. The daily commute of individuals between two neighbourhoods is modelled along the network link connecting the two nodes, which allows to quantify the disease spread between individuals from different neighbourhoods. Ultimately, this model captures meta-population dynamics in a multi-suburb setting affected by an epidemic spread at a greater scale. To date, these models have not yet considered intervention (e.g., vaccination) options, and the corresponding social interaction across the populations.
To incorporate the imitation dynamics, modelled game-theoretically, within a multi-suburb model representing mobility, we propose a series of integrated vaccination-focused SIR-network models. This allows us to systematically analyze how travelling patterns affect the voluntary vaccination uptake due to adoption of different imitation choices, in a large distributed population. The developed models use an increasingly complex set of vaccination strategies. Thus, our specific contribution is the study of the vaccination uptake, driven by imitation dynamics under a voluntary vaccination scheme, using an SIR model on a complex network representing a multi-suburb environment, within which the individuals commute between residential and work areas.
In Section 2.2, we present the models and methods associated with this study: in particular, we present analytical derivations of the basic reproduction number R 0 for the proposed models using the Next Generation Operator Approach [22], and carry out a comparative analysis of R 0 across these models. In Section 3, we simulate the epidemic and vaccination dynamics over time using the proposed models in different network settings, including a pilot case of a 3-node network, a 3000-node Erdös-Rényi random network and a real-world commuting network in Greater Sydney area generated from 2016 Australian census data. The comparison of the produced results across different models and settings is carried out for the larger network, with the focus on the emergent attractor dynamics, in terms of the proportion of vaccinated individuals. Particularly, we analyze the sensitivity of the individual strategies (whether to vaccinate or not) to the levels of disease prevalence produced by the different considered models. Section 4 concludes the study with a brief discussion of the importance of these results.

Basic Reproduction Number R 0
The basic reproduction number R 0 is defined as the number of secondary infections produced by an infected individual in an otherwise completely susceptible population [23]. It is well-known that R 0 is an epidemic threshold, with the disease dying out as R 0 < 1, or becoming endemic as R 0 > 1 [17,24]. This finding strictly holds only in deterministic models with infinite population [25]. The topology of the underlying contact network is known to affect the epidemic threshold [26]. Many disease transmission models have shown important correlations between R 0 and the key epidemic characteristics (e.g., disease prevalence, attack rates, etc.) [27]. In addition, R 0 has been considered as a critical threshold for phase transitions studied with methods of statistical physics or information theory [28].

Vaccination Model with Imitation Dynamics
Imitation dynamics, a process by which individuals copy the strategy of other individuals, is widely used to model vaccinating behaviours incorporated with SIR models. The model proposed in [3] applied game theory to represent parents' decision-making about whether to get their newborns vaccinated against childhood disease (e.g., measles, mumps, rubella, pertussis). In this model, individuals are in a homogeneously mixing population, and susceptible individuals have two 'pure strategies' regarding vaccination: to vaccinate or not to vaccinate. The non-vaccination decision can change to the vaccination decision at a particular sampling rate, however the vaccination decision cannot be changed to non-vaccination. Individuals adopt one of these strategies by weighing up their perceived payoffs, measured by the probability of morbidity from vaccination, and the risk of infection respectively. The payoff for vaccination ( f v ) is given as, and the payoff for non-vaccination ( f nv ), measured as the risk of infection, is given as where r v is the perceived risk of morbidity from vaccination, r nv is the perceived risk of morbidity from non-vaccination (i.e., infection), I(t) is the current disease prevalence in population fraction at time t, and m is the sensitivity to disease prevalence [3]. From Equations (1) and (2), it can be seen that the payoff for vaccinated individuals is a simple constant, and the payoff for unvaccinated individuals is proportional to the extent of epidemic prevalence (i.e., the severity of the disease).
It is assumed that life-long immunity is granted with effective vaccination. If one individual decides to vaccinate, s/he cannot revert back to the unvaccinated status. To make an individual switch to a vaccinating strategy, the payoff gain, ∆E = f v − f nv must be positive (∆E > 0). Let x denote the relative proportion of vaccinated individuals, and assume that an unvaccinated individual can sample a strategy from a vaccinated individual at certain rate σ, and can switch to the 'vaccinate strategy' with probability ρ∆E. Then the time evolution of x can be defined as: where δ = σρ can be interpreted as the combined imitation rate that individuals use to sample and imitate strategies of other individuals. Equation (3) can be rewritten so that x only depends on two parameters: where κ = δr v and ω = mr nv /r v , assuming that the risk of morbidity from vaccination, r v , and the risk of morbidity from non-vaccination (i.e., morbidity from infection), r nv , are constant during the course of epidemics. Hence, κ is the adjusted imitation rate, and ω measures individual's responsiveness towards the changes in disease prevalence. The disease prevalence I is determined from a simple SIR model [23] with birth and death that divides the population into three health compartments: S (susceptible), I (infected) and R (recovered). If a susceptible individual encounters an infected individual, s/he contracts the infection with transmission rate β, and thus progresses to the infected compartment, while an infected individual recovers at rate γ. It is assumed that the birth and death rate are equal, denoted by µ [29]. Individuals in all compartments die at an equal rate, and all newborns are added to the susceptible compartment unless and until they are vaccinated. Vaccinated newborns are moved to the recovered class directly, with a rate µx. Therefore, the dynamics are modelled as follows [3]: where S, I, and R represent the proportion of susceptible, infected, and recovered individuals respectively (that is, S + I + R = 1), and x represents the proportion of vaccinated individuals within the susceptible class at a given time. The variables S, I, and R, as well as x, are time-dependent, but for simplicity, we omit subscript t for these and other time-dependent variables.
Model (5) predicts oscillations in vaccine uptake in response to changes in disease prevalence. It is found that oscillations are more likely to occur when individuals imitate each other more quickly (i.e., higher κ). The oscillations are observed to be more volatile when people alter their vaccinating behaviours promptly in response to changes in disease prevalence (i.e., higher ω). Overall, higher κ or ω produce stable limit cycles at greater amplitude. Conversely, when individuals are insensitive to changes in disease prevalence (i.e., low ω), or imitate at a slower rate (i.e., low κ), the resultant vaccinating dynamics converge to equilibrium [3].

Multi-City Epidemic Model
Several multi-city epidemic models [17][18][19] consider a network of suburbs as M nodes, in which each node i ∈ V represents a suburb, where i = 1, 2, · · · M. If two nodes i, j ∈ V are linked, a fraction of population living in node i can travel to node j and back (commute from i to j, for example for work) on a daily basis. The connectivity of suburbs and the fraction of people commuting between them are represented by the population flux matrix φ whose entries represent the fraction of population daily commuting from i to j, φ ij ∈ [0, 1] (Equation (6)). Note that φ ij = φ ji , thus φ is not a symmetric matrix. Figure 1 shows an example of such travel dynamics. Trivially, each row in φ, representing proportions of the population of a suburb commuting to various destinations, has to sum up to unity: However, the column sum measures the population influx to a suburb during a day, and therefore depends on the node's connectivity. Column summation in φ does not necessarily equal to unity.
It is also necessary to differentiate between present population N p j and native population N j in node j on a particular day. N p j is used as a normalizing factor to account for the differences in population flux for different nodes: Assuming the disease transmission parameters are identical in all cities, the standard incidence can be expressed as [17]: where I l and S l denote respectively the number of infective and susceptible individuals in city l, and λ l is the average number of contacts in city l per unit time. Incorporating travelling pattern defined by φ, yields [19]: The double summation term in this expression captures the infection at suburb j due to the encounters between the residents from suburb i and the residents from suburb k (k could be any suburb including i or j) occurring at suburb j, provided that suburbs i, j, k are connected with non-zero population flux entries in φ.
R 0 can be derived using the Next Generation Approach (see Appendix A for more details). For example, for a multi-city SEIRS model with 4 compartments (susceptible, exposed, infective, recovered), introduced in [17], and a special case where the contact rate λ j is set to 1, while β is identical across all cities, R 0 has the following analytical solution: where 1/d, 1/ , and 1/γ denote the average lifetime, exposed period (i.e., the latent period of individuals becoming infectious) and infective period respectively. It can be seen that this solution concurs with a classical SEIRS model with no mobility. When → ∞, Equation (11)  In this study, we only consider vaccinations that confer lifelong immunity, mostly related to childhood diseases such as measles, mumps, rubella and pertussis. Vaccinations against such diseases are often administered for individuals at a young age, implying that newborns (more precisely, their parents) often face the vaccination decision. These are captured by the first extension above. However, during an outbreak, adults may face the vaccination decision themselves if they were not previously vaccinated. For example, diseases which are not included in formal childhood vaccination programs, such as smallpox [30], may present vaccination decisions to the entire susceptible population. These are captured by the second extension. The third extension introduces a small fraction of susceptible population as committed vaccine recipients, i.e., those who would always choose to vaccinate regardless of the prevalence. The purpose of modelling committed vaccine recipients is to demonstrate how successful immunization education campaigns could affect vaccination dynamics.
Our models divide the population into many homogeneous groups [31], based on their residential suburbs. Within each suburb (i.e., node), residents are treated as a homogeneous population. It is assumed that the total population within each node is conserved over time. The dynamics of epidemic and vaccination at each node, are referred to as 'local dynamics'. The aggregate epidemic dynamics of the entire network can be obtained by summing over all nodes, producing 'global dynamics'.
We model the 'imitation dynamics' based on the individual's travelling pattern and the connectivity of their node defined by Equation (6). For any node i ∈ V, let x i denote the fraction of vaccinated individuals in the susceptible class in suburb i. On a particular day, unvaccinated susceptible individuals (1 − x i ) commute to suburb j and encounter vaccinated individuals from node k (where k may be i itself, j or any other nodes) and imitate their strategy. However, this 'imitation' is only applicable in the case of a non-vaccinated individual imitating the strategy of a vaccinated individual (that is, deciding to vaccinate), since the opposite 'imitation' cannot occur. Therefore, in our model, every time a non-vaccinated person from i comes in contact with a vaccinated person from k, they imitate the 'vaccinate' strategy if the perceived payoff outweighs the non-vaccination strategy.
Following Equation (4) proposed in [3], the rate of change of the proportion of vaccinated individuals in i, that is, x i , over time can be expressed by: The players of the vaccination game are the parents, deciding whether or not to vaccinate their children using the information of the disease prevalence collected from their daily commute. For example, a susceptible individual residing at node i and working at node j uses the local disease prevalence at node j to decide whether to vaccinate or not. If such a susceptible individual is infected, that individual will be counted towards the local epidemic prevalence at node i.
We measure each health compartment as a proportion of the population. Hence, we define a ratio, p l , as the ratio between present population N p l and the 'native' population N l in node l on a particular day as: Our main focus is to investigate the effects of vaccinating behaviours on the global epidemic dynamics. To do so, we vary three parameters: • Individual's responsiveness to changes in disease prevalence, ω • Adjusted imitation rate, κ • Vaccination failure rate, ζ while κ and ω have had been considered in [3], we introduce a new parameter ζ as the vaccination failure rate, ζ ∈ [0, 1] to consider the cases where the vaccination may not be fully effective. The imitation component (Equation (12)) is common in all model extensions. While the epidemic compartments vary depending on the specific extension, Equation (12) is used consistently to model the relative rate of change in vaccination behaviours.

Vaccination Available to Newborns Only
Model (14) captures the scenario when vaccination opportunities are provided to newborns only: Unvaccinated newborns and newborns with unsuccessful vaccination µ[ζx i + (1 − x i )] stay in the susceptible class. Successfully vaccinated newborns, on the other hand, move to the recovered class µ(1 − ζ)x i . Other population dynamics across health compartments follow the model (5).
Since S + I + R = 1(Ṡ +İ +Ṙ = 0), model (14) can be reduced to: Model (15) has a disease-free equilibrium (disease-free initial condition) We can now obtain R 0 for the global dynamics in this model by using the Next Generation Approach, where R 0 is given by the most dominant eigenvalue (or 'spectral radius' ρ) of FV −1 , where F and V are M × M matrices, representing the 'new infections' and 'cases removed or transferred from the infected class', respectively in the disease free condition [22,32]. As a result, R 0 is determined as follows (see Appendix A for detailed derivation): where Clearly, the solution of R 0 depends on S 0 i , I 0 k , epidemiological parameters (i.e., β, γ, µ), and mobility matrix φ. Here, we obtain an analytical expression of R 0 by studying a special case where all nodes are uniform at disease free equilibrium. Proposition 1. At the disease free equilibrium, x 0 has two solutions (x 0 = 0 or x 0 = 1) if x 0 and S 0 are uniform across all nodes (see proof in Appendix B).
Using Proposition 1, F can be simplified as: We can now prove the following simple but useful proposition (see proof in Appendix B).

Proposition 2.
Matrix G, defined as follows: is a Markov matrix.
As a Markov matrix, G always has the most dominant eigenvalue of unity. All other eigenvalues are smaller than unity in absolute value [33]. Now the next generation matrix K can be obtained as follows: From Proposition 2 and Equation (17), R 0 can be obtained as: Noting Proposition 1, there are two cases: x = 0 and x = 1. When nobody vaccinates (x 0 = 0,S 0 = 1), the entire population remains susceptible, which reduces model (15) to a canonical SIR model without vaccination intervention and R 0 returns to β γ+µ . On the other hand, if the whole population is vaccinated (x 0 = 1), the fraction of susceptible population only depends on the vaccine failure rate ζ, resulting in: Clearly, fully effective vaccination (ζ = 0) would prohibit disease spread (R 0 = 0). Partially effective vaccination could potentially suppress disease transmission, or even eradicate disease spread if R 0 < 1. If all vaccinations fail (ζ = 1), model (15) concurs with a canonical SIR model without vaccination, which also corresponds to a special case in Equation (11) where → ∞.

Vaccination Available to the Entire Susceptible Class
If the vaccination opportunity is expanded to the entire susceptible class, including newborns and adults, the following model is proposed based on model (15): Model (26) has a disease-free equilibrium (S 0 Substituting S 0 i using Equation (16) yields where I is a M × M identity matrix. The next generation matrix K can then be obtained as: Using Proposition 2, R 0 can be obtained as: When nobody vaccinates (x 0 = 0), model (26) concurs with a canonical SIR model with R 0 reducing to β γ+µ . When the vaccination attains the full coverage in the population (x 0 = 1), the magnitude of R 0 depends on the vaccine failure rate ζ: Equation (32) Given that µ 1, R 0 is well below the critical threshold: If 0 < ζ < 1, by comparing Equations (25) and (32), we note that R 0 becomes smaller, provided ζ > µ, that is:

Vaccination Available to the Entire Susceptible Class with Committed Vaccine Recipients
We now consider the existence of committed vaccine recipients, x c , as a fraction of individuals who would choose to vaccinate regardless of payoff assessment [34,35] (0 < x c 1). We assume that committed vaccine recipients are also exposed to vaccination failure rate ζ and are distributed uniformly across all nodes. It is also important to point out that the fraction of committed vaccine recipients is constant over time. However, they can still affect vaccination decision for those who are not vaccinated, and consequently, contribute to the rate of change of the vaccinated fraction x. Model (26) can be further extended to reflect these considerations, as follows: Model (36) has a disease-free equilibrium: Substituting S 0 i , using Equation (37), yields: where I is a M × M identity matrix.
The next generation matrix K can now be obtained as: Proposition 3. At the disease free equilibrium, x 0 has two solutions x 0 = x c = 0 or x 0 + x c = 1 if x 0 , x c and S 0 are uniform across all nodes.
R 0 can be obtained for this case as: If nobody chooses to vaccinate (x 0 = x c = 0), R 0 can be reduced to β γ+µ . Conversely, if the entire population is vaccinated (x 0 + x c = 1), Equation (41) reduces to Equation (32), and R 0 is purely dependent on the vaccine failure rate ζ.

Model Parameterisation
The proposed models aim to simulate a scenario of a generic childhood disease (e.g., measles) where life-long full immunity is acquired after effective vaccination. The vaccination failure rate, ζ, is set as the probability of ineffective vaccination, to showcase the influence of unsuccessful vaccination on the global epidemic dynamics. In reality, the vaccination for Measles-Mumps-Rubella (MMR) is highly effective: for example, in Australia, an estimated 96% of vaccines that were administered are successful in conferring immunity [36].
The population flux matrix φ is derived from the network topology, in which each entry represents the connectivity between two nodes: if two nodes are not connected, φ ij = 0; if two nodes are connected, φ ij is randomly assigned in the range of (0, 1]. The population influx into a node i within a day is represented by the column sum ∑ M j=1 φ ji in flux matrix φ. The entries of φ are determined by network topologies. The parameters used for all simulations are summarized in Table 1. Same initial conditions are applied to all nodes. Parameters ω and κ are calibrated based on values used in [3]. Here, we aim to investigate how vaccination behaviours affect overall epidemic dynamics. To do so, we mainly study the effects of two parameters on all models: ω, responsiveness to changes in disease prevalence; and κ, imitation rate (i.e., how quickly individuals imitate each other). Our simulations were carried out on three networks: • a pilot case of a network with 3 nodes (suburbs), • an Erdös-Rényi random network [38] with 3000 nodes (suburbs), and • the commuting network in Greater Sydney with 311 suburbs (nodes) [39,40].
Each of the first two networks was used in conjunction with the following three models of vaccination behaviours: • vaccinating newborns only: using model (15) • vaccinating the susceptible class: using model (26), • vaccinating the susceptible class with committed vaccine recipients: using model (36).
The third network is a real-world network and is considered a more realistic representation of travelling patterns. It was only used in conjunction with model (15) where only newborns are vaccinated as the existing immunisation programmes in Australia typically administer measles vaccines to newborns.
In the pilot case of 3-node network, two network topologies are studied: an isolated 3-node network where all residents remain in their residential nodes without travelling to other nodes (equivalent to models proposed in [3]), and a fully connected network where the residents at each node commute to the other two nodes, and the population fractions commuting are symmetric and uniformly distributed (Figure 2).
We then consider a suburb-network modelled as an Erdös-Rényi random graph (with number of nodes M = 3000) with average degree k = 4, in order to study how an expansive travelling pattern affects the global epidemic dynamics. (See Figure 10 for the degree distribution of the network used.) Other topologies, such as lattice, scale-free [41,42], small-world networks and other real-world networks [20,43], can easily be substituted here, though in this study our focus is on an Erdös-Rényi random graphs.

Vaccinating Newborns Only
When there is no population mobility, we observe three distinctive equilibria (Figure 3; dotted lines): a pure non-vaccinating equilibrium (where x f = 0, representing the final condition, and ω = 1000), a mixed equilibrium (where x f = 0, ω = 2500), and stable limit cycles (where x f = 0, ω = 3500). This observation is in qualitative agreement with the vaccinating dynamics reported by [3].
When commuting is allowed, individuals commute to different nodes, and their decision will no longer rely on the single source of information (i.e., the disease prevalence in their residential node) but will also depend on the disease prevalence at their destination. As a consequence, the three distinctive equilibria are affected in different ways (Figure 3; solid lines). The amplitude of the stable limit cycles at ω = 3500 is reduced as a result of the reduced disease prevalence. It takes comparatively longer (compare the dotted lines and solid lines in Figure 3) to converge to the pure non-vaccinating equilibrium at ω = 1000, and to the endemic, mixed equilibrium at ω = 2500 with high amplitude of oscillation at the start of the epidemic spread. As ω can be interpreted as the responsiveness of vaccinating behaviour to the disease prevalence [3], if individuals are sufficiently responsive (i.e., ω is high), the overall epidemic is suppressed more due to population mobility (and the ensuing imitation), as evidenced by smaller prevalence peaks. Conversely, if individuals are insensitive to the prevalence change (i.e., ω is low), epidemic dynamics with equal population mobility may appear to be more volatile at the start, but the converged levels of both prevalence peak and vaccine uptake remain unchanged in comparison to the case where there is no population mobility. We further investigated how the vaccination failure rate ζ affects the global epidemic dynamics. If, for example, only a half of the vaccine administered is effective (ζ = 0.5), as shown in Figure 4b-d, the infection peaks arrive sooner, for all values of ω. As a result of having the earlier infection peaks, the individuals respond to the breakout and choose to vaccinate sooner, causing the vaccine uptake to rise. This seemingly counter-intuitive behaviours have also been reported in [44]. When ζ = 0, the prevalence peaks take longer to develop and the extended period gives individuals an illusion that there may not be an epidemic breakout, and consequently encourages 'free-riding' behaviour. In the case where half of the vaccinations fail, epidemic breaks out significantly earlier at lower peaks, an observation that is beneficial to encourage responsive individuals to choose to vaccinate. Although some (in this case half) of the vaccines fail, a sufficiently high vaccine uptake still curbs prevalence peaks and shortens the length of breakout period. In terms of the final vaccination uptake, the vaccine failure rate predominantly affects behaviours of responsive individuals (when ω is high, such as 3500). In this case, instead of the stable limit cycles observed when ζ = 0, an endemic, mixed equilibrium is reached when the vaccination failure rate is significant. Final vaccination uptake is impacted little by vaccination failure rate when individuals are insensitive to changes in disease prevalence (when ω is lower in value, such as 2500 or 1000). These observations are illustrated in Figure 4.

Vaccinating the Entire Susceptible Class
If (voluntary) vaccination is offered to the susceptible class regardless of the age, the initial condition of x = 0.95 represents the scenario that the vast majority of population are immunized to begin with, and the epidemic would not breakout until the false sense of security provided by the temporary 'herd immunity' settles in. This feature is observed in Figure 5a,b: the vaccination coverage continues to drop at the start of the epidemic breakout, indicating that individuals, regardless of their responsiveness towards the prevalence change, exploit the temporary herd immunity until an infection peak emerges. Since vaccination is available to the entire susceptible class, a small increase in x could help suppresses disease prevalence. However, when vaccination is partially effective (i.e., ζ = 0.5), the peaks in vaccine uptake, x, are no longer an accurate reflection of the actual vaccination coverage, and such peaks therefore may not be sufficient to adequately suppress infection peaks. It can also be observed that vaccinating the susceptible class encourages non-vaccinating behaviours due to the perception of herd immunity, and therefore pushes the epidemic towards an endemic equilibrium, particularly when sensitivity to prevalence is relatively high ( mid and high ω), as the responsive individuals would react promptly to the level of disease prevalence, altering their vaccination decisions. However, no substantial impact is observed on those individuals who are insensitive to changes in the disease prevalence ( low ω).

Erdös-Rényi Random Network of 3000 Nodes
We now present the simulation results on a much larger Erdö-Rényi random network of M = 3000 nodes, which more realistically reflects the size of a modern city and its commuting patterns. It was observed by previous studies [45] that such a larger system requires a higher vaccination coverage to achieve herd immunity, and thus curbs 'free-riding' behaviours more effectively.

Vaccinating Newborns Only
In this case, three distinct equilibria are observed (as shown in Figure 6a) for three values of ω: a pure non-vaccinating equilibrium at ω = 1000 and two endemic mixed equilibria at ω = 2500 and ω = 3500 respectively, replacing the stable limit cycles at high ω previously observed in the 3-node case. As expected, when individuals are more responsive to disease prevalence (higher ω), they are more likely to get vaccinated. The expansive travelling pattern also somewhat elevates the global vaccination coverage level (compared to the 3-node case), particularly in the case that the individuals show a moderate level of responsiveness to prevalence (ω = 2500), and shortens the convergence time to reach the equilibrium. However, for those who are insensitive to the changes in disease prevalence (ω = 1000), the more expansive commuting presented by the larger network does not affect either the level of voluntary vaccination, or the convergence time in global dynamics-these individuals remain unvaccinated as in the case of the smaller network. It is also found that the vaccination coverage is very sensitive to the disease prevalence change at the larger network since the disease prevalence peaks are notably lower than the 3-node case counterparts. If half of the vaccines administered are unsuccessful (ζ = 0.5), the impact of early peaks on global dynamics is magnified in a larger network (Figure 6b-d), leading to a shorter convergence time, although the final equilibria are hardly affected compared to the 3-node case as shown in Figure 4. It is observed that the shape of oscillation, in terms of amplitude and period, varies with different values of ω. Figure 6a (also Figure 4 for a much smaller network) show that lower ω is related with higher damping ratio in the oscillatory dynamics of the disease prevalence I and the vaccination coverage x. A closer inspection reveals that the oscillations observed are non-harmonic with the period reducing after every cycle (Figure 7a,b). Also, period is positively correlated with the value of ω, an observation supported by some preliminary analysis as reported in Appendix C, indicating that the time interval between prevalence peaks is longer if population is sufficiently responsive to prevalence change. Such a conclusion is not affected by the vaccine failure rate ζ, although the period is noticeably smaller at the start of the outbreak when ζ is larger (ζ = 0.5) (Figure 7c,d). Furthermore, it is worth noting that the oscillatory properties may also be affected by epidemiological parameters, such as transmission rate β, a shown in Appendix C. To verify the analytical derivation of R 0 , we investigated the relationship between epidemic and vaccination dynamics and the basic reproduction ratio R 0 . According to Equation (24), higher vaccination coverage leads to a reduction of R 0 . Figure 8 confirms this analytical dependency by showing that higher vaccination coverage (higher x), reduces R 0 and consequently leads to smaller cumulative prevalence, provided individuals are sufficiently responsive towards prevalence change (i.e., ω = 2500 and ω = 3500). When population is insensitive to prevalence change (ω = 1000), the vaccination dynamics always converge to the non-vaccination equilibrium, corresponding to the case where x = 0 in Equation (24). Cumulative prevalence I tot is obtained by integrating the prevalence over the simulated time frame. Note that different ω settings correspond to different ranges for R 0 due to the different vaccination coverage, x, at their respective endemic equilibria. Note that in (b) the case for ω = 1000 is not shown because it is trivially zero for all values of R 0 .
We also compared the epidemic dynamics in terms of the adjusted imitation rate, represented by the parameter κ, as shown in Figure 9 (recall that the value of imitation rate used in our simulations, unless otherwise stated, is κ = 0.001 as reported in Table 1). Here we present results where this value of imitation rate is compared with a much smaller value of κ = 0.00025. In populations where the imitation rate is comparatively high (i.e., higher κ), the oscillations in prevalence and vaccination dynamics appear earlier in time with larger amplitude, although the converged vaccine uptake and disease prevalence are relatively unaffected by the value of κ. For a higher κ, the convergence to equilibria is quicker. These observations are consistent with results reported by earlier studies [3]. A smaller value of κ (κ = 0.00025) also altered behaviours of those individuals who are insensitive to prevalence change. Instead of converging to the pure non-vaccination equilibrium, the dynamics converge to a mixed endemic state, meaning that when imitation rate is very low, some individuals would still choose to vaccinate even when they are not very sensitive to disease prevalence. For such a large network, vaccinating behaviour of individuals also depends on the weighted degree (number and weight of connections) of the suburb (node) in which they live. This is illustrated in Figure 10. For individuals living in highly connected suburbs (nodes), there are many commuting destinations, allowing access to a broad spectrum of information on local disease prevalence. Therefore, it is not surprising that we found that individuals living in 'hubs' are more likely to get vaccinated, particularly when the population is sensitive to prevalence (ω is high), as shown in Figure 10a. This observation is in accordance with previous studies [6]. Note that in in Figure 10a, nodes with degree k ≥ 10 are grouped into one bin to represent hubs, as the frequency count of these nodes is extremely low. Note also that there is a positive correlation between the number of degrees and the volume of population influx as measured by the sum of proportions from the source nodes (Figure 10b), indicating that a highly connected suburb has greater population influx on a daily basis.
To verify that the above mentioned simulation experiments reasonably model real-world scenarios, we performed additional experiments using a real-world commuting network of Greater Sydney generated from the 2016 Australian census data [39,40,46,47]. In contrast with agent-based models [46,47], in which commuting individuals interact at their workplaces or schools only with other commuters, rather than with the residential population, we consider a general case where commuting and residential populations are well-mixed at every node. This assumption provides a generic insight into how empirical networks affect epidemic and vaccination dynamics. As illustrated in Figure 11d, the real commuting network is substantially denser (M = 311, k ≈ 150) in comparison with the studied Erdös-Rényi random network. While the three distinctive equilibria still hold, it takes longer for insensitive individuals (i.e., ω = 1000) to reach non-vaccinating equilibrium (Figure 11b), an observation that is not previously noted in the 3-node pilot case and 3000-node Erdös-Rényi random network. In general, however, the dynamics produced on this real-world network concur with the results obtained for previously described simulation experiments.

Vaccinating Entire Susceptible Class
If (voluntary) vaccination is offered to the susceptible class regardless of age, we found that the global epidemic dynamics converge quicker compared to the similar scenario in the 3-node case. Only one predominant infection peak is observed, corresponding to the high infection prevalence around year 25, as shown in Figure 12. Oscillations of small magnitude are observed for both vaccine uptake and disease prevalence at later time-steps for middle or high values of ω (as shown in insets of Figure 12). These findings also largely hold if half of the vaccines administered are ineffective (i.e., ζ = 0.5), although the predominant prevalence peak and the corresponding vaccination peak around year 25 are both higher than their counterparts observed for the case where ζ = 0.
Furthermore, We found that employing a small fraction of committed vaccine recipients prevents a major epidemic by curbing disease prevalence. Such a finding holds for all ω (i.e., regardless of the population's responsiveness towards disease prevalence) as the magnitude of prevalence is too small to make a non-vaccinated individual switch to vaccinating ( Figure 13).

Discussion and Conclusions
We presented a series of SIR-network models with imitation dynamics, aiming to model scenarios where individuals commute between their residence and work, which is modelled by a commuting network where each node represents a suburb. These network models are able to capture diverse travelling patterns (i.e., reflecting local connectivity of suburbs), and different vaccinating behaviours affecting the global vaccination uptake and epidemic dynamics. We also analytically derived expressions for the basic reproduction number R 0 for the considered SIR-network models, and demonstrated how epidemics may evolve over time in these models.
We showed that the stable oscillations in the vaccinating dynamics are only likely to occur either when there is no population mobility across nodes, or only with limited commuting destinations. We observed that, compared to the case where vaccination is only provided to newborns, if vaccination is provided to the entire susceptible class, higher disease prevalence and more volatile oscillations in vaccination uptake are observed (particularly in populations which are relatively responsive to the changes in disease prevalence). A more expansive travelling pattern simulated in a larger network leads to the appearance of attractor dynamics in the relative proportion of vaccinated individuals, x, and the proportion of infected individuals, I, and the eventual convergence to the endemic, mixed equilibria, again if individuals are sufficiently responsive towards the changes in the disease prevalence. If individuals are insensitive to the prevalence, they are hardly affected by different vaccinating models and remain as unvaccinated individuals, although the existence of committed vaccine recipients noticeably delays the convergence to the non-vaccinating equilibrium. The presented models highlight the important role of committed vaccine recipients in actively reducing R 0 and disease prevalence, strongly contributing to eradicating an epidemic spread. Similar conclusions have been reached previously [34], and our results extend these to imitation dynamics in SIR-network models.
Previous studies drew an important conclusion that highly connected hubs play a key role in containing infections as they are more likely to get vaccinated due to the higher risk of infection in social networks [6,45]. Our results, verified by simulation experiments on a Greater Sydney commuting network, complement this finding by showing that a higher fraction of individuals who reside in highly connected suburbs (nodes) choose to vaccinate compared to those living in relatively less connected suburbs (nodes). These hubs, often recognized as business districts, also have significantly higher population influx as the destination for many commuters from other suburbs. Therefore, it is important for policy makers to leverage these job hubs in promoting vaccination campaigns and public health programs.
Overall, our results demonstrate that, in order to encourage vaccination behaviour and shorten the course of epidemic, policy makers need to carefully balance the following three considerations: ensuring a number of committed vaccine recipients exist in each suburb, utilizing the fact that people well-connected suburbs are more likely to vaccinate, and increasing individual awareness towards the prevalence change.
There are several avenues to extend this work further. This work assumes that the individuals from different suburbs (nodes) only differ in their travelling patterns, using the same epidemic and behaviour parameters for individuals from all nodes. Also, the same ω and κ are used for all nodes, by assuming that individuals living in all suburbs (nodes) are equally responsive towards disease prevalence and imitation. A greater level of accuracy in modelling can be achieved by establishing context-specific R 0 and imitating parameters for factors such as the local population density, community size [37], travelling rate [21], and suburbs' level of connectivity. For example, residents living in highly connected suburbs may be more alert to changes in disease prevalence, and adopt imitation behaviours more quickly. Bounded rationality can also be used to consider cases where individuals are not perfectly rational [48,49]. Different network topologies can also be used, particularly scale-free networks [41,42] where a small number of nodes have a large number of links each. These highly connected nodes are a better representation of suburbs with extremely high population influx (e.g., central business districts and job hubs). Dependency between epidemic dynamics and network topologies/properties are also of strong interest. Further investigation can be conducted by calibrating the network to real-world networks (other than the Greater Sydney commuting network that we studied here) to mimic a breakout in a targeted geographical region, such as a city, a suburb, or a statistical local area [46,47,50]. Previous studies also revealed that network topological metrics such as centrality measures, assortativity and robustness play an important role in the epidemic dynamics [51][52][53][54]. Further refinement on network-based modelling can also be attempted by introducing multiplex networks [55] where epidemic spreading and social interactions rely on two separate networks. It may also be instructive to translate the risk perception of vaccination and infection into tangible measures to demonstrate the aggregate social cost of an epidemic breakout, and help policy makers to visualize the cost effectiveness of different vaccinating strategies and estimate the financial burden for public health care. Parameters can also be calibrated to model other diseases with rich epidemic data, such as 2009 H1N1 influenza pandemic and 2003 SARS epidemic [43]. These considerations could advance this line of research to more accurately reflect contagion dynamics in urban environments, and provide further insights to public health planning.  Acknowledgments: This research was supported by the Sydney Informatics Hub at the University of Sydney, through the use of High Performance Computing (HPC) services.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Next Generation Operator Approach
The proposed model (15) can be seen as a finely categorized SIR deterministic model, so that each health compartment (Susceptible, Infected and Recovered) has M sub-classes. Let us denote I = {I 1 , I 2 , ..., I M } to represent M infected host compartments, and y to represent 2M other host compartments consisting of susceptible compartments y s = {y s1 , y s2 , ..., y sM } and recovered compartments y r = {y r1 , y r2 , ..., y rM }.
where i ∈ (1, M) and j ∈ (1, M). F i is the rate at which new infection enters infected compartments and V i is the transfer of individuals out of or into the infected compartments.
When close to disease-free equilibrium where S = S * = 1, the model can be linearized to: The next generation matrix, K, is then given by: where each entry K ij represents the expected number of secondary cases which an infected individual imposes on the rest of the compartments. F and V are given as follows: The basic reproduction number R 0 is given by the most dominant eigenvalue (or 'spectral radius' ρ) of K [22,31,32], and therefore:

Appendix B. Propositions and Proofs
Proposition A1. At the disease free equilibrium, x 0 has two solutions (x 0 = 0 or x 0 = 1) if x 0 and S 0 are uniform across all nodes.
Continuing with column 1, the column sum is: Similarly, the sums of each column are equal to 1, with all entries being non-negative population fractions. Hence, G is a Markov matrix.
Proposition A3. At the disease free equilibrium, x 0 has two solutions x 0 = x c = 0 or x 0 + x c = 1 if x 0 , x c and S 0 are uniform across all nodes.
Proof. In analogy with Proposition A1, with committed vaccine recipients, at the disease free equilibrium, S 0 i and x 0 i are uniform across all nodes, and therefore: Hence, under disease-free condition where I 0 = 0 andẋ i = 0, Equation (36) becomes: Equation (A15) leads to two solutions: x 0 + x c = 0 or x 0 + x c = 1.

Appendix C. Preliminary Analysis of the Oscillatory Behaviour of Epidemic and Vaccination Dynamics When Vaccinating Newborns
We found that a two-term exponential function ( f (x) = ae bx + ce dx ) is a satisfactory fit for Figure 7 as shown in Figure A1 and Table A1. Note that we do not fit a function for the case of ω = 1000 because of insufficient data points. Table A2 provides some preliminary insight on how ω affects the length of period. The length of period is relatively unchanged for epidemic and vaccination dynamics despite different oscillatory profiles. ω is positively correlated to the period, indicating that the time interval between prevalence peaks is longer if population is sufficiently responsive to prevalence change.  Table A1. Coefficients of fitting functions are summarised in Table A2.  We undertook an additional simulation experiment with a low transmission rate (β = 0.75), corresponding to a case where the disease is less infectious. Figure A2a,b show that the mixed, endemic equilibria at mid and high ω (ω = 2500 and ω = 3500) observed in Figure 6, in this case, are replaced by stable limit cycles, while the convergence to the pure, non-vaccination equilibrium at low ω (ω = 1000) is relatively unaffected. The reduced β also delays the first prominent peak with noticeably longer period between two peaks, as shown in Figure A2. Despite different oscillatory profiles, we note that higher ω generally corresponds to longer period, indicating that more responsive individuals extend the period between two epidemic peaks due to the higher vaccination coverage.