Use of an individual-based model of pneumococcal carriage for planning a randomized trial of a vaccine

For encapsulated bacteria such as Streptococcus pneumoniae, asymptomatic carriage is more common and longer in duration than disease, and hence is often a more convenient endpoint for clinical trials of vaccines against these bacteria. However, using a carriage endpoint entails specific challenges. Carriage is almost always measured as prevalence, whereas the vaccine may act by reducing incidence or duration. Thus, to determine sample size requirements, its impact on prevalence must first be estimated. The relationship between incidence and prevalence (or duration and prevalence) is convex, saturating at 100% prevalence. For this reason, the proportional effect of a vaccine on prevalence is typically less than its proportional effect on incidence or duration. This relationship is further complicated in the presence of multiple pathogen strains. In addition, host immunity to carriage accumulates rapidly with frequent exposures in early years of life, creating potentially complex interactions with the vaccine’s effect. We conducted a simulation study to predict the impact of an inactivated whole cell pneumococcal vaccine—believed to reduce carriage duration—on carriage prevalence in different age groups and trial settings. We used an individual-based model of pneumococcal carriage that incorporates relevant immunological processes, both vaccine-induced and naturally acquired. Our simulations showed that for a wide range of vaccine efficacies, sampling time and age at vaccination are important determinants of sample size. There is a window of favorable sampling times during which the required sample size is relatively low, and this window is prolonged with a younger age at vaccination, and in a trial setting with lower transmission intensity. These results illustrate the ability of simulation studies to inform the planning of vaccine trials with carriage endpoints, and the methods we present here can be applied to trials evaluating other pneumococcal vaccine candidates or comparing alternative dosing schedules for the existing conjugate vaccines. Author Summary Streptococcus pneumoniae, a bacterium carried in the nasopharynx of many healthy people, is also a leading cause of bacterial pneumonia, sepsis, and ear infections in children aged five years and younger. Vaccines targeting select strains of S. pneumoniae have been effective, and the development of new vaccines, particularly those that target all strains, can further lower disease burden. For clinical trials of these vaccines, the number of study participants needed depends on the expected effect of the vaccine on a conveniently measured outcome: asymptomatic carriage. The most economical way to test a vaccine for its effect on carriage is by measuring prevalence at a specific time, and comparing vaccinated to unvaccinated participants. The relationship between incidence (or duration) and prevalence is complex, and changes with time as children develop natural immunity. We explored this relationship using a mathematical model. Given a vaccine efficacy, our computer simulations predict that fewer study participants are needed if they are vaccinated at a younger age, taken from a population with intermediate levels of transmission, and sampled for carriage at a certain time window: 9 to 18 months after vaccination. Our study illustrates how simulation studies can help plan more efficient vaccine trials.

rapidly with frequent exposures in early years of life, creating potentially complex interactions 23 with the vaccine's effect. We conducted a simulation study to predict the impact of an 24 inactivated whole cell pneumococcal vaccine-believed to reduce carriage duration-on 25 carriage prevalence in different age groups and trial settings. We used an individual-based 26 model of pneumococcal carriage that incorporates relevant immunological processes, both 27 vaccine-induced and naturally acquired. Our simulations showed that for a wide range of 28 vaccine efficacies, sampling time and age at vaccination are important determinants of 29 sample size. There is a window of favorable sampling times during which the required sample 30 size is relatively low, and this window is prolonged with a younger age at vaccination, and in a 31 trial setting with lower transmission intensity. These results illustrate the ability of simulation 32 studies to inform the planning of vaccine trials with carriage endpoints, and the methods we 3 present here can be applied to trials evaluating other pneumococcal vaccine candidates or 34 comparing alternative dosing schedules for the existing conjugate vaccines. 35

54
For encapsulated bacteria such as Streptococcus pneumoniae [1], Haemophilus 55 influenzae [2], and Neisseria meningitidis [3], asymptomatic carriage in the human upper 56 respiratory tract is a precursor to mucosal or invasive disease. The population of bacteria in 57 9 two settings that differed in their transmission intensity: the higher transmission setting had an 144 under-five carriage prevalence of 66%; the lower transmission setting, 55%. 145 For the higher transmission setting, we ran 50 simulations of the vaccine trial using 146 different random seeds and recorded the carriage prevalence every month (defined as 30 147 days), starting from birth to 24 months after vaccination (Fig 1). For both infants and toddlers, 148 all vaccine efficacies led to reductions in prevalence throughout the follow-up period. Higher 149 efficacies resulted in greater reductions in carriage. However, that marginal benefit attenuated 150 with time as both controls and vaccinees acquired more natural immunity from carriage 151 episodes. Similar patterns were observed in the toddler trials, but with smaller reductions in 152 prevalence (Fig 2A-C). For the infants, the prevalence in the control and vaccine arms followed non-monotonic 176 trajectories over the course of the follow-up period. In the infants, the median prevalence in 177 the control arms started at 74% at 2 months of age, peaked at 91% at 8 months of age, and 178 then declined (Fig 2A-C, Fig 1A). The timing of the peak is consistent with previously 179 reported data from Kilifi,Kenya [31]. In the vaccinated infants, the median prevalence peaked 180 at the same time, at 8 months of age for the 3 c.e. vaccine efficacy, or slightly earlier, at 5 181 months of age for the 5 c.e. and 10 c.e. wSP vaccine efficacies (Fig 2A-C, blue). For the 182 toddlers, who are vaccinated later in life at 12 months of age, the age-specific prevalence in 183 both the control and vaccine arms steadily declined across the 24-month follow-up period ( Fig  184

2A-C, purple). 185
From the joint trajectory of the control and vaccine arm prevalence over the follow-up 186 period, we determined how the sample size required for a two-sample test of equal proportion 187 varied with sampling time. We assumed a 5% type I error probability, 80% power, and 188 balanced arms, and use the term "sample size" to refer to the combined size of both arms. In 189 infants, for all vaccine efficacies, the median sample size decreased dramatically-almost 190 ten-fold or more-in the period 3 to 9 months post-vaccination, plateaued, and then started 191 increasing around 18 months post-vaccination. In toddlers, the median sample size over time 192 was also U-shaped, reaching a minimum at 9 months post-vaccination before increasing ( Fig To examine the impact of transmission intensity in the population on carriage prevalence  197 in the trial, we also ran 50 simulations of the vaccine trial in the lower transmission setting. As 198 in the higher transmission setting, all vaccine efficacies resulted in reductions in carriage 199 prevalence at all sampling times. The prevalence peak previously observed in infants was 200 delayed, due to the slower acquisition of non-serotype-specific immunity in a lower 201 transmission setting (Fig 1). Thus, the prevalence trajectories in controls and vaccinees 202 followed non-monotonic trajectories in both infants and toddlers (Fig 3A-C). In the infant 203 arms, the kink in the prevalence trajectory between 9 and 12 months post-vaccination was 204 due to the change in age-specific contact patterns as the participants aged into the next age 205 group (Fig 3A-C, Table S1). 220 period 3 to 9 months post-vaccination, and reached similar minimums. In the infant arms, the 223 total sample size remained close to the minimum until the end of the 24-month follow-up 224 period. In the toddler arms, the median sample size increased slightly near the end of the 225 follow-up period. However, this rebound was considerably smaller than in the higher 226 transmission setting, and the median sample size at 24 months post-vaccination was roughly 227 five-to six-fold smaller. The sample sizes for the infant and toddler arms were more similar 228 than in the higher transmission setting, particularly for later sampling times (Fig 3D-F). 229

230
Using a computational, individual-based transmission model of pneumococcal carriage, 231 we estimated that a vaccine that enhances the immune response by an amount 232 corresponding to 3, 5, or 10 carriage episodes could reduce age-specific carriage prevalence 233 up to 7%, 10%, and 17%, respectively, compared to control in a setting similar to that of the 234 wSP vaccine trial in Kenya, but that the magnitude of the reduction would depend strongly on 235 the age at which participants were sampled. We found, however, that larger reductions could 236 be observed if the same trial were performed in infants, in a lower-transmission setting, or 237 both. Altogether, this analysis indicated that an infant trial conducted in a lower-transmission 238 setting for a vaccine simulating 3, 5, or 10 exposures could be adequately powered with fewer 239 than 800, 330, or 110 participants respectively, if the sampling window were chosen to be 15 240 to 24 months post-vaccination. Suboptimal choices of setting, age group, and sampling time 241 could multiply the required sample size by a factor of ten or more.
used to explain serotype diversity and explore serotype replacement following the introduction 244 of conjugate vaccines. With modifications, this model is also well suited to address our 245 modeling questions, because it incorporates many processes, epidemiological and 246 immunological, that complicate the relationship between the efficacy of a vaccine believed to 247 reduce carriage duration but not risk of acquisition, and its effect on carriage prevalence. Our 248 extensions-an algorithm to fit the model to specific epidemiological settings and the ability to 249 randomize trial participants to different vaccine interventions-allow this model to be used for 250 vaccine trial planning. 251 Our simulated vaccine trials show that sampling time and participant age greatly influence 252 the number of participants needed to detect a protective effect of a vaccine whose effect is 253 accelerating the development of immunity against carriage duration, as the wSP vaccine and 254 perhaps other protein-based vaccines targeting carriage are expected to do. Across different 255 combinations of vaccine efficacies and participant ages, the required sample size reached a 256 minimum approximately 9 months post-vaccination before rebounding in later months. This 257 favorable sampling time is consistent with simulation results by Scott et al., who explored 258 similar questions, but more generally and for vaccines whose primary effect is on acquisition 259 rather than duration, and using a compartmental transmission model [15]. This timing is also 260 consistent with what Auranen et al., who explored pneumococcal trial design questions with a 261 Markov transition model, suggest: waiting at least twice the average carriage duration after 262 immune response before sampling [32]. 263 In our simulations, the U-shaped trajectory of sample size over the follow-up period 264 indicates a window of favorable sampling times, when the sample size is relatively small as longer, when trial participants were younger, and when the transmission level was lower. In 267 these scenarios, natural immunity is weaker initially or develops more slowly, and thus 268 immune enhancement by the vaccine is more apparent. This intuition is what our simulation 269 study attempts to quantitate, in terms of sample size, for different trial conditions. 270 Certain model assumptions may affect our conclusions. Our formulation of vaccine 271 efficacy requires estimating the acquisition rate of exposure-dependent immunity. Direct 272 estimates of vaccine efficacy against carriage, when they become available, can be used 273 instead. We assume that the vaccine shortens only future carriage episodes, but not ones 274 already present at the time of vaccination. Since the intrinsic duration of the fittest serotype is 275 five months, this assumption would delay the vaccine's effect on carriage prevalence, and 276 thus, our reported favorable sampling times. This delay would affect infants more than 277 toddlers, as they are more immunologically naïve and experience longer carriage durations. 278 Auranen et al., in their study, report that sampling time is determined by the rate of clearance 279 rather than rate of acquisition, which reinforces the importance of determining whether a 280 vaccine accelerates the clearance of pre-existing carriage episodes [32]. Another important 281 assumption is that exposure, rather than age alone, is responsible for the progressive 282 shortening of carriage episodes as an individual gets older. If immune maturation due to 283 calendar age, rather than or in addition to increased exposure, actually reduces carriage 284 duration, then that would bolster the case for younger trial participants. Regardless of age at 285 vaccination, the favorable sampling windows will likely be shortened as well. Our simulation 286 framework can be easily updated should future evidence suggest revisiting these In its current form, our current simulation framework is already adaptable enough to 289 examine a variety of scenarios. The ability to tailor simulations to specific settings is 290 particularly useful-vaccine trials take place in countries with different age and serotype 291 distributions, and Phase I/II and Phase III trials of the same vaccine may be conducted in the 292 different locations. While we present results for a vaccine against carriage duration, we can 293 also model vaccine protection against acquisition, and specify whether a vaccine effect is 294 serotype-specific. The analysis presented here can be easily repeated, without changes to 295 the source code, for trials involving polysaccharide conjugate vaccines, which protect against 296 acquisition [4] and whose protection is serotype-specific [10], and novel vaccines with both 297 polysaccharide and protein antigens [33], which may elicit a combination of serotype-specific 298 and cross-reactive responses against carriage. The general population can also be 299 vaccinated. Hence, our framework can be used to simulate trials-such as those comparing 300 dosing schedules-that take place in countries with existing vaccination programs. In addition 301 to planning future trials, our simulation framework can be used to examine completed trials. 302 For completed trials with carriage endpoints that have not found a statistically significant 303 vaccine effect, such as a recent phase II trial of a protein and polysaccharide-based vaccine 304 in Gambian infants [33], simulation studies such as this can help assess whether inadequate 305 power is a compelling explanation. 306 The analysis presented in this paper does not consider the effect of vaccination on 307 carriage density or other factors (apart from duration) that would affect the infectiousness of a 308 person who is vaccinated yet still becomes colonized. More generally, we do not consider the 309 impact of vaccination on transmission at all in our simulations: simulated trial participants are 310 computationally isolated from other hosts to approximate an individually randomized trial in which the participants are a negligible fraction of the population. However, our current 312 framework can also simulate roll-outs of vaccination programs in the simulated population, 313 where there is transmission between individuals, thus allowing the indirect effect of 314 vaccination to be included. Vaccines with direct effects against transmissibility, possibly via 315 reducing bacterial density in the nasopharynx, can be incorporated into our framework as 316 well, with minimal modifications to the source code. 317

Mathematical model 319
Pneumococcal transmission dynamic model. This simulation study was based on a 320 published individual-based model of pneumococcal carriage that incorporates many of the 321 complexities relevant to our modeling questions [30]. Briefly, hosts are exposed to and can be 322 colonized by multiple serotypes through age-specific contact with others. Serotypes differ in 323 their mean duration of colonization in a naive host ("intrinsic duration"), which ranges from 20 324 to 150 days [19,20], and in their ability to prevent other strains from colonizing the same host. 325 These phenotypes are positively correlated-i.e. fitter serotypes have longer intrinsic 326 durations and are more likely to prevent concurrent colonizations-through their dependence 327 on a serotype-specific fitness parameter. Hosts acquire immunity through colonizations. 328 Clearing a colonization results in serotype-specific (anti-capsular) immunity that reduces risk 329 of acquisition of the same serotype. Each clearance, of any serotype, enhances non-330 serotype-specific immunity that reduces the mean duration of carriage episodes. wSP vaccine effect. The wSP vaccine was modeled as accelerating the acquisition of 332 non-serotype-specific immunity that reduces carriage duration. As in Cobey et al. [30], the 333 duration of a carriage episode is drawn from an exponential distribution with a mean given by 334 ! ", $ % = ! min + (! -− ! min ) exp −3$ % , 1 335 where " is the serotype carried, $ % is the number of cleared carriage episodes (of any 336 serotype), ! min is the minimum mean duration, and ! -is the intrinsic duration of serotype ". 337 The exposure-dependent development of non-serotype-specific immunity is captured in the 338 exponential decay term in Equation 1. Each cleared colonization is immunizing, but with 339 diminishing returns, and brings the mean duration closer to the minimum mean duration. For 340 a vaccinated individual, the mean duration is given by 341 ! ", $ % = ! ", $ % + $ 5 , 2 342 where $ 5 is a positive constant characterizing the strength of the vaccine effect. Thus, the 343 wSP vaccine can be thought of as boosting the non-serotype-specific immunity by an 344 additional $ 5 cleared colonizations, and we can express its efficacy in terms of "colonization 345 equivalents" or "c.e." We considered three different values of $ 5 : 3, 5, and 10. The duration of 346 each carriage episode was determined at the time of colonization, and hence, the vaccine did 347 not affect colonizations already present on the day of vaccination. For simplicity, we assumed 348 that full efficacy is achieved immediately upon receipt of a single dose. 349 Vaccine trials. To the original transmission model, we added the ability to simulate 350 vaccine trials. Each trial arm was characterized by the number of participants, the enrollment colonizations do not contribute to the force of colonization for the main population, but their 354 exposures and risk of colonization were equivalent to those of the same age in the main 355 population. This implementation design ensured that their colonization histories remain 356 representative of participants within the main population, while affording two advantages: 1) 357 We can have an arbitrarily large number of trial participants without skewing the 358 epidemiological dynamics of the population, and 2) participants can be "enrolled" simply by 359 birthing them into the simulation, without skewing the age structure of the population. 360 Alternatively, we could have achieved these properties by simulating a large enough 361 population such that the trial participants are a negligible fraction and thus do not create 362 appreciable herd immunity in the population-the case in most real-world individually-363 randomized vaccine trials. However, that approach would have been considerably more 364 computationally intensive. 365 Simulations. Simulations were initiated with hosts of different ages and no colonizations. 366 The number of hosts was kept constant throughout a simulation. Every simulation was run 367 first for 50 years to allow the age distribution of the population to stabilize, after which 368 colonizations were seeded in the population and the simulation was run for another 50 years 369 to allow the epidemiological dynamics to equilibrate. At this point, the simulated vaccine trial 370 was initiated. For simplicity, all participants were birthed into the trial on the same calendar 371 day. To reduce sampling noise, each trial arm had 5000 participants, 100-fold more than the 372 trial arms in the Kenyan wSP study [28]. The participants were followed for five years and the 373 carriage prevalence in each trial arm was recorded every 30 days. These carriage 374 prevalences were then used as "true prevalences" to calculate the sample size needed to 375 compare between arms, based on a two-sample test for equal proportions and assuming a combined size of both arms. All combinations of vaccine efficacies (3,5,10 c.e. and control) 378 and ages at vaccination (60 and 360 days) were represented in each simulated trial (for a 379 total of 8 arms), allowing us to control for transmission in the main population when 380 comparing between arms. For computational speed, the main population was set at 25 381 thousand individuals. For each parameter set, we conducted 50 simulations runs -enough so 382 that trends could be distinguished from stochastic variation between simulations, but not too 383 many as to require an unreasonable amount of computation time.  Table S1. The age structure in the model is 396 described in more detail in Text S1. We fixed the non-serotype-specific immunity acquisition 397 rate so the simulated age-specific carriage durations are consistent with the age-specific rates 398 of clearance in Kenyan toddlers estimated by Abdullahi et al. [40] (Fig S3). The serotype fitness parameters were fit to serotype-specific carriage prevalences from a cross-sectional 400 study in Kilifi from 2006Kilifi from to 2008, before the introduction of the conjugate vaccine PCV10. 401 We chose to fit using only pre-PCV10 data. Trying to reproduce changes in serotype 402 distribution due to PCV10 would have introduced additional complications, while being 403 unlikely to yield further insight into our modeling questions given that the wSP vaccine is 404 expected to act in a serotype-agnostic manner [41]. A mathematical description of the fitting 405 algorithm can be found in Text S2 and the fitted serotype fitness parameters are listed in 406 Table S2. 407 For the lower transmission setting, we used a smaller overall contact rate, so the 408 simulated carriage prevalence at 12 months of age resembles preliminary estimates from a 409 study in Indonesia [42], the proposed site for a follow-up wSP vaccine efficacy trial (Fig S3). 410 To facilitate comparisons between settings, we kept the same age distribution, age-specific 411 mixing pattern, and fitness parameters used in the higher transmission setting. A summary of 412 the model parameters and their values can be found in Table 1. 413

Sensitivity analyses 423
To isolate the effect of transmission intensity in our main analyses, we had used the same 424 age-specific mixing pattern-based on Kenya contact survey data [39]-in both the higher and 425 lower transmission settings. Real-world vaccine trials, however, will take place in the context 426 of different mixing patterns, or may be planned in the absence of reliable social contact data. 427 To examine the robustness of our findings to the pattern of age-specific mixing, we repeated 428 our analyses assuming random mixing between individuals, i.e., equal contact rate for all 429 pairs of individuals. We re-fit the model to the observed Kenya carriage survey data [31], and 430 ran a set of 50 simulations. With a random mixing pattern, there was a slightly higher carriage 431 prevalence in trial participants during the first two years of follow-up. However, the total 432 sample sizes, in both magnitude and trend across sampling time, remained similar to those analyses, we chose values that were small enough to allow simulations to finish reasonably 436 quickly, and reduced the effect of simulation variability by running multiple simulations and 437 considering sample median. To assess whether the sample median may be biased, we 438 performed univariate sensitivity analyses of the population and trial arm size. Specifically, 439 within the higher transmission setting, we varied population size between 10K, 25K, and 50K 440 individuals (not including trial participants), with the trial arm size fixed at 5K. We also varied 441 the trial arm size between 2.5K, 5K, or 10K participants, with the population size fixed at 25K. 442 Note that the middle values, a population size of 25K and a trial arm size of 5K, were the ones 443 used in the main analyses. Twenty-five simulations were run for each set of parameter 444 values. Varying the population and varying the trial arm size did not appreciably alter the 445 sample median of the simulated carriage prevalences (Fig S5). Larger population sizes led to 446 smaller variability between simulations, which is expected given the stochastic nature of 447 transmission in the model (Fig S5A, B). Larger trial arm sizes did not reduce variability, 448 suggesting that the epidemiological dynamics in the general population are driving the 449 variability in the trial arm prevalences, at least for the trial arm sizes examined (Fig S5C, D). 450    Nohynek H, Mäkelä PH, Lucero MG.  Text. Model age structure The age distribution and age-specific contact rate of hosts is important to consider in pneumococcal transmission modeling, since carriage prevalence varies with age [1,2], as does frequency of contact with other age groups [3,4].
The age distribution of the simulated hosts was matched to the 2015 age distribution in Kenya, based on data from the United Nations World Population Prospects [5]. The number of simulated hosts was constant, and for a fixed-sized population, we can set its age distribution by choosing the correct lifespan distribution: For a simulated host, the probability of living exactly n years is calculated as the difference in the number of n- year old people and n + 1-year old people, divided by the total number of people. For this method to be valid, the age distribution must be monotonically decreasing, i.e. there cannot be more people in an older age class as compared to any younger age class. This is the case for Kenya's age distribution in 2015. The World Population Prospects data was given in 5-year age classes, which we linearly interpolated to obtain 1-year age classes. The oldest age class in the data was 100 years or greater; in our model, we assume that the maximum lifespan is 101 years.

S2 Text. Model fitting algorithm
To simulate specific epidemiological settings, we implemented an algorithm that fit model parameters to given serotype-specific carriage prevalences, e.g. prevalences from survey data. In our model, the prevalence of each serotype is determined primarily by its fitness parameter and the overall contact rate shared by all serotypes. The fitness parameter can take values, possibly non-integral, from 1 to ! " , the number of serotypes, Lower values correspond to better fitness. Lowering the fitness parameter results in two phenotypic changes-longer colonization duration and enhanced competitive abilitythat both increase prevalence. Hence, there is a monotonic relationship between a serotype's fitness parameter and its expected carriage prevalence, and this allows us to tune the fitness parameters in a straightforward manner.
The algorithm iteratively updates its estimate of the serotype fitness parameters. Let the current estimate at the start of iteration # be denoted by the vector $ % , indexed by serotype. We run a simulation using $ % . For serotype &, let ' " % be its average prevalence over the last 25 simulation years, ' " be its observed prevalence, and ( " % = ' " % − ' , be