Modeling the effect of observational social learning on parental decision-making for childhood vaccination and diseases spread over household networks

In this paper, we introduce a novel model for parental decision-making about vaccinations against a childhood disease that spreads through a contact network. This model considers a bilayer network comprising two overlapping networks, which are either Erdős–Rényi (random) networks or Barabási–Albert networks. The model also employs a Bayesian aggregation rule for observational social learning on a social network. This new model encompasses other decision models, such as voting and DeGroot models, as special cases. Using our model, we demonstrate how certain levels of social learning about vaccination preferences can converge opinions, influencing vaccine uptake and ultimately disease spread. In addition, we explore how two different cultures of social learning affect the establishment of social norms of vaccination and the uptake of vaccines. In every scenario, the interplay between the dynamics of observational social learning and disease spread is influenced by the network’s topology, along with vaccine safety and availability.


Introduction
Human herding behavior or cascading is a convergence in opinion driven by social learning [1,2].Such herding could also be induced by payoff externalities, sanctions, preference interaction, direct communication, and/or observational influence [2].Observational influence results from integrating learned behavior or perception of others' opinions with one's own opinion.Observational social learning causes a subtle pressure on people to conform [3,4,5].That social learning takes place through different channels of information sharing, and through observation and/or perception of public choices [6].However, boundedly rational observational learning occurs when there is incomplete or insufficient information on the behavior of others, [7].
Social learning has been shown to play a vital role in people's decision-making even in the presence of information [8].Parents learn about the vaccination choices/attitudes of other parents and look for consensus vaccination signals [9].When parents share their feelings about the vaccination of their children while others observe these opinions, pressure is placed on other parents [10,11].Mixed messages and signals of parents' choices can cause problems for other parents when interpreting them.Incomplete information on the vaccine and the opinions of friends on social networks, creating boundedly social agents, can cause problems in vaccine acceptance and therefore in vaccine uptake and disease spread [12].
People and households are connected to each other through different types of networks.From a graph theory point of view, there are several types of network models, e.g.Erdős-Rényi (random) network model (ERN) and Barabási-Albert network model (BAN), see [13].Different network models represent various real-life systems.Many real-life networks are scale-free networks (SFN) in which the degree distribution is a power law with exponent 2 < γ ≤ 3, see [14].
Here, we consider bilayer networks in which two same-type overlapping networks and their mutual influence are of interest.The first network is the physical network through which faceto-face contacts take place and pediatric disease transmission occurs.The second network is a social (bidirectional) weighted network through which information and opinions about the vaccine are shared, shaping parents' decisions.Pediatric disease spreads on the physical network within and between households, while information and opinion sharing and perception transpire on the social network of households, parents in that case.We assume that the two networks overlap with parents linked to a larger number of other parents, relatives, coworkers, and distant friends, whose children might or might not be connected on the physical network.We introduce some adjustments to network models to reflect some important aspects of household networks.First, household networks should account for the number of children in households.Second, households without children must not have physical connections.That is, the expected degree of a household should increase with the number of occupying children, on average.
In this paper, we present a new model of parental decision-making to protect their children against a measles-like disease that spreads over household networks.This model considers boundedly rational observational social learning using a Bayesian aggregation formula.This Bayesian formula is different from that of Muller-Frank and Neri [7], who presented a quasi-Bayesian model of boundedly rational observational learning in a general context.We show that our model can give rise to social norms and at the same time encompasses other selection models like voting and DeGroot.Our new model considers socially bounded agents (parents who value children the most) who possess imperfect information about the vaccination choices of their network's neighbors.
We postulate that those agents can only perceive a correct or wrong message, probably, because of fear of retribution or confusion.There is a chance that the agents will send a correct message of their opinion with probability q and a message with probability 1−q, see [33].Using that model, we study the cascading opinion of vaccination in a boundedly rational observational social learning and compare it to other models of social pressure [29,34], see also [12].Moreover, we study the influence of those types of signal games on the spread of vaccine opinion and disease in social and physical networks in the presence of resource limitations represented in vaccine efficacy and accessibility, as well as vaccine safety.Finally, we study the effect of the presence of two cultures of social learning on the establishment of a social norm and so on vaccine uptake.

Networks.
To model the spread of the disease, we use an agent-based network model whose nodes are households N .Households are occupied by a number of children (C i , such that 0 ≤ C i ≤ n C , for i = 1, . . ., N ) who are connected through a physical network.The number n c is the maximum number of children per household.Parents are connected through a different (social, internet and physical) bi-directional weighted network via which they exchange opinions and information, and observe choices.Parents in a household without children could still be connected to other parents and shape their opinions.We use two types of networks: Erdős-Rényi (random) network model (ERN) and Barabási-Albert network model (BAN), see [13].In physical networks, we postulate the degree of nodes to be proportional to the number of children in the households.Meanwhile, the parent's network overlaps with the children's network via random rewiring of the children's network but with wiring probabilities that are greater than severance probabilities.We assume that the parent's network is weighted bidirectional with weights given by the (learning) probabilities q j,i , for i, j = 1, . . ., N and (i, j) is a social network link.See supplementary material (SI.Model) for complete information about the networks.An instance of a histogram for each random network is shown in Fig. S1.

Birth process.
A birth process is postulated to depend on the number of children occupying a household.The probability of a new pregnancy is modeled using a logistic function with a median of C * such that it decreases as the number of children in a household approaches n c .See supplementary material (SI.Model).A pregnancy lasts for a period of 280 days.The delivery of a newborn updates the number of children, but not the number of links in the children's network.We consider miscarriages and children's deaths to be rare and so are not included.

Disease spread.
We assume that a new measles-like disease (vaccine-preventable, pediatric) is spreading between children within and between households.We assume a mean length m p days of incubation period with a maximum of ℓ days.A new infection in a household occurs with a probability 1 where β is the probability of infecting a child in another household (through the physical network), β h is the probability of infecting a sibling within the same household, and I(i) is the total number of infected siblings in the same household.The number n I (i) is the number of infected children connected to the household i through the network.The number n I (i) is divided by C i to approximate the probability of infection formula based on the assumption that the children in a household have the same number of friends on average.We assume that the epidemic starts with I 0 infected children in different households and is randomly chosen.

Vaccination decision-making.
Parents are randomly allocated into three types: never-vaccinators, who oppose vaccination at all times, non-vaccinators, and vaccinators.We assume a small percentage of parents are nevervaccinators who are still sharing their opinions.Parents in the household i shape their subjective decision to vaccinate their children based on the reward of vaccination, given by π i = α i I − γ i A, where I and A are the total number of adverse events infected and vaccine, respectively, up to the decision-making.Parents who experience an adverse event due to vaccinating their children are switched to never-vaccinators in the sequel.The parameter α i is the degree of relevance of disease infectiousness and γ i is the degree of relevance of vaccine adverse events to the subjective opinion of the household i.The probability of accepting vaccination against the disease is given by pr i = 1/(1 + exp(−π i )) for the household i and is equal to zero for never-vaccinators.Probabilities pr i = 1/(1 + exp(−α i I 0 )) are used to generate the initial stance of parents in the household i toward vaccination.

Observational social learning.
Let the learning probability q i,j that the household i has the correct perception or learning about the household's j opinion or stance about vaccination.Meanwhile, 1−q i,j is the probability that the household i will make the wrong perception of the household opinion j'.Such social learning does not have to be symmetric, that is, q i,j ̸ = q j,i .For instance, followers of a celebrity learn from the celebrity more than in the other way around.Reciprocal or symmetric social learning involves q j,i = q i,j .Let household i have a set of vaccinator neighbors in the social network N V (i) with cardinality n V (i) and have a set of non-vaccinator's social network neighbors N N (i) with cardinality n N (i).Let also the total number of neighbors be n S (i) = n V (i)+n N (i), where N S (i) = N V (i)∪N N (i).Then parents in the household i make the decision to vaccinate their children based on the following posterior probability: or what we call the Bayesian aggregation rule in observational social learning.The rationale of the formula in equation (2.1) is that the prior probability of vaccination pr i is updated by independent information collected or perceived from neighbors in the network.A vaccinator is perceived to have that opinion with probability q i,j , and a non-vaccinator is perceived to have the opinion to vaccinate with probability 1 − q i,j .
Nondirectional social learning, in the sense of outward uniformity, means q j,i = q j,k =: q j for all i ̸ = k and for all j.
In case that q j = q for all j, then for 0 < q < 1.An uninformative probability q = .5results in no social influence on the parent, since P S (i) = pr i .
The model in equation (2.1) can be re-written as where π i = logit(pr i ) and It is also a model of bounded rationality.Putting P S(i) in the Boltzmann distribution form with temperature τ = .5,reveals a bounded rational model supported by physics [35].In that case, the reward for vaccination is given by is the induced social pressure to vaccinate.DeGroot model of selection [14], can also be seen in the term ∆ i if the neighbor j ∈ N S (i) is given a weight of logit(q j,i ) and using the discrete opinions to be valued as +1 for vaccination and −1 for no vaccination.This makes a stochastic DeGroot model of selection a special case of the Bayesian aggregation rule in equation (2.1).
Similarly, when q i = q and pr i = pr for all i as in the model of equation (2.2), the model can be re-written as where π = logit(pr) and δ i = n S (i) logit(q) is the degree of injunctive social norm practiced by household i. See [29,34,12] for the case δ i =: δ for all i, in which case which leads to the voting model of selection using G := 2 n V (i) n S (i) − 1 to decide the winner strategy when G > 0 and loser when G < 0. In other words, the latter model is also a special case of (2.1) in case of homogeneous observational learning.But since δ is nonnegative, the equation (2.5) can only reveal the human behavior when q ≥ 0.5 and so the model in (2.5) is less attractive than that in (2.1).
On each day, the parental position on vaccination is updated randomly based on the probabilities P S and vaccinators are chosen with the probability ρ of vaccinating all their children where ρ is the probability of getting access to vaccination based on the resources.

Epidemiological measures.
To analyze the effect of the probability of learning q j,i on the opinion of parents and on the spread of disease, we use a number of epidemiological measures: the size of epidemics, the peak of the epidemic, the uptake of vaccines, the number of vaccinators, and the basic reproduction number R 0 that we use only for the calibration of ERN.We use R 0 as an epidemiological measure and not as a threshold, see [36].The size of the epidemic is the total number of children infected at the end of the epidemic.The uptake of the vaccine is the total number of children who are vaccinated.The ending number of vaccinators is used to measure whether the opinion of vaccination becomes a consensus at the end of the epidemic.The basic reproduction number R 0 is defined to be the average number of secondary cases in a completely susceptible population.We use that definition to build an algorithm to estimate the value of R 0 (see Algorithm 1 in Supplementary Material SIII, Methods).In that algorithm, we use Bayes' theorem to calculate the probability that infection happens due to contact with the index case, which is then used to calculate the mean number of infections.To find the grand mean, we averaged over several simulations of disease transmission and then over the N households, which can include the index case, then finally over simulations of the various networks.

Model simulation.
The model is implemented using stochastic simulation for 100 runs to examine the effect of the probability of social learning q j,i on the uptake of vaccines and the spread of pediatric diseases.We assume that q j,i are uniformly distributed on q ± 0.05 for a prespecified value of q where 0 < q < 1.The stance towards vaccination and the disease states of infected children are updated at the beginning of each time step (day).Multiple infections can occur on the same day in the network and the numbers n I (i) and I(i) are updated every day for all i, i = 1, 2, . . ., N .An infected child on the j th day after infection jumps either to the end of the incubation period (recovers) or remains infected with a transition to the following day with probabilities given by the truncated exponential distribution (see supplementary material SI, Model).
Our simulation codes rely on the NumPy-compatible CuPy Python library [37], accelerated with NVIDIA CUDA [38] for parallel calculations in Graphical Processing Units (GPUs).Most of the calculations were performed on a CentOS workstation with 8 NVIDIA Tesla (Kepler) K80 GPU cards, each of which has 2496 CUDA cores and 12GB memory.For more information, see Supplementary Material SV, Codes of Simulations.
The eight GPU cards were used to distribute the 100 runs of stochastic simulations, with each card running one network simulation at a time.While the epidemic process can be run immediately after the network generations, we typically run the two processes separately.After network generation, the networks were saved to a hard drive, and at the beginning of the epidemic process, these networks were read from the hard drive.
• Average generation time of one network, including the savings of the networks to hard drive: -Erdő s-R'enyi Network (ERN): 21 seconds of wall clock time.
-Barab\'asi-Albert Network (BAN): 116 seconds of wall clock time.This network has more serial steps than the ERN, since the households are added one by one to the network.
• Average epidemic run time over one network: -Erdő s-R'enyi Network (ERN): 13 seconds of wall clock time.
-Barab\Asi-Albert Network (BAN): 11 seconds of wall clock time.The shorter time is due to the fact that the sizes of the epidemics are smaller than those on the ERNs.

Parameter values.
Model parameterization is done using the literature, calibration, and guesstimation.We use a number of N = 100, 000 households in a moderate-sized city with a random number of children in each household of a mean equal to two and half children.We assume that the mean degree in the children's network is 40 and the mean degree of the parents' network is 60 in the ERN.We assume that the epidemic starts with I 0 = 10 initially infected children that are randomly dispersed among N households.We postulate a disease of a mean incubation period of 11 days and a maximum of 16 days, [39, p.8].We find the values for β by calibration using the R 0 values between 12 and 18, [39, p.8].We assume that a fraction of 5% of the population will refuse to vaccinate at all (nevervaccinators), due to medical or ideological reasons.A full table of the definition of parameters and their values can be found in the supplementary material SII and Table S1.

Results
Using voting models of selection, the effect of the degree of the injunctive social norm or peer pressure δ on epidemic sizes and their peaks as well as the vaccine uptake is almost not noticeable for selected values of δ in [0.025, 0.225], Fig. 1.That is true for the case of Erdős-Rényi (random) network model (ERN) in Fig. 1 (a), (b) and (c), and clearer for Barabási-Albert network model in Fig. S1 (d), (e) and (f).The degree of injunctive social norm works differently on the two types of networks on vaccine uptake and on the sizes and peaks of epidemics as indicated by the results of the simulations shown in top panels versus bottom panels of Fig. 1.Irrelevant to the group pressure, vaccine uptake on Barabási-Albert networks is larger than on Erdős-Rényi networks.In case of Erdős-Rényi (random) network model (ERN) and using the general Bayesian aggregation rule in equation (2.1), we find richer dynamical behaviors than those when voter model with an injunctive social norm or peer pressure δ is used.First, let us assume that there are enough vaccines to vaccinate one child in every 100 children every day.When P adv = .0001,as q increases, the pressure imposed on parents increases and leads to higher uptake of vaccines and, consequently, to a smaller size and peak of epidemics, Figs. 2 (a), (b), and (c).But when the probability of an adverse event increases, and so more adverse events occur, as the probability of perceiving the correct position increases, vaccine uptake levels drop.When P adv = .001,the uptake of the vaccine will continue to increase with the values of q < .5, a pattern that changes for q > .5,Figs. 2 (d), (e), and (f).A higher probability of adverse events of P adv = .01,results in a decrease in vaccine uptake as the value of q increases beyond q ∼ .2.In that case, both the size and the peak of the epidemics increase with the increase in the probability q, Fig. 2   Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .01.
Those patterns change when we assume that vaccines are scarce, for instance, when it is only possible to vaccinate one child in every 1000 children every day.First, the uncertainty of the outcomes increases, compare panels of Fig. 2 to Fig. 3.The less vaccine available, the less number of adverse cases will appear.It takes a higher probability of an adverse event to effectively motivate parents to refuse vaccination of their children, Fig. 3 (b), (e) and (h).It is noticeable that the peak of the epidemics in that case do not change significantly with the change in the learning probability.Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .001.
In case of Barabási-Albert network model (BAN) and using the Bayesian updating rule in equation (2.1), the sizes of the epidemics are smaller than those in ERNs, compare panels of Fig. 2 to Fig. 4, and panels of Fig. 3 to Fig. 5.That is against the known fact that disease spreads faster on BANs.But that might be due to the increase in vaccine uptake that we observe in BANs in contrast to ERNs.The difference between BANs and ERNs is also true for the values of q.When the probability of adverse event increases to P adv = .01,vaccine uptake levels dropped from those levels when P adv = .0001and P adv = .001,compare panel (h) of Fig. 4 with panels (b) and (e).However, in that case, the learning probability causes a pressure to increase vaccine uptake even if it causes more adverse events; compare panel Fig. 4 (h) to panel Fig. 2 (h  Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .01.
Again, those patterns disappear when we assume that vaccines are scarce, and they are only available for one child in every 1000 children every day.Like in the ERNs, the outcome uncertainty increases, compare panels of Fig. 4 to Fig. 5.But in contrast to the ERNs, the learning probability does not have a significant influence over the vaccine uptake nor the size and peak of epidemics; compare panels of Fig. 3 to Fig. 5.  Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .001.
The parameter planes in Fig. 7 for the relevance of the disease α and the relevance of the adverse event of the vaccine γ to the rational choice component and for different values of q show consistent patterns with the simulations in Fig. 2. That means that for diseases and opinions spreading on ERN, the learning effect effectively suppresses the rational perception of the parent about the reward of vaccination.Parameter planes in Figs.S2 and S3 show consistent patterns as well.
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 q 4.4  6. Parameter planes of α against q and γ against q for epidemic sizes, the total number of vaccinated children and the peak of the epidemic on Erdős-Rényi (random) network model (ERN).In all simulations, the median value of the simulations is used to plot the parameter planes that are performed at P adv = .0001and ρ = .01.
Similarly, the parameter planes in Fig. 8 for the relevance of the disease α and the relevance of the vaccine's adverse event γ to the rational choice component and for different values of q show consistent patterns with the simulations in Fig. 4. Qualitatively, for diseases and opinion spreading on BAN, vaccine uptake and epidemics can have a symmetric reaction to high and low learning probability.Quantitatively, the parents' rational perception of the payoff of vaccination show influence over vaccine uptake and epidemics.As perceived disease risk increases, vaccine uptake increases, and epidemic sizes and peaks decrease, Fig. 8 (a), (b), and (c).And as the perceived risk of adverse events from vaccine increases, vaccine uptake decreases, and epidemic sizes and peaks increase, Fig. 8 (d), (e), and (f).The parameter planes in Figs.S4 and S5 also show consistent patterns. .Parameter planes of α against q and γ against q for epidemic sizes, total number of vaccinated children and the peak of the epidemic on Barabási-Albert network model (BAN).In all simulations, the median value of the simulations is used to plot the parameter planes that are performed in P adv = .0001and ρ = .01.
The probability of learning of parents can lead to a vaccination consensus in the case of Erdős-Rényi (random) network model (ERN) and Barabási-Albert network model (BAN) shown in Figs. 8 (a) and (b) as the number of vaccinators on the last day of the epidemic.However, a moderate to high learning probability is required in the case of ERN.The pattern stays exactly the same as in Fig. 8 (a) for ERN when P adv = .01,data are not shown to avoid redundancy.It changes drastically, however, for BAN when P adv = .01,Fig. 8 (c).In the latter, fewer parents end up accepting to vaccinate their children with a decline over the increase in the learning probability.Simulations are done on Erdős-Rényi (random) network model (ERN) using P adv = .0001in (a).It was done on Barabási-Albert network model (BAN) when using P adv = .0001in (b) and using P adv = .01 in (c).In all of the simulations ρ = .01.
Up to this point, the simulations of the model assumed that the population had a homogeneous culture.That is, parents show their actual preference or strategy with probabilities q i,j ∈ (q−.05, q+ .05)for one fixed value of q.In the next part, we examine the effect of the population having a cultural attribute.That attribute gives rise to a heterogeneous population, in which the population consists of two groups that use probabilities in two different regions: .1±.05 and .9±.05.We choose .1 and .9 as extreme case.We call the group of the latter type, subpopulation with attribute 1.The proportion of those with attribute 1 can affect the fate of vaccination uptake and the epidemic.
In the case of an epidemic spreading on an ERN, having a small proportion of parents with attribute 1, results in a large epidemic size and its peak, as well as less vaccine uptake.This is apparent in the simulations shown in Fig. 9.In the case of an epidemic spreading on a BAN, heterogeneity in the population's culture does not show a significant effect on the epidemics or vaccine uptake; see simulations in Fig. 10.If the proportion of parents with attribute 1 is more than one-half, vaccination becomes a consensus in the case of ERN; see Fig. 11 (a).Vaccination becomes a consensus in all cases for BAN; see Fig. 11 (b). .In all of the simulations are done using P adv = .0001and ρ = .01.

Discussions and conclusion
In this paper, we introduce a Bayesian aggregation model for boundedly rational observational social learning as a method of decision-making concerning vaccination.The learning model is based on social observations of neighbors on the social network through which information is spreading.Other models in the literature are shown to be special cases of our model.Some models, such as the degree of injunctive social norm δ, in contrast, lack the scale that our model enjoys.Using our new model, we explore the influence of social learning on the development of consensus on a network.The model also explores how heterogeneity in the culture of observational social learning affects vaccine uptake.We used stochastic simulations of disease and information processes that spread over two overlapping networks.The processes took part in one of two types of networks: Erdős-Rényi (random) network model (ERN) and Barabási-Albert network model (BAN).The results of observational social learning on those different networks and its mutual influence on the disease spread on the overlapping network were dependent on the network type and on the vaccine safety and availability.In ERN, as adverse events become rare due to their small likelihood or inaccessibility to vaccines, the increase in pressure ensued through declaring the adopted strategy with higher probability, especially when q ≥ .5, will result in an increase in vaccine uptake and eventually vaccination as a social norm.(Notice that q = .5 is the rational agent case.)The vaccine uptake in BAN is relatively higher than in ERN, and vice versa for epidemic sizes.Furthermore, the higher the probability of learning when q ≥ .5 in BAN, the less vaccine uptake will be.
Based on simulations of epidemic processes on the ERN and BAN, degree distribution plays a role in the levels of vaccine uptake and parental acceptance of vaccination.The uptake of the vaccine in the case of BAN is higher than in ERN.That might be for two reasons: As the disease is known to spread quickly in BAN, the accumulation of the number of cases will be rapid, increasing the subjective probability of vaccinating.The second reason is that the large number of neighbors will increase the inclination toward one of the two opinions even with small values of q.
Vaccine availability and accessibility interact with vaccine adverse events in a way that could influence parental opinion and thus vaccine uptake.In ERN, the more the vaccine is available and the less the chances of its adverse events, it will help to increase the uptake of the vaccine and establish the vaccination as a consensus.In BAN, that paradigm changes as the degree distribution has a heavy tail.
Mixed populations, with two different cultures of sharing and perceiving their opinion about vaccination, could have a significant effect on vaccine uptake.Having a fraction of the population with low learning probability q 1 in a population with q 2 > q 1 could decrease the resulting level of vaccine uptake and increase the size of the epidemic compared to a homogeneous population with the same learning probability q 2 .Social studies using surveys and behavioral game experiments might consider the personal characteristics of each parent that give rise to directional learning probabilities q i,j and q j,i .Population surveys can then be used to predict vaccine uptake levels.To elicit an increase in vaccine uptake, it is not enough to consider only the degree of parent link in the information network to spread awareness.Further efforts to promote social information exchange in social norm interventions will induce prosocial decisions about vaccination [3,4,5,40].According to our model, those efforts can lead to a consensus on the vaccination opinion and increase vaccine uptake even in the presence of a fraction of never-vaccinators and amid the vaccine's lack of safety and low availability.

Supplementary Material
A complete description of the model and methods, parameter values, and auxiliary figures is given in the Supplementary Online Material.

SI. Model
In this section, we give a complete description of physical and social networks and the dynamics of disease and information spread.

Children physical network
The first network N C connects children to households, so it is the network in which disease transmission occurs from child to child.Links between children represent face-to-face possibilities through school, clubs, communal activities, etc.These links are assumed to be static; that is, they do not change over time once they are established.That would be reasonable since the modeled outbreaks span a few weeks.At this level, we only allow transmission through the network, and no other transmission routes (such as environmental ones) are allowed.In that network, we assume that nodes are formed from a number of N households that are linked through physical transmission connections for each child (agent) in the household.The set of neighbors of the household i on this network is denoted by N P (i) and their size is n P (i).
1. Formation.A number of children C i (possibly zero) is randomly assigned to each household i, for i = 1, . . ., N , using a binomial distribution with parameters n c and p c .Here, we consider two types of networks: 1) Erdős-Rényi (random) network (ERN): Two households i and j are physically connected with the probability p • C i C j , which is proportional to the geometric mean of the number of children in both households.That probability could be interpreted through Newton's law of gravity and deserves future wider investigation.A household with zero children will not have physical connections.A regular ERN is a special case where C i = 1 for all i.The parameter p is selected to be less than 1/n c .
2) Barabási-Albert Network (BAN): In this network, we connect the first two households and then connect each of the following households, iteratively, according to the weighted degree distribution of the preceding households.That is, the household, say number j (j = 3, 4, . . ., N ), is connected to one of the households 1, 2, . . ., j − 1 according to their weighted degree distribution , k = 1, 2, . . ., j − 1, where d(k) is the degree of the household k.A household with zero children will not have physical connections.A regular BAN is a special case where C k = 1 for all k.See Fig. S1, for simulation instances of the overlapping ERNs and overlapping BANs. 2. Birth process.We assume that families decide to have a new baby (start a new pregnancy) with a probability dependent on the number of children per household.We assume that the probability of a new birth in a household i, if not already in pregnancy, is given by .
After 280 days C i is updated to C i + 1.The parameter σ is the population-level birth rate.
The parameter k measures the sensitivity of the probability to number of children in the household and C * is the median for that probability.We estimated k so as to have that probability approaching zero as C i approaches n C .Humans' gestation period is assumed to be 280 days.After delivery, parents either vaccinate or do not vaccinate their new babies.The first 280 days of the simulation of the model are discarded as a burn-in period.3. Disease transmission within and between households.Initially, a number of children I 0 are randomly and uniformly selected to be infected.A new infection happens in household, say number i, on any single day due to transmission within the household or between households with probability 1 where β is the probability to infect a child in another household (through the physical network), β h is the probability to infect a sibling within the same household, and I(i) = ℓ k=1 I(i, k) is the total number of infected siblings in the same household where I(i, k) is the number of infected siblings at day/stage k of the incubation period.The number n I (i) = j∈N P (i) I(j) is the number of infected children connected to that household i through the network.The number n I (i) is divided by C i in the probability of infection formula to reflect the approximating assumption that the children in a household have the same number of friends on average.Multiple infections can happen on the same day in the network and the numbers n I (i) and I(i) are updated every day for all i, i = 1, 2, . . ., N .4. Disease progression.The incubation period is assumed to be distributed according to a discretized exponential distribution with mean length of m p days and maximum ℓ days.An infected child on day j after infection moves either to the end of the incubation period (recovers) or stays infected into the following day with probabilities Q j and 1−Q j , respectively, such that ∞ j f (x)dx for j = 0, 1, . . ., ℓ − 1, and Q ℓ = 1, where f (x) is the probability density function of the exponential distribution.

Parental social network
The second network N P overlaps with the children's physical network and assumes larger social/information links between parents in households.Through social networks, parents share information and perceive opinions to and from network neighbors who are not necessarily physical neighbors but possibly friends, relatives, colleagues, online friends, etc.The set of neighbors of household i on this network is denoted by N S (i) and their size is n S (i).To form an overlapping social network with the children's network; links in the children's physical network are retained randomly with a probability of p re .Other links are added to any one of the updated non-connected households with probability p ad .

SII. Parameter Values
The model simulation is run for a large number of times for different sets of values of selected parameters.A parameterization of the model is given in Table 1.The parameters were assigned using the literature, calibration and guesstimation.The time unit is given in days.   .Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Erdős-Rényi (random) network model (ERN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .001and ρ = .01. Figure S3: Fig. S3.Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Erdős-Rényi (random) network model (ERN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .01and ρ = .01. .Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .001and ρ = .01. .Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .01and ρ = .01.

SIV. Estimating the Basic Reproduction Number R 0
The basic reproduction number is defined to be the average number of secondary infections caused by an index case introduced to a completely susceptible population till its recovery.To calculate the basic reproduction number as it is defined for the disease process over a network N C of households, we carry out algorithm 1.It is based on running the disease process for many number of repetition for an index case in a household that is selected uniformly from the network.Then the expected number of cases caused by the index case is estimated by the mean of binomial distributions whose probabilities are found using Bayes' theorem.
In particular, the probabilities P i (t) and P j (t) are the probabilities that the index case is the cause of infection of a new case within the same household i and in the neighbors j, respectively.Then the mean number of infections among IC i (t) and IC j (t) that are caused by the index case are given by IC i (t)P i (t) and IC j (t)P j (t), respectively.

SV. Coding of Simulation.
Our numerical simulations rely on the NumPy-comptatible Python library CuPy [37], accelerated with NVIDIA CUDA [38] for parallel calculations on Graphical Processing Units (GPUs).
Algorithm 1 To estimate the basic reproduction number using simulations of disease processes on networks.Input: number of simulation runs of the disease process L, number of simulated networks W , inputs for process in table S1 Output: R 0 begin for w = 1, 2, . . ., W ▷ Generate a network N C w of N households and number of children in households {C i (0) : i = 1, . . ., N }.Also, generate the starting stand of households about vaccination.
for i = 1, 2, . . ., N , such that C i (0) > 0 for ℓ = 1, 2, . . ., L ▷ Select one child in household i to be the index case ▷ Run the disease process till recovery of the index case, say till day T , which will change from one run to another.
▷ Record the number of new cases in household j, IC j (t), and the total number of prevalent children in household j less the index case, I j (t), for all j ∈ {i} ∪ N P (i).Record the total number of prevalent children in the outer neighbors k of the household j, I j,k (t), for t = 1, 2, . . ., T for all j ∈ N P (i), and k ∈ N P (j) − {i} ▷ Calculate: The probability that an infection in household i happened due to the index case in household i P i (t) = hβ hβ + (1 − (1 − hβ) I i (t) ) + j∈N P (i) (1 − (1 − β) I j (t)/C j (t) ) and the probability that an infection in household j happened due to the index case in household i P j (t) = β β + (1 − (1 − β) I i (t) ) + (1 − (1 − hβ) I j (t) ) + k∈N P (j)−{i} (1 − (1 − β) I j,k (t)/C j,k (t) ) where C j,k (t) are the number of children in the outer neighbors k of the household j.
▷ Calculate: R 0 (w, i, ℓ) = We introduce a sparse storage format for the adjacency matrices of the networks, where we record only the position of the nonzero elements in the upper triangular part by a single index number (see Fig. S6).With k denoting the index of the upper triangular entries in column-major order and (i, j) denoting the corresponding row and column indices the following relationships hold: The indexing formulas (5.1)-( 5.3) provide a fast and efficient storage and calculations through parallel implementations on the GPUs.

Figure 3 . 1 :
Figure 3.1: Fig. 1.Simulations of sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic for different values of δ in equation (2.5).Simulations are done on Erdős-Rényi (random) network model (ERN) in (a), (b), and (c), and on Barabási-Albert network model (BAN) in (d), (e) and (f).In all of the simulations P adv = .0001and ρ = .01.

Figure 3 . 2 :
Figure 3.2: Fig. 2. Simulations of sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Erdős-Rényi (random) network model (ERN) for different values of q.Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .01.

Figure 3 . 3 :
Figure 3.3: Fig. 3. Simulations of sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Erdős-Rényi (random) network model (ERN) for different values of q.Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .001.

Figure 3 . 4 :
Figure 3.4: Fig. 4. Simulations of sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN) for different values of q.Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .01.

Figure 3 . 5 :
Figure 3.5: Fig. 5. Simulations of sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN) for different values of q.Simulations are done using P adv = .0001in (a), (b), and (c), P adv = .001in (d), (e), and (f), and P adv = .01 in (g), (h) and (i).In all of the simulations ρ = .001.

Figure 3 . 6 :
Figure 3.6: Fig. 6.Parameter planes of α against q and γ against q for epidemic sizes, the total number of vaccinated children and the peak of the epidemic on Erdős-Rényi (random) network model (ERN).In all simulations, the median value of the simulations is used to plot the parameter planes that are performed at P adv = .0001and ρ = .01.

Figure 3 . 7 :
Figure 3.7: Fig.7.Parameter planes of α against q and γ against q for epidemic sizes, total number of vaccinated children and the peak of the epidemic on Barabási-Albert network model (BAN).In all simulations, the median value of the simulations is used to plot the parameter planes that are performed in P adv = .0001and ρ = .01.

Figure 3 . 8 :
Figure 3.8: Fig. 8. Simulations of the final total number of vaccinators on Erdős-Rényi (random) network model (ERN) and Barabási-Albert network model (BAN) for different values of q.Simulations are done on Erdős-Rényi (random) network model (ERN) using P adv = .0001in (a).It was done on Barabási-Albert network model (BAN) when using P adv = .0001in (b) and using P adv = .01 in (c).In all of the simulations ρ = .01.

Figure 3 . 9 :
Figure 3.9: Fig. 9. Simulations of sizes of epidemic, total number of vaccinated children, and the peak of the epidemic on Erdős-Rényi (random) network model (ERN) for different proportions of attribute 1. Simulations are done using ρ = .01 in (a), (b) and (c), and ρ = .001in (d), (e) and (f).In all of the simulations P adv = .0001.

Figure 3 . 10 :
Figure 3.10: Fig. 10.Simulations of sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN) for different proportions of attribute 1. Simulations are done using ρ = .01 in (a), (b) and (c), and ρ = .001in (d), (e) and (f).In all of the simulations P adv = .0001.

Figure S2 :
Figure S2: Fig. S2.Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Erdős-Rényi (random) network model (ERN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .001and ρ = .01.

Figure S4 :
FigureS4: Fig.S4.Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .001and ρ = .01.

Figure S5 :
Figure S5: Fig. S5.Parameter planes of α against q and γ against q for sizes of the epidemic, the total number of vaccinated children, and the peak of the epidemic on Barabási-Albert network model (BAN).In all of the simulations, the median value of the simulations are used to plot the parameter planes that are performed at P adv = .01and ρ = .01.
Fig. S6.One-index representation of network models.

Table S1 :
Model parameters, their description and base values.