Mathematical assessment of the impact of cohort vaccination on pneumococcal carriage and serotype replacement

ABSTRACT Although pneumococcal vaccines are quite effective in reducing disease burden, factors such as imperfect vaccine efficacy and serotype replacement present an important challenge against realizing direct and herd protection benefits of the vaccines. In this study, a novel mathematical model is designed and used to describe the dynamics of two Streptococcus pneumoniae (SP) serotypes, in response to the introduction of a cohort vaccination program which targets one of the two serotypes. The model is fitted to a pediatric SP carriage prevalence data from Atlanta, GA. The model, which is rigorously analysed to investigate the existence and asymptotic stability properties of the associated equilibria (in addition to exploring conditions for competitive exclusion), is simulated to assess the impact of vaccination under different levels of serotype-specific competition and illustrate the phenomenon of serotype replacement. The calibrated model is used to forecast the carriage prevalence in the pediatric cohort over 30 years.


Introduction
Nasopharyngial carriage of (or carriage with) Streptococcus pneumoniae (SP) is the primary transmission mechanism and precursor for pneumococcal disease (PD). SP, first isolated in 1881, is a gram-positive diplococcus with over 90 distinct capsular serotypes [32]. With the human nasopharynx as its main reservoir (mostly of young children), SP can spread from person-to-person [5]. Asymptomatic carriers of SP are responsible for most transmission [2,58]. The colonizing bacteria can also spread to the surrounding tissue, or invade the bloodstream of an individual, and can cause mild infection in the upper respiratory tract (acute otitis media and sinusitis) or more severe form of PD (such as pneumonia, meningitis and bacteremia) [56]. PD can be invasive (bacteremia, meningitis and bacteremic pneumonia) when SP infects a normally sterile body site [63] or non-invasive (non-bacteremic pneumonia, mastoiditis, conjuctivitis, acute otitis media or sinusitis). exclusive compartments of individuals who are susceptible (S), vaccinated (V), exposed (E), infectious (I), or recovered (R), range from adopting a deterministic to stochastic to individual-based formalism (see the review article by Lochen and Anderson [37], which also emphasizes that, even within the deterministic modelling paradigm, varying formalisms, such as SIS, SIR, SEI, SI and SIRS have been used to model pneumococcal carriage and disease [13-15, 33-35, 41, 54, 58, 62] (although the vast majority adopts an SIS structure)). This diversity of modelling paradigms and structures affirm the complexity of the natural history of pneumococcal carriage and scarcity of data on recovery from pneumococcal carriage and disease.
Lipsitch [34] presented an SIS model, with two SP serotypes and co-colonization, for analysing the impact of vaccination on competitive exclusion and coexistence of the serotypes. Their study shows that the presence of an NVT synergizes with the vaccine to reduce VT transmission. The model in [34] was used in a subsequent study to highlight the importance of serotype replacement and the need to further investigate the contrasting role of this phenomenon (where serotype replacement diminishes the impact of vaccines by allowing the NVT related PD to rise as the VT related PD is reduced by the serotype-specific vaccines) [35]. Melegaro et al. [41] developed and used an SIS model to predict the impact of PCV7 on the IPD incidence. In particular, this study showed that although elimination of the VT with PCV7 is achievable, serotype replacement with NVT poses a significant challenge or setback to attaining the optimal benefits of the vaccine. Bottomley et al. [8] presented an SIR model of SP carriage with multiple serotypes to predict the prevalence of VT and NVT following the PCV7 introduction. This model incorporates serotype-specific immunity, where a proportion of infected individuals is assumed to develop lifelong immunity to the serotype after surviving the infection. Sutton et al. [55] developed and analysed an age-structured SEIS model of PD with vaccination. Their model, which assumes a seasonal infection rate, was used to assess the impact of vaccination strategies at a community level.
Generalizing the model of Melegaro et al. [41], Choi et al. [14] proposed a model incorporating carriage with VT, NVT, and a VT-NVT co-colonization to predict the vaccine impact on VT and NVT IPD incidence. The structure of the model does not explicitly incorporate natural immunity. Similarly, the model in Choi et al. [15] splits the serotypes into groups and estimates the parameters governing infection transmission and competition effect between the groups. The effect of switching from PCV7 to PCV13, and that of no vaccination, on IPD cases is assessed. Cobey and Lipsitch [13] used their model to investigate how weak serotype-specific immune response can support coexistence of diverse serotypes.
Gaivao et al. [29] constructed a model to examine competition and facilitation among two co-colonizing serotypes. Global behaviour of the model is mathematically analysed and conditions for endemic persistence are explored. The model does not include vaccination or any other intervention. A two-strain epidemic model was constructed and analysed in [33] incorporating super-infection, where individuals infected with one strain can get infected with the other, replacing the first strain. It was shown, through rigorous global asymptotic stability analysis, that if regulated by super-infection, the disease dynamics can lead to serotype replacement even with a perfect vaccine. Co-colonization is not considered. Vaccination is assumed to protect completely against a target serotype (VT) and partially, fully, or not at all against a second serotype. The model also does not allow for natural immunity after carriage. It establishes, using an individual reproduction number for each serotype, that the presence of an NVT synergizes with the vaccine to reduce the VT transmission. Some modelling studies considered age structure in the dynamics of pneumococcal carriage, disease and vaccination. For instance, Snedecor et al. [54] used an age-structured model to evaluate the herd immunity effects of infant vaccination with PCV7 in the other age groups. Similarly, Wu et al. [65] developed an age-structured model to carry out the long-term cost-effectiveness analysis of PCV13 in Taiwan. Wasserman et al. [68] used an age-structured vaccination model, explicitly incorporating disease compartments, to assess the increase in IPD burden after a change in the dosing schedule in the United Kingdom.
The current study is based on the design, analysis and simulations of a novel dynamic model for gaining realistic insight into pneumococcal carriage and the transmission dynamics of two SP serotypes (one serotype is assumed to be VT and the other is assumed to be NVT) in the presence of a pneumococcal vaccine in a human population. The likelihood of co-colonization with more than two serotypes is insignificant (see [62] and references therein). Most (90%) of the samples detected with carriage in the study reported in Brugger et al. [7] contained two strains. Thus the model to be developed in this study will assume co-colonization with two SP serotypes. Furthermore, in the study reported in [7], half of the co-colonizations were found to be either VT only (20%) or NVT only (30%). Therefore, the proposed model assumes co-colonization with mixed VT-NVT types. That is, it is assumed that the vaccine provides positive efficacy against serotype i (i.e. i is VT) and no protection against serotype j (i.e. j is NVT). Similar co-colonization (between the mixed VT and NVT) is assumed in [41]. This categorization of SP serotypes into VT and NVT i and j respectively is suited to address serotype replacement arising from routine vaccination (the NVT serotype increases in prevalence to fill the niche vacated by the VT following the vaccination). The novel model to be developed also accounts for partial immunity after carriage, in addition to incorporating carriage with an individual serotype.
Being the main reservoir for transmission, the asymptomatic carriage states are the target of the conjugate vaccines [35,59] and are more frequent than disease states [29,35]. Since the contribution of symptomatic individuals to the transmission and overall mortality is negligible, they are not considered in this study. The model to be developed in this study tracks the natural history of carriage. In particular, individuals who had no prior history of carriage will be distinguished from those who have recovered from previous carriage with a certain serotype (or both serotypes), hence acquiring partial natural immunity against such serotype(s). This reduces the risk of subsequent carriage with the same serotype(s). The vaccine protection can wane over time.

Model formulation
The multi-serotype vaccination model for SP transmission dynamics is formulated by subdividing the total human population into mutually exclusive compartments based on the natural history of the disease (i.e. non-colonization, carriage with a single SP serotype i or j, co-colonization with both serotypes, recovery from carriage, etc.) and vaccination status (i.e. vaccinated or unvaccinated). Serotype-i is assumed to be VT, while j is NVT. Current carriage with SP serotype i(j) is represented using a superscript notation in state variables, vaccine protection is represented using the subscript notation vi, and natural immunity due to exposure history of serotype i and j is represented using the subscript notation ni and nj, respectively. Finally, a subscript v alone (and not followed by i) is used to distinguish the class of individuals who, although were vaccinated in the past, have experienced SP acquisition (consequently, the vaccine is assumed to provide no protection against that particular serotype). It is further assumed that co-colonization is acquired sequentially (and not simultaneously). Elbasha and Galvani [20] made a similar assumption for a model with multiple serotypes of human papillomavirus.
The total human population at time t, P(t), is subdivided into sub-populations of unvaccinated non-colonized (N u (t)), colonized with serotype k (C k u (t)) (where k ∈ {i, j}), co-colonized (C ij u (t)), those who recovered from carriage and having immunity against serotype k (R unk (t)), those who recovered from, and have natural immunity against, both serotypes i and j (R uninj (t)), and those with subsequent carriage with i(j) but having natural immunity against j(i) . The compartments representing corresponding vaccinated sub-populations are defined similarly and denoted with the subscript v (instead of u). The equations of the model are given by (A.1) in Appendix A. The basic qualitative features of the model are rigorously assessed in Appendix B. A schematic diagram of the model is depicted in Figure 1, and the state variables and parameters of the model are described in Tables 1 and 2, respectively.
The model (A.1) is an extension of numerous models for pneumococcal carriage that incorporate a vaccine, such as those in [8,14,15,33,34,41,54,55,58,62,65,68], by inter alia: (i) Explicitly tracking the natural history of carriage (the model in [29] considers carriage and co-colonization, but does not incorporate a vaccine) (ii) Allowing for the waning of vaccine-induced immunity over time (this was not accounted for in the aforementioned modelling studies, except in [15,33,41]) (iii) Allowing temporary partial natural immunity following carriage (this was not accounted for in the aforementioned modelling studies, except in [58,62,65]).   The model (A.1) is formulated based on the following key underlying assumptions: (a) The population is well-mixed, so that every member of the cohort can acquire infection from any other member of the cohort (i.e. we assume a homogeneously-mixed population).
(b) The average waiting time in each epidemiological compartment is exponentially distributed. (c) The vaccine does not provide therapeutic benefit against pneumococcal carriage. (d) Exposure to carriage with a certain serotype induces partial natural immunity against subsequent carriage acquisitions with the same serotype.
Numerical simulations of the model (A.1) will be carried out using the ranges and baseline parameter values tabulated in Table 2. The ranges and baseline values of many of the parameters of the model are obtained from the published literature, as described in Section 5.1.

Qualitative analysis of the model without vaccination
Before analysing the dynamics of the vaccination model (A.1), it is instructive to study the dynamics of the special case of the model in the absence of vaccination. The aim of this exercise is to determine whether adding vaccination to a base model for the dynamics of two serotypes of SP in a cohort alters its qualitative features. In other words, we seek to answer the question whether there will be dynamical features present in the vaccination model that may not be present in the reduced version of the model that does not include vaccination. This is explored below. It should be noted, first of all, that setting the vaccinerelated parameters (φ, ε and ω) to zero in the model (A.1) gives lim , R vninj (t)) = (0, 0, 0, 0, 0, 0, 0, 0).
The reduced version of the model (A.1), with these vaccination compartments set to zero, is henceforth called the reduced model.

Existence and asymptotic stability of equilibria
In this section, the existence and asymptotic stability of the various equilibria (carriagefree, boundary and coexistence) of the reduced model will be explored (theoretically and/or numerically). In this study, the coexistence equilibrium is defined as an equilibrium where some individuals carry serotype-i alone, some other individuals carry serotype-j alone, and some individuals are co-colonized with both serotypes.
It can be shown, using standard linearization around the CFE or the next generation operator method [17,18], that the associated reproduction number of the reduced model is given by where, R i 0 = cβ i γ i +μ and R j 0 = cβ j γ j +μ are the constituent basic reproduction numbers for serotype i and j, respectively (with each constituent reproduction number computed using the next generation operator method). R i 0 and R j 0 represent the average number of secondary carriage cases generated by a typical individual carrying the corresponding serotype(s) during their carriage period, when introduced into a completely non-colonized population in the absence of a pneumococcal vaccine [3,17,18,30]. The result below follows from Theorem 2 in [17]. The epidemiological implication of Theorem 3.1 is that, for the reduced model, a small influx of SP-colonized individuals into the population will not generate a large SP outbreak in the population if R 0 < 1. In other words, for a small initial number of SP-colonized individuals, the SP serotypes can be effectively controlled in the population if the basic reproduction number (R 0 ) can be brought to (and maintained at) a value less than unity.
To ensure that such effective control is independent of the initial number of SP-colonized individuals, a global asymptotic stability must be established for the CFE. This is explored below.
The proof of Theorem 3.2, based on using a comparison theorem [50,53], is given in Appendix C. The epidemiological implication of Theorem 3.2 is that, for the case of the reduced model, bringing (and maintaining) the associated reproduction number (R 0 ) to a value less than unity is necessary and sufficient for the effective control (or elimination) of both SP serotypes in the population. The result of Theorem 3.2 is numerically illustrated in Figure 2(a).
It is worth mentioning that heterogeneity exists in the values of R 0 based on geographical location. In particular, Gjini et al. [28] estimates R 0 for SP to be approximately 3.42 (95% CI: 2.96,3.96) in a daycare setting in Portugal. In a separate study, Gjini et al. [26] estimated a range of R 0 between 1.23 and 4.77 for SP in a similar setting across different countries. Lipsitch [34] uses a range of R 0 between 1.6 and 4 to analyse different carriage prevalence scenarios.

Boundary equilibria
The reduced model has two boundary equilibria (where only one of the two serotypes (the VT serotypy i or the NVT serotype j) exists at steady-state), as explored below. (i) Serotype-i-only boundary equilibria We explore conditions for the existence of the serotype-i-only boundary equilibria based on whether or not R i 0 is greater, less or equal to unity, as below.
All the rate parameters are in years.
Case 1: R i 0 > 1. For this case, the reduced model has a unique serotype-i-only boundary equilibrium, given by where

It follows from the expressions in (3) that the serotype-i-only boundary equilibrium (E 01 ) exists if
In this case, the reduced model has a serotype-i-only boundary equilibrium of the form where Since However, since η i is necessarily less than unity (see Section 2), it follows that this boundary equilibrium is not realistic for the reduced model. In other words, the reduced model does not have a serotype-i-only boundary equilibrium whenever R i 0 < 1. Case 3: R i 0 = 1. For this case, the reduced model has the following serotype-i-only boundary equilibrium: where It can be seen from (7) that . Thus, since η i is necessarily less than unity, it follows that the reduced model does not have a serotype-i-only boundary equilibrium when R i 0 = 1. In summary, it follows from the above analyses that the reduced model (noting that η i < 1) has a unique serotype-i-only boundary equilibrium (E i 01 ), which exists only if R i 0 > 1.
(ii) Serotype-j-only boundary equilibrium. Using the approach above, it can similarly be shown that the reduced model has the following unique serotype-j-only boundary equilibrium (which exists only when R j 0 > 1): where The theoretical results derived in this section are summarized below.

Theorem 3.3: The reduced model has a unique serotype-i-only (serotype-j-only) boundary equilibrium only if
, and no boundary equilibrium otherwise.

Asymptotic stability of boundary equilibria: special case
The asymptotic stability of the unique serotype-i-only (or serotype-j-only) boundary equilibrium of the reduced model will now be explored for a special case where the following simplifying assumptions hold: A1: Co-colonized individuals can only recover from both serotypes (i.e. the cocolonization ij) and not from any of the individual serotypes (i or j). That is, we set r 1 = r 2 = 0. Similar assumption was made in [20]. A2: No heterogeneity in recovery from single serotype carriage (i.e. γ i = γ j = γ ij = γ ). A3: No heterogeneity in the modification parameter for the reduction of carriage rate with serotype k due to prior carriage of serotype k (k = i, j). That is, we set η i = η j = 1.
The version of the reduced model with Assumptions A1-A3 is referred to as the simplified model. It can be seen that the reproduction number of the simplified model is given by Furthermore, the simplified model has a unique serotype-i-only (serotype-j-only) boundary equilibrium, denoted by E i 02 (E j 02 ), given, respectively, by where with . To investigate the local asymptotic stability of the two boundary equilibria of the simplified model, we compute the eigenvalues of the associated Jacobian matrix evaluated at each boundary equilibrium. In particular, the eigenvalues of the Jacobian matrix (of the simplified model) evaluated at the serotype-i-only boundary equilibrium (E i 02 ) are from which it follows that, for R i 0s > 1 (needed for the existence of the boundary equilibrium), the first eight eigenvalues (e i k ; k = 1, 2, . . . , 8) are real and negative. Further, e i 9 < 0 if and only if Similarly, it can be shown that the eigenvalues of the Jacobian matrix (of the simplified model) evaluated at the serotype-j-only boundary equilibrium (E j 02 ) are These results are summarized in Theorems 3.4 and 3.5.

Theorem 3.4: The serotype-i-only boundary equilibrium (E i
The theoretical results in Theorems 3.4 and 3.5 establish the conditions for the phenomenon of competitive exclusion in the reduced model. The serotype with the higher reproductive number (greater than unity) drives the other to extinction if the reproduction number of the second serotype is bounded by the threshold given in Theorems 3.4 or 3.5, whichever is applicable. These theoretical results are illustrated in Figure 2(b) and (c), showing initial solutions of the simplified model converging to the serotype-ionly or serotype-j-only boundary equilibrium (in line with the two theorems). It is worth mentioning that, if R i 0s = R j 0s ≥ 1 and θ i = θ j = 0, a stable boundary equilibrium that contains individual serotype-i and individual serotype-j (without co-colonization) exists. This scenario is depicted in Figure 2(e). Figure 3(b)-(e) show the phenomenon of competitive exclusion when one reproduction number is greater than unity and the other satisfies the inequality given in Theorems 3.4 or 3.5, whichever is applicable. Figure 3(a) shows the phenomenon of competitive exclusion in the absence of serotype-specific competition. For example, in the absence of competition presented by serotype i, setting θ i = 1 in the inequality in Theorem 3.4 yields the strict threshold R j 0s < 1, leading to the extinction of serotype j if R i 0s > 1. This corresponds to the much narrow stability region of serotype i-only equilibrium in Figure 3(a). The competitive exclusion of serotype i by j can be described similarly.
In this case, a boundary equilibrium containing serotype-i and serotype-j (without co-colonization) exists and is stable. The rate parameters are in years.

Coexistence equilibrium of the model without vaccination: numerical illustration
Extensive numerical simulations were carried out to determine whether or not the simplified model has at least one coexistence equilibrium (i.e. an equilibrium where both the VT and NVT and/or their co-colonization exist in the population). These simulations suggest the following conjecture.

Conjecture 3.1:
The simplified model has a locally-asymptotically stable coexistence equilibrium whenever the following conditions hold: carriage-free, boundary and coexistence equilibria) are locally-asymptotically stable for various levels of the serotype-specific competition parameters (θ i and θ j ). It should be recalled that setting θ i (θ j ) = 1 implies that serotype i(j) poses no competition to the (other) serotype j(i). This case (with θ i = θ j = 1 and R k 0s , k = i, j) is depicted in Figure 3(a). Similarly, the case with θ i (θ j ) < 1 implies that current carriage with serotype i(j) reduces the likelihood of acquiring carriage with the (other) serotype j(i). Furthermore, this figure shows that when the reproduction number of each of the two serotypes is less than unity, the carriage-free equilibrium becomes locally-asymptotically stable. Competitive exclusion occurs whenever the reproduction number of one serotype exceeds unity, while that of the other is less than the threshold given in Theorems 3.4 or 3.5 (in this case, the former serotype dominates and drives the latter to extinction).
The scenario where the competition parameter θ i is reduced from the value unity to 0.1 (i.e. carriage with serotype i reduces the likelihood of acquiring carriage with serotype j by 10%), while keeping θ j at unity, is depicted in Figure 3(b). This figure shows that while the stability region for serotype i (i.e. the stability region for the serotype-i-only boundary equilibrium) increases in size, the stability region of the coexistence equilibrium narrows. Similar dynamics were observed when θ j is decreased from unity to 0.15, while keeping θ i at unity (Figure 3 c). When both θ i and θ j are decreased, for instance to 0.1 and 0.15 respectively, the window for the stability region of the coexistence equilibrium significantly narrows (Figure 3 d). In other words, there is now significant competition (and less synergy) between the two serotypes. Finally, when θ i = θ j = 0, the coexistence equilibrium does not exist. Furthermore, in this setting, when min{R i 0s , R j 0s } ≥ 1, a stable equilibrium that contains serotype-i and serotype-j (without co-colonization) exists when (R i 0s , R j 0s ) is on the straight line R i 0s = R j 0s in the R i 0s -R j 0s plane (Figure 3 e). In summary, this study shows regions, in R i 0s and R j 0s parameter space, for the elimination, or coexistence, or existence of the two individual serotypes only. The size of these regions is greatly affected by the values of the competition parameters, θ i and θ j (for serotypes i and j, respectively). In particular, the size of the stability region of the coexistence equilibrium (when each of the reproduction numbers exceed unity) decreases with decreasing values of the competition parameters (below unity).

Qualitative analysis of the model with vaccination
The local asymptotic stability of the CFE of the vaccination model (A.1) will now be analysed. The CFE of this model is given by and N # v = φ ω+μ . It should be noted that since μ(ω + μ) > 0, it follows that N # u > 0 and N # v > 0. It should be observed that It can be shown, using the next generation operator method [17,18], that the reproduction number of the vaccination model, denoted by R v , is given by The result below follows from Theorem 2 in [17].
The epidemiological implication of Theorem 4.1 is that, for the cohort vaccination model (A.1), a small influx of SP-colonized individuals into the community will not generate a large outbreak of SP in the community if the vaccination program can bring (and maintain) the reproduction R v to a value less than unity. Using the parameter values in Table 2, the value of the basic reproduction number of the model (

Numerical simulations
The vaccination model (A.1) will now be simulated to assess the community-wide impact of implementing the cohort vaccination strategy (using the PCV13 vaccine), as well as to measure the impact of serotype-specific competition on the burden of the VT and NVT serotypes in the community. The model will be simulated using the parameter values and ranges in Table 2 (which are based on data for SP dynamics in the USA), unless otherwise stated. Further, for all the simulations the PCV13 vaccine will be introduced in the community when the population is already at an endemic steady state for SP. Before discussing the specifics of the two sets of simulations to be carried out, it is important to recall the concepts of cohort vaccination and serotype-specific competition, as briefly described below.
Cohort vaccination is generally defined as the immunization of the birth cohort. In the context of our study, cohort vaccination refers to the implementation of a pediatric vaccination program against SP. In the USA, anti-SP cohort vaccination consists of four PCV13 doses, administered to children at ages 2,4,6 and 12-15 months. In other words, the associated coverage parameter of the cohort vaccination model (A.1), denoted by φ, represents the proportion of children within this age cohort (henceforth referred to as the pediatric cohort) that received all four PCV13 doses each year.
The competition parameter θ i (θ j ) measures the assumed reduction in the likelihood of being colonized with serotype j(i) after already been colonized with serotype i(j). Competition may manifest itself via reduction in the rate at which new serotypes are acquired or through accelerated recovery of co-colonized individuals [5]. Hence, it is instructive to determine the impact of the competition parameters on the burden of the two serotypes in the community, particular as it relates to its important implications in serotype replacement. The modelling study by Bottomley et al. [8] predicts the persistence of low prevalence serotypes facilitated by the vaccine through reduction in serotype-specific competition (such serotypes would otherwise have been driven to extinction). Melegaro et al. [40] noted that serotypes which do not offer strong competition, but either persist longer or have high transmission potential, may replace other serotypes.

Justification of estimation of baseline values of model parameters
In this section, the estimate of the baseline values of the parameters in Table 2 will be justified, using available data and published literature. Specifically, numerous studies have generated SP carriage prevalence data for specific age groups, communities/regions or for specific serotypes. For example, Desai et al. [16] published carriage data for commonly carried SP serotypes in children in Atlanta, GA. Similarly, Balsells et al. [6] presented a systematic review and meta-analysis of SP carriage data among children under 5. Grant et al. [25] studied carriage prevalence in all ages in American Indians.
The recruitment rate parameter ( ) is estimated to be 3.7 × 10 6 per year. This is supported by CDC's provincial data for births (Vital Statistics Rapid Release) for 2018 [9], which showed that there were 3,788,235 births in the USA during 2018. The natural death parameter (μ) is estimated based on the average lifespan (1/μ) of the USA, which is approximately 78.6 years [10]. Thus μ = 1/78.6 per year. The carriage clearance rate parameter (γ ) is estimated based on the data reported by Bottomley et al. [8], Melegaro et al. [41] and Cobey et al. [13]. In particular, Bottomley et al. [8] estimated the per capita carriage clearance rate (with the 95% CI) for the low, medium and high transmission serotypes to be 0.049 (0.035,0.068), 0.039 (0.031,0.049), 0.031 (0.025,0.038) per day, respectively. Since the model in the Bottomley et al. study [8] does not consider co-colonization, the recovery rates are for individual serotypes (hence, the above estimates correspond to a range for γ i and γ j to be 11-18 per year).
Melegaro et al. [41] estimated the carriage duration (varies by host age) to be 72 days (for 0-1 year old children), 28 days (for 2-4 years of age), 18 days (for 5-17 years of age), 17 days (for those 18 years of age and older). Thus the Melegaro et al. study suggest a range of γ i and γ j to lie between 5 and 20 per year. The intrinsic durations of carriage are reported in Cobey et al. [13], quoting some studies, as lying between 25 and 225 days across serotypes. Hence, the Cobey et al. study suggests a range for γ i and γ j to lie between 1.6 and 14.6 per year. Considering all these estimates, we set γ i to be 9.9 per year and γ j , the clearance rate of the NVT j, to be 8.9 per year (both of which lie within the ranges described above). Furthermore, as noted by Gaivao et al. [29] (citing other studies), the clearance rate of cocolonization is lower than that of a single carriage (especially if co-colonization worsens the health of the host). Thus we estimate γ ij (also considering the range reported in the above mentioned studies) to be 7 per year. The parameters r 1 and r 2 , for the proportion of co-colonized individuals who recovered from carriage with serotype i or j, respectively, are estimated to be 0.4 and 0.3.
The estimate for the contact rate parameter (c) is adapted from Prem et al. [51], which reports the average number of contacts per individual (divided into 5-year age groups) per day to be 12.35, and the base-case probability of transmission of a single SP serotype per contact (β i , β j ) is assumed to be (0.0041,0.0034). A hazard ratio of 0.6 (0.39-0.9) (i.e. a colonized individual has an average of 60% chance of getting colonized again with the same serotype) was reported in Granat et al. [27]. Hence, we set the modification parameter for the reduction of carriage with serotype k (η k ; k = i, j) due to prior recovery from the same serotype to be 0.6 (0.39-0.9). The estimates for the efficacy of the PCV7 vaccine against the acquisition of carriage range from 0.5 in [8] to 0.756 in [41]. Rinta-Kokko et al. [52] reported an efficacy of 0.6 with further serotype-specific values, namely 0.76, 0.60 and 0.39, respectively, against acquisition of the vaccine types 6B, 9V and 23F, and 0.42 (95% CI 0.24−0.56) against all vaccine-type acquisition. Consequently, we set the vaccine efficacy against carriage acquisition (ε) to be 0.7. Sutton et al. [55] estimated the duration for loss of vaccine-induced protection to be 0.015 ( 0-0.03) per year. Similarly, Melegaro et al. [41] estimated the duration of vaccine protection to be 8.3 years (95%CI 5−20)., while Flasche et al. [23] estimated this duration (against carriage and disease) to be 5.8 (2.1-11.2) years. Consequently, following Melegaro et al. [41], we set the value of the vaccine waning rate (ω) to be ω = 1/8.3 per year (which also falls within the range of values in Flasche et al. [23]).
The cohort vaccination coverage parameter (φ) is varied across parameter space (for the various scenarios we simulated in this study). Further, since vaccinated colonized individuals are expected transmit carriage at a reduced rate, in comparison to unvaccinated colonized individuals, we set the associated modification parameter (α v ) to be 0.5. The modification parameter for the assumed increase in the recovery rate of vaccinated colonized individuals, in relation to unvaccinated colonized individuals (θ v ), is set at 2. Using data on unvaccinated children colonized with SP serotypes, Auranen et al. [5] showed a 10-fold reduction in the probability of acquisition of another SP serotype in the colonized children. Consequently, we set the competition parameter for VT (θ i ) to be 0.1 (i.e. children colonized with the VT are 90% less likely to acquire carriage with the NVT). Similarly, the competition parameter for the NVT (θ j ) is estimated to be 0.15 (within the 95% CI (0.05 -0.15) of the estimate for θ i ). It should be mentioned that, using data for 27 SP serotypes in Kenyan children aged 3-59 months, Lipsitch et al. [36] estimated the competition parameters to be as high as 0.52 (95% CI: 0.37 -0.63).

Effect of cohort vaccination coverage on dynamics of the SP serotypes
The cohort vaccination model (A.1) will now be simulated to assess the impact of the routine immunization against the vaccine serotype of SP in the cohort group. In particular, simulations will be carried out to assess the impact of variable levels of the cohort vaccination coverage on the dynamics of both the vaccine (VT) and non-vaccine (NVT) SP serotypes in the community. The simulations are aimed at exploring the impact of the cohort vaccination on curtailing the prevalence of the VT serotype (i), as well as in exploring how vaccination indirectly affects the prevalence of the NVT serotype (j). For these simulations, the vaccination model (A.1) is ran for a period of 100 years, using various values of the cohort vaccination coverage parameter (φ). It should be recalled that the cohort vaccination coverage parameter φ represents the fraction of the pediatric cohort in the community that is vaccinated. For these simulations, we consider the following three arbitrarily-chosen coverage levels of the cohort vaccination strategy implemented in the community: 1. Low coverage level of the cohort vaccination program: here, we set φ = 0.4 (i.e. 40% of children in the pediatric cohort are vaccinated each year). 2. High coverage level of the cohort vaccination program: here, we set φ = 0.8 (i.e. 80% of children in the pediatric cohort are vaccinated each year). It should be mentioned that this high level of cohort vaccination coverage is representative of the more realistic childhood vaccination in industrialized countries such as the USA, where the pediatric pneumococcal vaccine coverage, using PCV13, for 2017 in some US states (noticeably Arizona, Georgia, Minnesota, New Jersey and New York) was around 80% [12].
Simulations of the model for the worst-case scenario (i.e. with no vaccination, so that φ = 0) show that the two streptococcus serotypes (VT and NVT) remain at constant prevalence level at 10.3% and 0.93%, respectively (Figure 4 a), but with the prevalence of the VT exceeding that of the NVT. In this setting, the co-colonization remains constant at 0.036% (Figure 4 a). This is owing to the fact that, for these simulations, the constituent reproduction numbers of the cohort vaccination model, are ordered as When the cohort vaccination program is implemented at the aforementioned low (40%) coverage level, the simulation results obtained show prevalence in VT (NVT) decreases (increases) from the initial prevalence 10.3% (0.93%) to 1.1% (3.7%) (Figure 4 b). Further, the prevalence of the co-colonization decreases from initial prevalence of 0.036% to 0.015% (Figure 4 b). In this setting, initial prevalence of VT is higher than initial prevalence NVT, however, the prevalence of NVT exceeds the prevalence of VT eventually. This clearly illustrates the presence of the phenomenon of serotype replacement in the cohort vaccination model against the two serotypes of SP (A.1).
When the cohort vaccination coverage is increased to the high coverage level of the cohort vaccination program (i.e. φ = 0.8), the prevalence of the VT (NVT) decreases (increases) from 10.3% (0.93%) to 0.5% (4.04%) and prevalence of co-colonization decreases from 0.036% to 0.0072% (Figure 4 c). In summary, our simulations show that, while the use of cohort vaccination always reduces the prevalence of both the VT serotype (i), such strategy always increases the prevalence of the NVT serotype (j). In this setting, the phenomenon of serotype replacement is observed where the prevalence of NVT exceeds VT.

Effect of serotype-specific competition on the effectiveness of cohort vaccination
The vaccination model (A.1) is now simulated to assess the impact of the competition parameters, θ i and θ j , on the dynamics of the two serotypes in the community. This allows us to measure the impact of such competition on the effectiveness of the cohort vaccination program implemented in the community. The aim is to explore scenarios for serotype dominance and/or co-existence under this serotype competition and cohort vaccination setting. For these simulations, we consider the following three arbitrarily-chosen levels of the serotype-specific competition parameters (it should be mentioned that, in the absence of evidence of stronger or weaker competition offered by VT versus NVT, we choose to set θ i = θ j ): 1. Low serotype-specific competition: here, we set 0.7 ≤ θ i = θ j ≤ 1 and consider no (φ = 0), low (φ = 0.4) and high (φ = 0.8) coverage levels of the cohort vaccination program. 2. Moderate serotype-specific competition: this is defined by setting 0.3 < θ i = θ j < 0.7, combined with and no (φ = 0), low (φ = 0.4) and high (φ = 0.8) coverage levels of the cohort vaccination program. 3. High serotype-specific competition: in this case, we set 0 ≤ θ i = θ j ≤ 0.3, with no (φ = 0), low (φ = 0.4) and high (φ = 0.8) coverage levels of the cohort vaccination program.
Simulations for the case when the low serotype-specific competition scenario is combined with a low coverage of cohort vaccination, depicted in Figure 5(b), show that the prevalence of the VT, NVT and co-colonization decreases (from initial prevalence 15.1%, 11.68%, and 5.94%, respectively) to prevalence to 1.9%, 5.4% and 0.24%, respectively. When the low serotype-specific competition scenario is combined with a high coverage of cohort vaccination, the prevalence of the VT, NVT and co-colonization decrease (from initial prevalence 15.1%, 11.68% and 5.94%, respectively) to prevalence to 0.65%, 4.7% and 0.07%, respectively (Figure 5 c). Furthermore, in this setting, for low and high coverage of cohort vaccination, the prevalence of the NVT is higher than that of the VT (even though initial prevalence of VT is higher than that of NVT), showing a serotype replacement (Figure 5 b and c).
When the level of serotype-specific competition is moderate (i.e. we set θ i = θ j = 0.5) is combined with a low coverage of cohort vaccination, our simulations show that the prevalence of VT, NVT and co-colonization decreases (from initial prevalence 10.91%, 5.3%, 0.96%, respectively, to 1.4%, 4.4% and 0.09%, respectively , Figure 5 e). If high coverage level of the cohort vaccination is used, VT decreases to 0.58%, NVT slightly decreases to 4.3% and co-colonization decreases to 0.036% (Figure 5 f). In this setting, for low and high coverage of cohort vaccination, the prevalence of NVT is higher than that of the VT (showing a low or moderate level of serotype replacement).
It is worth noting from the above simulations for the low and moderate levels of the serotype-specific competition scenarios that the two serotypes (VT and NVT) always coexist. Further, low prevalence levels of co-colonization were observed when no or low coverage level of the cohort vaccination is used (as demonstrated in Figures 5 a and b). It Other parameter values used are as given in Table 2. The reproduction numbers for these simulations are given as follows: (i) no cohort vaccination (φ = 0): The rate parameters are in years.
is intuitive to expect the prevalence of co-colonization when low and moderate serotypespecific competition scenarios are used. This is because these scenarios are associated with higher values of the competition parameters (θ i and θ j ), which, in turn, increase the likelihood of carriage with the other serotype in individuals that were already colonized. In other words, low and moderate levels of serotype-specific competition favours coexistence, rather than serotype removal.
Although a few studies, such as those by Lipsitch et al. [36] (using data from Kenya) and Bottomley et al. [8] (using data from The Gambia), gave high estimates for the competition parameters (0.52 and 0.63 (0.407, 0.975), respectively), majority of the published studies on the dynamics of SP serotypes (particularly those based on using data from European countries) estimated very low values for the competition parameters. For example, using data from a longitudinal study of children in three Danish day-care centers, Auranen et al. [5] estimated the competition parameters to lie in the range 0.05-0.15. Similarly, the study by Gjini et al. [28], based on using data from Portugal, estimated the competition parameters to be in the range 0.04-0.09 (clearly showing a very strong serotype-specific competition in their data set). Numminen et al. [49] used data from Norway to estimate the competition parameters to lie in a range of 0.06-0.14. Finally, Gjini [24] estimated (based on the posterior distributions) ranges of 0.06-0.11 for Portugal, 0.03-0.05 for Norway, and 0.04-0.05 for Greece).
We also simulated the vaccination model (A.1) for the case when the serotype-specific competition is high (i.e. we simulated the model (A.1) with θ i = θ j = 0.2). For low coverage of cohort vaccination, the results obtained show that the prevalence of the VT (NVT) decreases (increases) from initial prevalence 10.12% (1.3%) to 1.14% (3.8%) (Figure 5 h). For this setting, the prevalence of the co-colonization decreases from initial prevalence 0.086% to 0.025%. For high coverage level of the cohort vaccination, VT (NVT) decreases (increases) to 0.5% (4.1%) and the co-colonization decreases to 0.012% (Figure 5 i). In this case, the prevalence of NVT much higher than that of the VT (clearly demonstrating serotype replacement).
In summary, the simulations in Figure 5 show the following main results: (i) When there is no cohort vaccination, the prevalence of co-colonization decreases as the level of the serotype-specific competition increases. Moreover, the prevalence of co-colonization decreases as level of cohort vaccination coverage increases. (ii) The values of the constituent reproduction numbers (R i v and R j v ) of the vaccination model (A.1) do not change with changing values of the competition parameters (θ i and θ j ; since these parameters do not appear in the expressions of the constituent reproduction numbers). This is because the two parameters are associated with cocolonization, which is secondary (and not primary) infection. Hence, they do not appear in any of the next generation matrices (used to compute the reproduction number of the vaccination model (A.1)). (iii) For each level of the cohort vaccination program (low or high), the total prevalence (the sum of the prevalence of VT, NVT and co-colonization) remains the same (at all time t) for all levels of the serotype-specific competition. (iv) The prospect of vaccination-induced serotype replacement (where the NVT replaces the VT, but without causing it to go extinct) is significantly increased with increasing level of both serotype-specific competition and cohort vaccination coverage.

Discussion and conclusion
Pnuemococcal disease (PD), a major respiratory disease caused by nasopharyngial carriage of any of the over 90 serotypes of the bacterium Streptococcus pneumoniae (SP), continues to cause major public health burden globally. The populations at highest risk of PD, namely children under the age of 5 and adults of age 65 or older, continue to suffer high mortality and morbidity annually. Vaccination, using any of the currently available pneumococcal conjugate (PCV) and polysaccharide (PPV) vaccines, is the most effective intervention against PD, and many countries around the world have included pneumococcal vaccine in their regular age-based immunization programs. Each vaccine is serotype specific. That is, each vaccine targets specific SP serotypes, with varying efficacy across the targeted serotypes. The non-vaccine SP serotypes (NVT) have emerged (with increasing prevalence) due to the lack of resource competition with the vaccine serotype (VT), owing to the fact that the incidence and prevalence of the latter have greatly declined following the introduction of the aforementioned vaccine into human populations. This dynamic phenomenon, of the increase in the prevalence of NVT (following the dramatic decline of the VT due to mass immunization), is called serotype replacement. This study presented a novel mathematical model for gaining insight into the population-level impact of pneumococcal vaccination on the transmission dynamics of SP carriage in a cohort. The resulting dynamic model, which monitors the temporal dynamics of two SP serotypes (VT and NVT), takes the form of a deterministic system of nonlinear differential equations. The model allowed for the assessment of the carriage history of the two SP serotypes in the presence of a vaccine. The SP carriage history was needed to enable the proper accounting of the development of natural immunity against SP serotypes due to past exposure. This immunity impacts subsequent colonizations after recovery. The model allowed for co-colonization with the two serotypes, which was assumed to occur only via the sequential carriage of individuals already colonized with the other serotype (i.e. the model does not allow for simultaneous carriage with the two SP serotypes).
The special case of the model, in the absence of vaccination (termed the reduced model), was rigorously analysed to gain insight into its dynamical features. Our analyses revealed that the carriage-free equilibrium of the reduced model is globally-asymptotically stable whenever the associated basic reproduction number of the model (denoted byR 0 ) is less than unity. The epidemiological implication of this result is that the two SP serotypes can be eliminated from the cohort if the reproduction threshold (R 0 ) can be brought to (and maintained at) a value less than unity. In other words, this study shows that bringing (and maintaining) the threshold quantity R 0 to a value less than unity is necessary and sufficient for the elimination of the two serotypes in the cohort. To the authors' knowledge, this may be the very first study to establish necessary and sufficient conditions for the theoretical elimination of the two SP serotypes in a cohort, by proving the global asymptotic stability property of the associated carriage-free equilibrium of the model. Some modelling studies (e.g. [28]) have established local stability of the CFE, and, in the case of two serotypes, local stability of the boundary-and coexistence equilibria. It is worth mentioning, however, that the model in [28] also considers the possibility of both serotypes being VT or both NVT.
Conditions for the existence of the boundary equilibria (corresponding to the VT-only or NVT-only) and a coexistence equilibrium (where both the VT and NVT, and/or the co-colonization co-exist) of the reduced model (without vaccination) were theoretically derived. The VT-only and NVT-only boundary equilibria were shown to be locallyasymptotically stable under certain conditions. The implication of this result is that, for the version of the model without vaccination, the epidemiological phenomenon of competitive exclusion can occur if one serotype has a reproduction number greater than unity, and the other less than a certain threshold (in this case, the former drives the latter serotype to extinction). Furthermore, when the two reproduction numbers of the reduced model are both greater than unity, the two serotypes (and their co-colonization) may co-exist. In this case, the serotype with the higher reproduction number dominates the other serotype (but does not drive it to extinction). This gives rise to the phenomenon of serotype replacement.
Our theoretical analysis, coupled with extensive numerical simulations, further showed that the level of serotype-specific competition (measured by the parameters θ i and θ j ) plays an important role in the dynamics of the SP serotypes and the co-colonization in the cohort. Specifically, we showed that the size of the stability regions for the coexistence equilibrium of the reduced model significantly decreases with increasing level of serotype-specific competition. In fact, the stability region disappears for high serotype-specific competition. For such high competition, the reduced model no longer has a coexistence equilibrium (only the carriage-free and the two boundary equilibria exist in this case). Thus this study shows that the mechanism of serotype-specific competition has the dual role of facilitating or impeding serotype coexistence, depending on the strength or level of the competition (lower level of competition supports coexistence, while higher competition level impedes coexistence).
The full model, which incorporates the routine pediatric/cohort vaccination (which entails the administration of four PCV13 doses, given to children in the pediatric cohort), was also rigorously analysed. In particular, it was shown that the carriage-free equilibrium of the full model was also locally-asymptotically stable whenever the associated reproduction number of the full model (denoted by R v ) is less than unity.
The model developed in this study can be extended in numerous ways. For instance, age structure and pneumococcal disease manifestations, which are essential for analysing the effectiveness of age-based immunization programs and formulating effective policy recommendations to reduce disease burden resulting from SP carriage, can be explicitly incorporated into the model. The model developed in this study is based on assuming a well-mixed population. This simplifying assumption (needed for mathematical tractability) constitutes a limitation of the model, and the results we presented could very well be affected if this assumption is relaxed. Our study can be extended by relaxing the homogenous-mixing assumption, allowing, instead, a heterogenous-mixing (i.e. preferential mixing) framework in a population that is also stratified by age. It is important to optimize the trade-off between vaccine coverage of the number of pathogenic serotypes (hence, maximizing the beneficial effects) and the risk of increase in the pneumococcal disease (PD) burden following the replacement of the VT with the NVT [35]. Hence, there is a need to explicitly include the dynamics of PD in the model (because the impact of vaccination on PD will vary with varying serotype distributions post-vaccination [35]). The dosing schedules of the pneumococcal vaccination program can also be explicitly incorporated in the model, in addition to accounting for the related costs and vaccine coverage/compliance associated with the vaccine dosing schedules (some pneumococcal vaccine modelling studies, such as those in [14,15,19,55,68], have accounted for dosing schedule and its overall impact on the dynamics of the disease). The suggested extensions can further contribute to the efforts to realistically assess the community-wide effectiveness of pneumococcal vaccination programs.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
One of the authors (ABG) acknowledge the support, in part, of the Simons Foundation (Award #585022) and the National Science Foundation (Award #1917512). The authors are grateful to the anonymous reviewers for their constructive comments.

Appendix. Model equations
Recall that, in the formulation described in Section 2, the total human population at time t, denoted by P(t), is subdivided into the subpopulations of unvaccinated non-colonized (N u (t)), unvaccinated colonized with serotype i (C i u (t)), unvaccinated colonized with serotype j (C j u (t)), unvaccinated cocolonized with serotypes i and j (C ij u (t)), unvaccinated who recovered from carriage and having immunity against serotype i (R uni (t)), unvaccinated recovered from carriage and having immunity against serotype j (R unj (t)) and unvaccinated who recovered from, and have natural immunity against, serotypes i and j (R uninj (t)), those vaccinated non-colonized (N v (t)), vaccinated colonized with serotype i (C i v (t)), vaccinated colonized with serotype j (C j v (t)), vaccinated co-colonized with serotypes i and j (C ij v (t)), recovered from (and having natural immunity against) serotype i (R vni (t)), recovered from serotype j and having vaccine-derived immunity against serotype i (R vnj (t)), recovered from (and having immunity against) serotypes i and j (R vninj (t)), vaccinated and colonized with serotype i with natural immunity against serotype j (C i vnj ) and vaccinated colonized with serotype j having immunity against serotype i (C j vni ). Individuals acquire carriage with serotype i, j or the co-colonization at the rates λ i and λ j , respectively given by where c is the effective contact rate, β k (with k = i, j) is the probability of carriage with serotype i or j, respectively, per contact. Similarly, the parameter 0 < α v < 1 accounts for the assumed reduced infectiousness of vaccinated colonized individuals, in comparison to the corresponding unvaccinated colonized individuals. The model for the dynamics of two SP serotypes, in the presence of a routine cohort vaccination program, is given by the following deterministic system of nonlinear differential equations (where a prime represents differentiation with respect to time t) The flow diagram of the model (A.1) is depicted in Figure 1, and its associated state variables and parameters are described in Tables 1 and 2, respectively.
The equations for the infected components of the reduced model can be written in the form: where the next generation matrices, F and V, and the non-negative matrix, S, are given, respectively, by with N * u = μ . Since N u (t) ≤ N * u = μ , in u , it follows that matrix S is non-negative. Hence, Equation (C.3) can reduced to the following linear inequality: Since all the eigenvalues of F−V have negative real part for R 0 < 1 (from the local asymptotic stability result in Theorem 3.1), it follows that the linearized system of differential inequality (C.4) is asymptotically stable if R 0 < 1. That is,  (NVT) generated from a cohort of 6-59 months children in Atlanta, GA, in 6 study periods, 6 months apart, from July 2010 to June 2013 (the end of the first study period is considered as the initial time for the start of the simulations), as reported in [16].
is based on the serotype distribution given in [16] for the first study period). Based on this, the initial conditions for each cohort vaccination coverage level is given below: Simulations of the vaccination model (A.1) for the worst-case scenario (i.e. with no vaccination, so that φ = 0) show that the prevalence of the streptococcus serotype 19A (35B) increases (decreases) from the initial prevalence 7.3% (2.5%) and reaches 37.93% (0%) at steady state, but the co-colonization remains constant at 0% for the entire 30 years period (Figure A2 a). For these simulations the constituent reproduction numbers of the cohort vaccination model are ordered as . It is worth mentioning that even if both R i v and R j v are greater than unity, the solutions converge to serotype 19A-only boundary equilibrium. This suggests that the competitive exclusion results given in Theorems 3.4 and 3.5, for the reduced model discussed in Section 3, also extend to the full vaccination model (A.1).
When the cohort vaccination program is implemented at the aforementioned low (40%) coverage level, the simulation results obtained show that prevalence of serotype 19A (35B) increases (slightly decreases) from the initial prevalence 7.3% (2.5%) and reaches 16.06% (0%), but the co-colonization remains constant at 0% for the entire 30 years period (Figure A2 b).
For high coverage level of the cohort vaccination program (i.e. φ = 0.8), there is a decrease (increase) in the prevalence of the 19A (35B) from the initial prevalence 7.3% (2.5%) which reaches 0.89% (4.21%) at steady state, and the co-colonization increases to 0.013% (Figure A2 c). In this In summary, our simulations show that, while the use of cohort vaccination reduces the prevalence of 19A, high cohort vaccination coverage increases the prevalence of 35B. For high level of cohort vaccination coverage, the presence of the phenomenon of serotype replacement in the cohort vaccination model against the two serotypes of SP (A.1) was illustrated. Furthermore, when there is no, or there is low, coverage of cohort vaccination, the phenomenon of competitive exclusion is observed, where 19A (the serotype with higher reproduction number) drives serotype 35B to extinction (Figures A2(a) and (b)). However, it is worth highlighting that the competitive exclusion was not observed in the simulations carried out in the presence of cohort vaccination with high coverage (Figure A2 c).