On the impact of masking and blocking hypotheses for measuring the efficacy of new tuberculosis vaccines

Data analyzed: the Brazilian BCG-REVAC clinical trials

BCG-REVAC consisted of a set of cluster-randomized trials involving more than 200,000 school-aged children in the Brazilian cities of Manaus and Salvador, whose principal aim was evaluating the effectiveness of BCG under different vaccination protocols.

The enrolled population of the study consisted of non-infected school children between 7 and 14 years old at the moment of randomization. Within this population, individuals presenting a positive BCG scar are separated from the rest, distinguishing, this way, the enrolled individuals who were vaccinated at birth from those who were not. Each group is then split into an intervention and a control group; individuals in the intervention group were vaccinated within the context of the trial. Summing up, there are 4 cohorts in each city: non-vaccinated (1), vaccinated after birth (2), firstly vaccinated at school age during the trial (3), and revaccinated, after a first dose applied after birth, in the trial too (4). Upon such classification of enrolled individuals in cohorts, the effectiveness of BCG vaccination strategies was measured by comparing the TB incidence rate within an end-point associated to active disease in the four cohorts, according to three different types of trials: Trial I: BCG at birth vs. no intervention (cohort 2 vs. 1). Trial II: BCG first dose at school age vs. no intervention (cohort 3 vs. 1). Trial III: revaccination at school age vs. first dose at birth only (cohort 4 vs. 2).

A model to describe BCG efficacy variation: masking, blocking and immunity waning

The six clinical trials conducted within the framework of BCG-REVAC study-three types of trials per two cities- output efficacies that span from 1% to 40% protection (see Fig. 1, red continuous lines). In order to explain this variability, we propose a model according to which the different protection levels found in each of the four cohorts in the study, schematically shown in Fig. 2, result from the interplay between the intrinsic vaccine efficacy, its temporal waning patterns, masking and blocking effects. These three mechanisms of vaccine protection shifts are ultimately responsible for vaccine’s performance variation, either in space or in time.

First of all, in absence of masking or blocking, a naive vaccinated individual will receive a protection level, right after vaccination, that we call e(0). As time after vaccination goes by, this protection level will wane up to e(t) < e(0), generally speaking. This implies that, if we deal with a population in which the incidence rate of new TB cases per unit time is equal to x; t years after vaccination, this rate is modified to (1 − e(t))x, provided that no additional effects take place. Taking that into account, a protective vaccine will have positive efficacy values e(t) ∈ (0, 1], being also possible for a (failed) vaccine to have a negative efficacy if it augments the disease risk among vaccinated individuals instead of reducing it. In our model, the time waning patterns of the intrinsic vaccine efficacy do not depend on the geographical area, but just on time since vaccination, which approximately is, in average, 4.5 years for school age vaccination (cohort 3) and 16 years for newborn vaccination (cohort 2), which implies the consideration of two intrinsic efficacy parameters: e(4.5) and e(16).

Besides vaccination, ES can also support protection against disease through the masking mechanism. The masking level, denoted by m, is a protection parameter formally equivalent to the intrinsic vaccine efficacy (thus verifying m ∈ (0, 1] for a protective effect, and negative otherwise), whose effects are suffered by initially naive, non-vaccinated individuals subject to ES. Thus, in principle, the longer the time an individual has been exposed to ES—i.e., the older the individual is at the moment of observation—the higher is the masking-related protection she might show. Masking is also a geography-dependent effect, since it depends on ES, which forces us to consider two masking parameters: m^M for Manaus and m^S for Salvador. The dependence of these parameters on age cannot be resolved, since all the cohorts analyzed in the study have approximately the same age.

Additionally, if e(t) describes the protection provided by the vaccine to a naive individual in absence of masking or blocking t years after vaccination, we also need to describe how this protection is modified if the vaccine is applied to non-naive subjects. If an individual’s immune system has been stimulated prior to vaccination (either by masking like in cohort 3, or by a previous vaccine, like in cohort 4 before the second dose), and consequently she is partially protected against the disease, it is unrealistic to assume that the full effect of the new dose is additive (Andersen & Doherty, 2005). Instead of that, our model considers that a vaccine dose applied on a previously protected individual will contribute, at most, up to resetting the initial protection levels e(0), provided that no blocking of the vaccine takes place. In cohort 3, this implies that, right after the school age vaccination, if the vaccine is not blocked (b = 0, see below), it will have a protective effect e′ that will be concurrent with the masking protection m so as to reduce the disease risk from [1 − m]x to [1 − e′][1 − m]x. Our estimation of e′ comes from assuming that such disease risk must equate what we would observe if a vaccine of full efficacy were applied on naive individuals, and observed 4.5 years later: (1)${\displaystyle \left[1-e\prime \right]\left[1-m\right]x=\left[1-e\left(4.5\right)\right]x\to e\prime =\frac{e\left(4.5\right)-m}{1-m}.}$ Similarly, the school-age dose at cohort 4, will add to the protection provided by the newborn dose e(16), diminishing the disease risk from [1 − e(16)]x to [1 − e″][1 − e(16)]x. To estimate e″, we assume that, if the second vaccine is not blocked, the disease risk achieved by both vaccines together [1 − e(16)][1 − e″]x is equivalent to the disease risk reached by the same vaccine, if applied on unprotected individuals, 4.5 years after vaccination: (2)${\displaystyle \left[1-e\left(16\right)\right]\left[1-e″\right]x=\left[1-e\left(4.5\right)\right]x\to e″=\frac{e\left(4.5\right)-e\left(16\right)}{1-e\left(16\right)}.}$ Finally, vaccine intrinsic efficacy can be blocked by prior ES; an effect that we model through the blocking probability b ∈ [0, 1], where b = 0 means that no blocking appears, while b = 1 stands for a totally blocked vaccine, meaning that vaccinated individuals would only have the protection level that they already had before vaccination. Blocking is also a geography-dependent factor, since it is considered a consequence of ES as well, which forces us to distinguish b^M and b^S for Manaus and Salvador, respectively. Unlike masking, blocking does not depend on the age of the individuals at the moment of observation, but on their age at the moment of vaccination. In this case we study cohorts vaccinated at two moments in life—at birth and at the beginning of the trials—being the first of these cases (the newborn vaccination) considered blocking-free, as it is assumed that when the vaccine is applied immediately after birth, there is no place for prior ES.

Taking all these effects into account, we are left with a set of six independent parameters $\overrightarrow{P}=\left\{e\left(4.5\right),e\left(16\right),m^M,b^M,m^S,b^S\right\}$ to describe the variability observed in the trials, either temporal or geographical, under the light of blocking and masking effects, concurrently. The temporal trends of the level of protection of each cohort are schematically shown in Fig. 3.

In the following, we will refer to this full model as model 1. In Fig. 2, we represent the variations on the disease rates provoked by each effect that takes place in each cohort according to model 1. Summing all the possible contributions to the development of active disease for each cohort, we derive the general disease rates characterizing each cohort of one city as follows: (3)${\displaystyle \begin{array}{c}{d}_1^l=\left(1-m^l\right)x\hfill \\ d_2^l=\left[1-e\left(16\right)\right]x\hfill \\ d_3^l=\left(1-b^l\right)\left[1-e\left(4.5\right)\right]x+b^l\left(1-m^l\right)x\hfill \\ d_4^l=\left(1-b^l\right)\left[1-e\left(4.5\right)\right]x+b^l\left(1-e\left(16\right)\right)x\hfill \end{array}}$ where the superscript indicates location, and x the incidence rate observed in the population.

From (3), it is immediate to derive the expressions for the observed efficacies $\bar{e}$ of each trial according to model 1, which read as: (4)${\displaystyle M_1:\left\{\begin{array}{c}{\bar{e}}_I^l=1-d_2^l∕d_1^l=\frac{e\left(16\right)-m^l}{1-m^l}\hfill \\ {\bar{e}}_{II}^l=1-d_3^l∕d_1^l=\frac{e\left(4.5\right)-m^l}{1-m^l}\left(1-b^l\right)\hfill \\ {\bar{e}}_{III}^l=1-d_4^l∕d_2^l=\frac{e\left(4.5\right)-e\left(16\right)}{1-e\left(16\right)}\left(1-b^l\right)\cdot \hfill \end{array}\right.}$ The system of Eq. 4 represents a full model for the vaccine efficacies observed during BCG-REVAC trials, which is based on the assumption that the sources of geographical variability for BCG’s performance are both masking and blocking effects. From the full model, two reduced versions can be conceived: a masking-free model (model 2 in the following) in which m^M = m^S = 0, and a blocking free model in which b^M = b^S = 0 (model 3). The efficacies associated to each trial, for models 2 and 3 straightforwardly read as follows: (5)${\displaystyle M_2:\left\{\begin{array}{c}{\bar{e}}_I^l=e\left(16\right)\hfill \\ {\bar{e}}_{II}^l=e\left(4.5\right)\left(1-b^l\right)\hfill \\ {\bar{e}}_{III}^l=\frac{e\left(4.5\right)-e\left(16\right)}{1-e\left(16\right)}\left(1-b^l\right)\hfill \end{array}\right.}$ (6)${\displaystyle M_3:\left\{\begin{array}{c}{\bar{e}}_I^l=\frac{e\left(16\right)-m^l}{1-m^l}\hfill \\ {\bar{e}}_{II}^l=\frac{e\left(4.5\right)-m^l}{1-m^l}\hfill \\ {\bar{e}}_{III}^l=\frac{e\left(4.5\right)-e\left(16\right)}{1-e\left(16\right)}\cdot \hfill \end{array}\right.}$ By considering these three models, our approach allows quantifying and comparing the plausibility of blocking and masking hypotheses to potentially explain the variation in BCG efficacy trials observed in the controlled setup conceived in the BCG-REVAC trials, taking into account the non-linearities associated to each mechanism, which play a central role in the derivation of Eqs. (4)–(6).

Models solution: parameters estimation and confidence intervals

In order to identify the set or sets of parameters yielding a best fit for the efficacies observed in BCG-REVAC trials, we compare the model prediction associated to any parameter set $\overrightarrow{P}$ to a set of empirical probability distributions derived from BCG-REVAC data. From each of the confidence intervals reported in Barreto et al., (2014) we build a two-piece normal distribution (Wallis, 2014) for each trial reported, centered in the reported values ${\left[{\bar{e}}_i^l\right]}_{\mathrm{BCG}-\mathrm{REVAC}}$ (for location l ∈ {Manaus,Salvador} and trial $i\in \left\{I,II,III\right\}$), and with asymmetric variances ${\left[{\sigma }_i^l\right]}_{\mathrm{BCG}-\mathrm{REVAC}}^{\pm }$ equal to one half the radius of the confidence intervals reported in Barreto et al., (2014), so preserving the confidence levels of the intervals reported (see Fig. 1).

Once the empirical distributions have been defined, for each possible set of parameters and for each of the six trials we define the Z-score associated to the model prediction as: (7)${\displaystyle Z_i^l\left(\overrightarrow{P}\right)=\left|\frac{{\left[{\bar{e}}_i^l\left(\overrightarrow{P}\right)\right]}_{\mathrm{mod}}-{\left[{\bar{e}}_i^l\right]}_{\mathrm{BCG}-\mathrm{REVAC}}}{{\left[{\sigma }_i^l\right]}_{\mathrm{BCG}-\mathrm{REVAC}}^{\pm }}\right|}$ where ${\left[{\sigma }_i^l\right]}_{\mathrm{BCG}-\mathrm{REVAC}}^{\pm }$ will take each of its two possible values depending on the sign of $\left({\left[{\bar{e}}_i^l\left(\overrightarrow{P}\right)\right]}_{\mathrm{mod}}-{\left[{\bar{e}}_i^l\right]}_{\mathrm{BCG}-\mathrm{REVAC}}\right)$. From $Z_i^l\left(\overrightarrow{P}\right)$, we define the corresponding p-values $p_i^l\left[Z_i^l\left(\overrightarrow{P}\right)\right]$ as the probability of the empirical distributions reproducing BCG-REVAC data to have a Z-score $\tilde{Z}$ so that $\left|\tilde{Z}\right|>\left|Z_i^l\left(\overrightarrow{P}\right)\right|$. This allows us to define the following likelihood function: (8)${\displaystyle L\left(\overrightarrow{P}\right)=\prod _{l,i}{p}_i^l\left[Z_i^l\left(\overrightarrow{P}\right)\right]}$ to maximize so as to identify the model’s parameters ${\overrightarrow{P}}^{\ast }$ more likely to yield the BCG-REVAC results. The global landscape of $L\left(\overrightarrow{P}\right)$ is explored using a hill-climbing algorithm designed to identify all possible local maxima in the space of parameters. Finally, a Levemberg-Marquardt algorithm is used to find a more accurate value of the global maximum, if the latter is unique.

In order to estimate the confidence interval associated to our model estimation, the following numerical procedure is performed. First, and starting from the maximum likelihood estimate ${\overrightarrow{P}}^{\ast }$, we move on each parameters’ axis until a value of L = 0.05 is reached in each case. We call this increment A_j (j ∈ [1, 6]) (see Fig. 4). These values are not symmetrical, again, and so we distinguish between $A_j^{+}$ and $A_j^{-}$. Using these asymmetric widths, we construct a two-piece normal distribution for every parameter [30], centered in ${\overrightarrow{P}}_o$ and having an asymmetric variance given by ${\sigma }_j^{\pm }=c{A}_j^{\pm }$, where c is a common modulation coefficient. Besides, the distribution is truncated at 1. Finally, we numerically estimate c by generating sets of points in the parameter space whose coordinates in each axis are obtained from the split normal distributions mentioned for an initial guess of c. Through an iterative process we search the value c = c^∗ for which a 95% of the points generated in the parameters space, yield efficacy estimations verifying $L\left(\overrightarrow{P}\right)>0.05$. Once we have found the optimal value of the scaling coefficient, the reported uncertainty of the j-th parameter corresponds to 95% CI given the distributions we have used.