Novel Ordered Stepped-Wedge Cluster Trial Designs for Detecting Ebola Vaccine Efficacy Using a Spatially Structured Mathematical Model

Background During the 2014 Ebola virus disease (EVD) outbreak, policy-makers were confronted with difficult decisions on how best to test the efficacy of EVD vaccines. On one hand, many were reluctant to withhold a vaccine that might prevent a fatal disease from study participants randomized to a control arm. On the other, regulatory bodies called for rigorous placebo-controlled trials to permit direct measurement of vaccine efficacy prior to approval of the products. A stepped-wedge cluster study (SWCT) was proposed as an alternative to a more traditional randomized controlled vaccine trial to address these concerns. Here, we propose novel “ordered stepped-wedge cluster trial” (OSWCT) designs to further mitigate tradeoffs between ethical concerns, logistics, and statistical rigor. Methodology/Principal Findings We constructed a spatially structured mathematical model of the EVD outbreak in Sierra Leone. We used the output of this model to simulate and compare a series of stepped-wedge cluster vaccine studies. Our model reproduced the observed order of first case occurrence within districts of Sierra Leone. Depending on the infection risk within the trial population and the trial start dates, the statistical power to detect a vaccine efficacy of 90% varied from 14% to 32% for standard SWCT, and from 67% to 91% for OSWCTs for an alpha error of 5%. The model’s projection of first case occurrence was robust to changes in disease natural history parameters. Conclusions/Significance Ordering clusters in a step-wedge trial based on the cluster’s underlying risk of infection as predicted by a spatial model can increase the statistical power of a SWCT. In the event of another hemorrhagic fever outbreak, implementation of our proposed OSWCT designs could improve statistical power when a step-wedge study is desirable based on either ethical concerns or logistical constraints.

1 Model description and equations

Ebola transmission model
We constructed a simple SEIR model that captures Ebola virus disease (EVD) natural history, incorporating transmission scenarios in the community, treatment units and during burial ceremonies . Eight major compartments in the model depicted in figure 3 (main text) reflect the dynamics of the disease. The state parameters are described below, and superscripts 'c' and 'etu' indicate state that corresponds to a host living in the community or in an Ebola treatment unit respectively.

• S -Susceptible to infection by EVD;
• E -Infected by EVD but do not have any symptom nor infectious to others; • F c e and F etu e -Infected by EVD and have early or "dry" symptoms (fever, myalgia, headache, presence of oropharyngeal lesions, nausea, abdominal pain, and rash) but not yet infectious to others; • I c and I etu -Infected by EVD and present late or "wet" symptoms (vomiting, diarrhea, coughing, or hemorrhage from any site e.g. epistaxis, hematemesis, hematochezia, melena) and infectious to others; • F u -Dead bodies due to EVD and infectious to others during funeral practices; • R -Recover from EVD; • N -Total host population.
The full set of dynamic equations is given by the following set of ordinary differential equations, where key intervention measure's parameters are time depending, and Table 1 (main text) provides the full list of parameters with descriptions and values.
The above equations describe the transmission process inside a "patch" (a region). The parameters β c (t), β etu (t), and β f (t) are the transmission coefficients in the community, in Ebola treatment units and during funerals respectively. α −1 is the natural EVD incubation time, φ −1 is the time from onset of "dry" symptoms to onset of "wet" symptoms, r −1 is the duration of infectiousness for survivors (duration of "wet" symptoms), and γ −1 is the time from onset of "wet" symptoms to death. θ(t) is the probability of been admitted in ETU's upon any symptoms. α −1 h is the time to hospitalization, φ −1 h is the time to onset of "wet" symptoms for those hospitalized with "dry" symptoms, r −1 h and γ −1 h are respectively the time till recovery or death after hospitalization with "wet" symptoms, f −1 is the duration of funeral practices. δ 1 (t) and δ 2 (t) are the case fatality rates in the community and in ETU's respectively, and k(t) is the probability of proper funerals when an Ebola patient die within the community. I etu 1 can be interpreted as the newly infected with "wet" symptoms in ETUs, and I etu 2 as those who enter the treatment units after been infected with "wet" symptoms for some time.

Gravity Model
Human mobility between the "patches" (regions) is described by the gravity model: where m(i, j) is the mobility rate from regions i to region j, ρ is a scalar multiplied by the gravity matrix representing between-patch coupling, P (i) is the population density of region i, D(i, j) is the distance between 2 "patches", n is a power that determines the strength of the dependence of transmission rate on distance, and T r is the total number of patches or regions. We assumed that individuals in the compartments (E), (F c e ) and (I c ) may travel at rate m between the regions. Therefore the force of transmission is also a function of the migration rate and the following additional equations for transmission between "patches" (or regions) can be written: where the subscript 'i' indicates the state or parameter value in the i th patch (i = 1, . . . , T r ), and T r represents the total number of regions considered. The transmission model for EVD natural history coupled with a gravity model give the metapopulation model that captures the spatio-temporal spread of EVD where disease spreads within and between "patches" (regions).

Initial Dynamics of Infection
We first calibrated the transmission model to capture the early disease dynamics in Sierra Leone (SL). We fitted the model to the cumulative cases count of EVD in SL from May 19 2014 to September 16, 2014 using the least square method to estimate β c (t), β etu (t), β f (t), and θ(t). We assumed that prior to the scale up of the interventions measures these parameters remained constant, and also during that time period when an Ebola patient within the community died we considered his probability k(t) to receive proper funerals (that is the corpus is well handle without risk of infection to susceptible during funeral practices) to be zero.

Changing Dynamics of Infection over Time
To model the effects of intervention and behavioral changes on the disease's dynamic, we assumed that the transmission coefficients began to decrease exponentially after a more vigorous international response in September 2014 such that: where β c 0 , β etu 0 , and β f 0 are respectively the transmission coefficients in the community, in ETUs and during funerals before intervention date τ ; q c , q etu and q f are respectively the exponential decay rates from the initial transmissions to the final transmissions coefficients β c 1 , β etu 1 , and β f 1 . We refitted the transmission model to incorporate data from September 2014 to mid January 2015, and using the least square method we estimated q c , q etu , q f , β c 1 , β etu 1 and β f 1 (Table 1 main text).
We found the case fatality ratio δ(t) in SL for the months of September through December 2014 to be best fitted to a cubic regression curve We modeled the probability of hospitalization as a function of the cumulative number of beds in ETUs, the number of hospitalized patients, and the Ebola prevalence within the community. Given that beside the available beds in treatment units other factors such as distance to treatment units, and the infected individual's willingness to seek hospitalization are to be considered for the probability of hospitalization, when there were more available beds in treatment units than infected cases in the community, we corrected our value by the reported proportion of hospitalized cases as of June 2015.  Figure 10 in the main text shows a comparison of the model projection of the order and timing of first case occurrence at the district-level in Sierra Leone compared to the WHO reported data. We also used the model to project the order and timing of first cases occurrence at the chiefdoms level in Sierra Leone, we assumed that infection between districts ceased after the onset of international response in the country in late September 2014. Figure S1 below shows a color map of the spatiotemporal trend in infection cases at the chiefdoms-level in Sierra Leone. Figure S1: Spatial model's projection of order and timing of first case occurrence in the affected chiefdoms of Sierra Leone.
We finaly used the spatial model to forecast infectious cases at the districts-level, the best fit parameters ( Table 1 in the main) for the Ebola natural history model obtained at the country level were used for each district in the metapopulation model and we rescaled the transmissions coefficients to qualitatively match the districts-level data. Figure S2 shows the model projection of the cumulative cases in some of the main affected districts of Sierra Leone.

Model Validation
To test the extent to which the spatial model's ordering of first cases depended on the disease's parameters used in the Ebola transmission model, we simulated a hypothetical outbreak of smallpox, chickenpox and measles in the same regions of SL. We considered a simple SEPIR model used by Burcu et al. [10], where susceptible individuals (S) when infected are moved into the exposed compartment (E) where they are not yet infectious to others, after an incubation period they become prodromal (P) with reduced infectivity. After certain days the prodromal show sign of first symptoms where they are moved to the infectious compartment (I) and they are assumed fully infectious to others, afterward they can either recover or die. We drawn the model's parameter values from published literature (Table S1). We then coupled the transmission model for each disease to the gravity model (with our initial estimate values for the gravity parameters ρ and n) and we reran the spatial model to obtain new ordering of first case occurrence. We measured the correlation between the spatial model's ordering of EVD outbreak versus smallpox, chickenpox, and measles outbreak by evaluating the Spearman correlation coefficients.

Vaccine Design Simulation and Calculation of Statistical Power
As explained in the main text, we simulated and assessed the statistical power of different type of vaccine trial designs: the standard stepped-wedge cluster trial (SWCT) and our proposed "ordered stepped-wedge cluster trial" (OSWCT). In the SWCT, we assumed that each week one unvaccinated cluster was randomly chosen from the control group to receive vaccination until all clusters have crossed over from the control to the treatment arm. Whereas, in the OSWCT, we based the ordering on either the observed incidence data, or the projected order of first case occurrence in each district, or based on the districts with projected highest incidence by our spatially structured Ebola mode. More specifically for instance for the peak-OSWCT (similar for first-oswct or data-OSWCT) at each time step of the trial we considered at most the four unvaccinated clusters that are predicted to have the highest incidence at that week, and among those four we randomly assigned one cluster to receive vaccination and the remaining 3 clusters served as control group. This step was repeated until all clusters crossed over from control to treatment group.
For all designs we assumed that the trial would be conducted in hight risk individuals such as healthcare workers and other front line workers, and that each simulated cluster had the same size, CS, or same number of individuals. To ensure that the clusters were geographically distinct with different infection risks, we chose one cluster to reside in each of the 14 districts of SL. To simulate the designs we followed Bellan et al [5], and we evaluated the cluster-level hazard by assuming that without effective vaccination a proportion p of the reported district cases would occur within the corresponding cluster. Thus we calculated the hazard in each cluster at each week of the trial as follow: H p (week) = (total reported district cases in a week) × p CS × 1 week We assumed that individuals in the clusters became infected following an exponential distribution waiting time, and each week each individual hazard was calculated as H p (week) × ε, where ε captures the expected variation in individual risk and ε ∼ ln N (1, 1). We assumed that clusters that crossed over to the treatment arm were fully vaccinated within a week with a vaccine efficacy of ve, and that it takes a delay di to develop protective immunity after vaccination. Consequently we assumed that the risk of infection was reduced by 1 − ve for individuals in the vaccinated clusters that had passed the delay di. Table S1 gives the values of the vaccine's parameters used in our baseline scenario, in sensitivity analyses we varied the trial start date, the proportion p, the vaccine efficacy ve, and the delay time until protective immunity di. For the trial start dates, we explored two scenarios, first when the trial begins early in the outbreak (before September 2014) before all districts were predicted to have their first case, and when the trial began later in the outbreak where districts have observed their first cases. For the later scenario, in the simulated trial we dropped the ordering design based on first case occurrence as all the districts would have already observed their first case. Bellan et al [5] have shown that when the simulated data for the SWCT is analyzed with a non-parametric based method (a permutation test) it preserved the targeted false-positive rate set at α = 0.05. We consequently analyzed our simulated data for all trial designs with a permutation test. From the simulated data we analyzed the regression model where Y it is the number of cases observed in cluster i at time t, X it is the vaccine status of cluster i at time t, CT it indicates if a cluster i has contributed or not cluster time from [t, t + 1) ,and P Y it is all person-time in cluster i from time t to t+1 (except during the time to develop protective immunity), and the model includes the autoregressive (AR1) covariance structure in C i for cluster-level variation. And the estimated vaccine efficacy was computed asve = 1 − exp(β vac ) To test the null hypothesis that vaccines has no effect on the outcome of vaccinated individuals, we kept the outcome value Y it for each cluster at each week and we randomly permuted 1000 times the time at which clusters received vaccination. We then re-analyzed the permuted data with the above regression model and for each permuted data we computed a correspondingβ p vac . To decide wether or not to reject the null hypothesis we calculated the Wald statistics [6] W p =β p vac − β vac se(β p vac ) , and we calculated the permutation p-value by obtaining the value that were more extreme than that of the original non-permuted data (2-side test).
We further assed the bias for each design by setting the vaccine efficacy ve = 0, and the figure bellow shows the distribution of the estimated vaccine efficacy for 1000 simulations of each design. We also calculated the corresponding type1error as shown in the main text. Figure S3: Distribution of the estimated vaccine efficacy, with true vaccine efficacy set to zero.
Similarly to the sensitivity analyses done in the main text for the effects of changes in trial start date, the expected proportion of district cases within clusters, and vaccine efficacy on the statistical power of each design; we also show here a corresponding analysis for changes in the time acquire immunity after vaccination. The vaccine efficacy was set to 90% with trial starting during the peak of the peak of the outbreak (weeks 28). Figure S4: Power by the required time to develop immunity to EVD after vaccination. SWCT: random ordering of clusters; data-OSWCT: clusters are ordered based on observed cases 2 weeks previously; peak-OSWCT: clusters are ordered based on the highest weekly projected incidence.