Biological view of vaccination described by mathematical modellings: from rubella to dengue vaccines

1.   Introduction

Rubella is an acute exanthematous viral infection of children and adults, with the illness characterized by rash, fever and lymphadenopathy, resembling a mild case of measles. Although many infections are sub-clinical, this virus has the potential to cause fetal infection, with resultant birth defects and various forms of arthritis. After an attack of rubella, lifelong protection against the disease develops in most persons. Rubella can be a distasterous disease in early gestation (in general, the younger the fetus when infected, the more severe the illness) and can lead to fetal death, premature delivery and an array of congenital defects, known as congenital rubella syndrome (CRS). The live-attenuate vaccine was licensed for use in the United States in 1969, and the rationale for use of the vaccine is to prevent CRS by control of postnatal rubella. Currently, RA27/3 vaccine is in use ^[1]. The theoretical maximum risk of congenital rubella could be as high as $2\%$, in contrast to a $20\%$ or greater risk after maternal rubella in the first trimester of pregnancy ^[2].

Classic dengue fever is an acute self-limited illness with fever, headache, arthalgia, myalgia, rash, lymphadenopathy and leukopenia caused by four distinct serotypes of dengue virus (DENV), a mosquito-borne flavivirus, the principal vector being Aedes aegypti. These serotypes (types $1$ to $4$) constitute a distinct antigenic complex within Flavivirus genus, and close antigenic relationship among four serotypes, but cross-protection in humans is incomplete and short-lived. Dengue hemorrhagic fever, or severe dengue (SD), is distinguished from classic dengue by plasma leakage, and simultaneous thrombocytopenia, resulting death in a proportion of cases. By the fact that $90\%$ of SD occur during secondary infections, the immune status of the host is believed to play an important role in determining the course of dengue infection. The presence of non-neutralizing antibodies to heterologous DENV may form virus-antibody complex that enhance dengue virus entry into Fc receptor-bearing lymphoid cells, known as antibody-dependent enhancement (ADE). It is possible that certain DENV serotypes or genotypes may have differing propensities to induce SD. Because of the potential for SD disease in individuals sequentially exposed to wild dengue viruses, it will be necessary to simultaneously immunize against all or multiple serotypes ^[3].

This work is a reflection about tetravalent dengue vaccination being currently applied in several countries in order to control dengue epidemics, despite the possibility of occurring SD. To achieve this goal, the experience of rubella vaccination in recent past years is evoked: vaccination strategies were designed in order to reduce CRS. The reflection is based on biological aspects and assumptions done by mathematical modellings. Preliminarily, however, we discuss the misunderstanding about the risk of epidemics being maintained even below a threshold promoted by some vaccination schemes.

The paper is structured as follows. In section 2, we formulate the mathematical models of directly or indirectly transmitted infections. In section 3, we discuss the results of the models, considering well known biological and mathematical aspects of rubella vaccine, and application to recently introduced dengue vaccine. Conclusion is given in section 4.

2.   Mathematical model

There are many strategies of vaccination administrated in a community. For instance, routine vaccination (vaccination over all time), mass vaccination campaign (a vaccination-day to reach as many people as possible) and a mixture of both. In mathematical modellings, routine vaccination is formulated by considering a vaccination rate applied over all time, while mass (or pulse) vaccination is described by a proportion of susceptibles being vaccinated in a specified time (in general, using Dirac delta function). Routine and pulse vaccinations against rubella were analyzed in ^[4,5], assuming that incubation and infectious periods obey exponential distribution. However, other authors considered specified periods of incubation and infectiousness by using delay equations. For instance, Song et al. ^[6] considered pulse vaccination and two time delays and obtained conditions for the global attractivity of infectious-free periodic solution. De la Sen et al. ^[7] also considered two time delays and dealt with a model considering regular and impulsive vaccination, and a similar model was adapted and applied to Ebola infection ^[8].

Here we formulate a mathematical model to describe infections transmitted directly or indirectly, taking as examples rubella and dengue. The model is autonomous by dealing with constant vaccination rate over all time.

2.1.   Direct transmission – Rubella

When infectious individuals eliminate virus in the environment, this virus can infect susceptible individuals. This kind of air-borne interaction is called directly transmitted infection. Many infections are controlled by administrating vaccines. Vaccines can be partially effective (may not induce immunity, a primary failure) and/or lose induced immunity after a period of time (waning of immunity, a secondary failure).

We describe quantitatively directly transmitted infections under vaccination. With respect to the rubella vaccine, we discuss the question about the reactivation of attenuated virus and, as a consequence, potential risk of CRS.

A simple mathematical model divides a population into five classes, named susceptible ($S$), vaccinated ($V$), exposed ($E$), infectious ($I$) and recovered ($R^{\prime }$) individuals, and the total population is assumed to be constant, given by $N = S+V+E+I+R^{\prime }$. The dynamics of the model is given by

where the transmission rate is $\beta$, $\mu _{h}$ is the mortality rate of humans, and the incubation and infectious periods are $\gamma _{h}^{-1}$ and $\sigma _{h}^{-1}$, respectively. With respect to vaccination, $\nu$ is the vaccination rate, $q$ measures an imperfect vaccine, with $0\leq q\leq 1$, and the rates at which immunity induced by both vaccine and natural infection wane are $\phi _{v}$ and $\phi _{r}$, respectively.

Suppose that we have $q = 0$ (perfect vaccine) and $\phi _{v} = \phi _{r} = 0$ (absence of waning of induced immunity by vaccine and infection), in which case we can join $V$ and $R^{\prime }$ in one compartment, resulting in

where $R = V+R^{\prime }$. This classical SEIR model shows the well known behavior with respect to equilibrium points: the disease free equilibrium is stable if $R_{\nu } < 1$; otherwise, a unique endemic equilibrium is stable, where the reproduction number $R_{\nu }$ is given by

with $R_{0}$ being the basic reproduction number. Notice that $\nu$ decreases $R_{\nu }$, and $R_{\nu } = 1$ when $\nu ^{c} = \mu _{h}\left(R_{0}-1\right)$, resulting in the eradication of disease whenever $\nu \geq \nu ^{c}$. The change of stability from disease free equilibrium to unique endemic equilibrium is called forward bifurcation.

Suppose now that $q = 1$, a completely inefficient vaccine. In this case, we can join $S$ and $V$ as one compartment, defining $S^{\prime } = S+V$, and Eq (2.1) becomes

which describes an infection inducing temporary immunity in the absence of vaccination ($R_{\nu }$ does not exist). Hence, there is only the basic reproduction number $R_{0}$, given by Eq (2.3), resulting in a forward bifurcation. Now, consider vaccination in an infection inducing temporary immunity, that is, subtracting $\nu S^{\prime }$ in the first equation, and adding the same term in the last equation. In this case, we have

where the basic reproduction number $R_{0}$ is given by Eq (2.3). Again, there is only a forward bifurcation. Details about special cases ($q = 0$ and $q = 1$ with vaccination) can be found in ^[9] and references therein.

2.2.   Indirect transmission – Dengue

When infectious individuals need a vector or intermediate host to infect susceptible individuals, this interaction is called indirectly transmitted infection. Dengue virus is a vector-borne disease, transmitted by mosquito Aedes. Severe dengue (SD) can be understood as a consequence of huge viremia occurring by secondary heterologous infection or initial inoculation. In the first route, due to cross-immunity, a secondary heterologous infection with different serotype of dengue virus can induce antibody dependent enhancement (ADE) phenomenon, resulting in huge viremia. By the fact that ADE phenomenon amplifies virus titre, the available tetravalent dengue vaccine may increase cases of SD. In the second route, due to increased temperature, infected mosquitoes having high virus load can inject huge amount of virus when biting.

The increased cases of SD due to ADE can be explained by a simple model of interaction between human and mosquito populations. In order to model mosquito population, the variables $L$, $M_{s}$, $M_{e}$ and $M_{i}$ are used, which represent aquatic phase of mosquitoes, and susceptible, exposed and infectious mosquitoes. The parameters of mosquito population are the transmission rate $\beta _{m}$, the oviposition rate $\phi$, the mortality rates of aquatic forms $\mu _{l}$ and adults $\mu _{m}$, the transition rate from aquatic to adult $\sigma _{l}$, the extrinsic incubation period $\gamma _{m}^{-1}$, the proportion of females $q$ and the carrying capacity $C$. For human population, the variables are $S$, $E$, $I$ and $R$, which represent susceptible, exposed, infectious and recovered humans. Among humans, the parameters are the transmission rate $\beta _{h}$, the mortality rate of humans $\mu _{h}$, the intrinsic incubation period $\gamma _{h}^{-1}$ and the infectious period $\sigma _{h}^{-1}$.

In the modelling, the four serotypes of DENV must be accounted. However, we present only one route of infection, that is, we restrict to the dynamics of infection with one serotype to individuals already infected by another serotype. For the susceptibility modelling, dengue transmission among sub-populations of mosquitoes is described by

and transmission among humans is described by

The additional parameter $\psi$ is the increased susceptibility (if $\psi > 1$), and $S$ is now representing individuals susceptible to new serotype of DENV (or, vaccinated) but maintain immunological memory against other serotype previously exposed (or, immunity induced by vaccine), and mosquitoes are transmitting the new serotype of DENV.

Instead of ADE, there is possibility of SD due to high viral load inoculated by infectious mosquitoes. For instance, the model proposed by Esteva and Yang ^[10] evaluated the effect of increased temperature on the SD. That model considered two sequential classes of infectious mosquitoes, those with lower and higher burdens of viremia. In this case, dengue transmission among mosquitoes is described by

and transmission among humans is described by

where $S$ and $M_{s}$ are susceptible individuals to all serotypes of dengue virus. The additional parameters $\sigma _{m}$ and $\psi$ are, respectively, the high loading virus rate ($\sigma _{m}^{-1}$ is the average period of time elapsed to infected mosquitoes harbor heavy viremia) and increased infectivity of mosquitoes harboring high viral load. Additional variables $M_{g}$ and $E_{g}$ are mosquitoes harboring high viral load and humans inoculated with high viral load. Hence, individuals with SD are a fraction of $E_{g}$.

Models formulated considering ADE (see Eq (2.4)) and high inoculation of virus by mosquitoes (see Eq (2.5)) are based on the occurrence of huge viremia in infected humans. But, they have subtle difference with respect to increased viremia. In ADE modeling, small inoculation of virus by mosquitoes (normally not resulting in infection) can result in infection due to the action of ADE, hence, the probability of infection can be considered increased. However, in initial inoculation modeling, due to high level of virus harbored by some mosquitoes, probability of infection is increased in comparison to mosquitoes harboring lower level of virus.

3.   Results

Initially, we need to establish that all vaccination schemes are safe with respect to primary and secondary failures of vaccine. In other words, there is not occurrence of increased endemic level due to wrongly-designed vaccination scheme. Further, we address, based on the foregoing mathematical models, the question of reactivation of attenuate rubella virus as a risk factor of CRS, and the question of dengue vaccine causing SD.

When a vaccine is imperfect ($0 < q < 1$) and both vaccine and natural infection immunities wane with time ($\phi _{v} > 0$ and $\phi _{r} > 0$), many authors established that endemic equilibria may arise even when $R_{\nu } < 1$, which behavior was called backward bifurcation. Among them we cite Arino et al. ^[11], Chengjun and Wei ^[12], Dushoff et al. ^[13] Jing and Deming ^[14], Kribs-Zaleta and Velasco-Hernández ^[15], Li et al. ^[16], Safi and Gumel ^[17], Sharomi et al. ^[18]. Summarizing, mathematical models established that some vaccination schemes will not eradicate the disease even when the vaccination effort surpassed below the threshold ($R_{\nu } < 1$), and may occur increased incidence of infection (catastrophic event). Our aim is to deny these affirmations taking into account the biological aspects, neglected by authors.

Equation (2.1) has the disease free equilibrium point denoted by $P^{0}$ with values

which is stable if $R_{\nu } < 1$ (see Appendix A for details), where

with $R_{0}$ being given by Eq (2.3). The endemic equilibrium, denoted by $P^{\ast }$, is given by

where $S^{\ast }$ is

and $I^{\ast }$ is positive solution of

with the coefficients being

and $\beta _{1}$ and $\beta _{0}$ are given by

From Eq (3.2), by assuming $R_{0} > 1$ (otherwise vaccination is meaningless), notice that $R_{\nu } < 1$ whenever $\nu > \nu ^{th}$, where $\nu ^{th} = \left(\phi _{r}+\mu _{h}\right) \left(R_{0}-1\right) /\left(1-qR_{0}\right)$, with $q < 1/R_{0}$. Observe that if a vaccine is highly imperfect, or $1/R_{0} < q < 1$, any vaccination scheme diminishes the number of cases, but never eliminates the disease.

It is clear that $I^{\ast } = 0$ is a solution of Eq (3.3), resulting in the disease free equilibrium $P^{0}$, and other solution comes from $a_{2}I^{2}+a_{1}I+a_{0} = 0$. The coefficient $a_{0}$ can be rewritten as

and it is easy to verify that there is a unique positive solution $I^{\ast }$ for $R_{\nu } > 1$, resulting in the endemic equilibrium $P^{\ast }$. Then anomalous epidemiological behavior can occur when $R_{\nu } < 1$ ($a_{0} > 0$). Let us analyze the existence of two positive solutions for $R_{\nu } < 1$, in which case appears backward bifurcation. There are two conditions to occur two positive solutions: (ⅰ) $a_{0} > 0$ and $a_{1} < 0$ and (ⅱ) $a_{1}^{2}-4a_{0}a_{2} > 0$. We analyze the weak condition (ⅰ) since the condition (ⅱ) is much more stringent than condition (ⅰ) and involves hard calculations ^[19].

Let us analyze the condition (ⅰ), which is satisfied if $\beta _{1} < N\beta < \beta _{0}$. These inequalities are true if $\beta _{1} < \beta _{0}$, or $\beta _{1}-\beta _{0} < 0$, that is,

where $Q(q) = b_{2}q^{2}+b_{1}q+b_{0}$, with the coefficients being

Hence, condition (ⅰ) is satisfied if $Q(q) < 0$, which occurs if (i.a) $b_{1} < 0$ and (i.b) $b_{1}^{2}-4b_{2}b_{0} > 0$. (It is easy to show that $Q(q)\geq 0$, for $q\leq 0$ or $q\geq 1$.) If both conditions are satisfied, then $Q(q)$ has two positive roots $q_{-}$ and $q_{+}$, with $q_{-} < q_{+}$, and $Q(q) < 0$ when $q_{-} < q < q_{+}$. This range $\left[q_{-}, q_{+}\right]$, recalling that $q_{+} < 1$ and $q < 1/R_{0}$ (in order to have $R_{\nu } < 1$), can be restricted to

where $\min \left\{ q_{+}, 1/R_{0}\right\}$ stands for the minimum value between $q_{+}$ and $1/R_{0}$. We analyze the weak condition (i.a) since the condition (i.b) is much more stringent and complex than condition (i.a). The condition $b_{1} < 0$ is satisfied if

where

Notice that $\phi _{\nu } > 0$ is true if (i.c) $M(\mu _{h}) > 0$, or $0\leq \mu _{h} < \mu _{h}^{c}$, where

and (i.d) $\phi _{r} > \phi _{r}^{c}$.

Let us interpret biologically the weak conditions to the appearance of backward bifurcation, by assuming $\sigma _{h} = \gamma _{h}$, which is very reasonable. As numerical example, we use $\mu _{h} = 0.014$ $years^{-1}$ and $\sigma _{h} = 52$ $years^{-1}$ for rubella ^[20]. First, $\mu _{h}^{c}\approx \sigma _{h}\left(1+\sqrt{2}\right)$, then $\mu _{h} < \mu _{h}^{c}\sim \sigma _{h}$ is satisfied, resulting in $M(\mu _{h}) > 0$, and the condition (i.c) situates in the biologically acceptable range. Now we evaluate the second condition (i.d), $\phi _{r} > \phi _{r}^{c}$, by assuming that $\sigma _{h} = \gamma _{h}\gg \mu _{h}$. In this case, we have $\phi _{r}^{c}\approx \mu _{h}$, resulting in $\phi _{r}\gtrsim \mu _{h}$ for condition (i.d), which means that the average life span of humans ($\mu _{h}^{-1}$) must be higher than the average period of time to lose the immunity induced by natural infection ($\phi _{r}^{-1}$), which seems to be biologically acceptable. (When $\phi _{r} = 0$ the backward bifurcation does not exist, see Eq (2.2).) Again, by assuming $\sigma _{h} = \gamma _{h}\gg \mu _{h}$, we yield

that is, the average period of time to lose the induced immunity by vaccine ($\phi _{\nu }^{-1}$) must be higher than the period of time to lose the immunity induced by natural infection ($\phi _{r}^{-1}$), which is biologically unfeasible.

Above results were assured by considering weak conditions (ⅰ) and (i.a). If more stringent conditions (ⅱ) and (i.b) are considered, more restrictive relationships between parameters $\phi _{\nu }$ and $\phi _{r}$ must be satisfied. We stress the fact that the condition (ⅰ) is satisfied if $q_{-} < q < q_{m}$, that is, if a vaccine is very efficient ($q < q_{-}$) or inefficient ($q > q_{m}$), then the backward bifurcation does not occur.

Therefore, we concluded that the backward (or catastrophic) bifurcation does not occur in a biologically range of model parameters (see Yang and Raimundo ^[21] for a discussion in tuberculosis infection), and all vaccination schemes are safe with respect to primary and secondary failures of vaccine. However, other biological features such as CRS (rubella) and SD (dengue) must be taken into account when designing vaccination scheme.

3.1.   Rubella vaccination – Past experience

Every vaccination strategy against rubella has the eradication of this infection as the main goal, but the minimization of CRS consists a secondary goal. The USA immunization policy aims to eradicate rubella by vaccinating essentially all girls and boys before they enter pre-school kindergarten or elementary school (this is a usual vaccination scheme applied to all other childhood infections). In the UK, however, the immunization strategy has been to vaccinate girls, and only girls, in the age range 10-15 years, thus taking advantage of the immunity acquired by natural infection in the first years of life. Another immunization scheme is combining such coverage with vaccination of as many 1- to 2-year-old boys and girls as possible ^[22].

The major concern about vaccine against rubella infection is the possibility of live-attenuate vaccine resulting in CRS in pregnant women. In Poland, $36$ unvaccinated but infected women were follow-up, and they concluded that CRS is a serious problem ^[23]. Immunization policies were elaborated taking into account the theoretical maximum risk of congenital rubella being as high as $2\%$ ^[2], $1.7\%$ ^[24] and $0.2\%$ ^[25].

In order to evaluate the risk of CRS due to vaccine, we consider a model described by Eq (2.2). In the analysis of the model, the rates $\mu _{h}$, $\gamma _{h}$, $\sigma _{h}$ and $\beta$ are assumed constant, but $\nu$ is allowed to vary with age, or, $\nu = \nu (a)$, where $a$ is the age of susceptibles being vaccinated. The vaccine is not applied to over all ages, hence the age depending vaccination rate for rubella infection is

where $\theta \left(x\right)$ is the step function, that is, $\theta \left(x\right) = 1$, if $x\geq 0$, otherwise, $0$; and $\nu _{1}$ and $\nu _{2}$ are constant vaccination rates, respectively, in the age intervals $\left[\alpha _{1}, \alpha _{2}\right]$ and $\left[\alpha _{3}, \alpha _{4} \right]$. Based on the SEIR model, different vaccination strategies are assessed.

In age-structured modelling, the time derivative becomes

If the contact rate is age-dependent, that is, $\beta = \beta (a, a^{\prime })$, where $a$ is the age of susceptibles and $a^{\prime }$ is the age of infectious, the time-dependent force of infection $\lambda (a, t)$ is defined by

(If $\beta$ is constant, the force of infection $\lambda$ is defined by $\lambda = \beta I$, where $I = \int_{0}^{\infty }I(a, t)da$.)

The steady state (letting zero to the partial time derivative, $\partial /\partial t = 0$) force of infection $\bar{\lambda}(a)$ of model (2.2) is defined by

where $\bar{I}(a)$ is the steady state age-structured infectious individuals.

One of the main results about vaccination modellings was the increase in the average age at first infection, resulting in increased number of CRS. Being the reduction of CRS an indirect result (the main goal is the eradication of rubella) of such immunizations, two parameters must be evaluated ^[22]. The first is the steady state average age at which individuals acquire infection (shortly, average age at infection), denoted by $\bar{A}$, which is defined by

where $\bar{S}(a)$ is the steady state age-structured susceptible individuals. The second is the ratio between the numbers of cases of CRS after and before immunization in the age interval $\left[a_{1}, a_{2}\right]$ (shortly, relative risk of rubella) at steady state, denoted by $\bar{\rho} \left(a_{1}, a_{2}\right)$, which is defined by

where $r(a)$ is the risk of CRS proportional to age-specific fertility (an example of such function can be found in ^[26]). The steady states $\bar{\lambda}_{0}(a)$ and $\bar{S}_{0}(a)$ are, respectively, age-specific force of infection and age-distributed number of susceptible individuals just before introduction of vaccination, that is, the solution of model (2.2) letting $\nu = 0$.

Additional to the asymptotic average age at which individuals acquire infection, given by Eq (3.7), we can compare the time evolution with respect to the steady state $\bar{A}_{0}$ before the introduction of vaccination. Hence, the time dependent relative average age at infection $\varphi (t)$ is defined by

where $\bar{A}_{0}$ is obtained by changing $\bar{\lambda}(a)$ and $\bar{S} (a)$ in Eq (3.7) by steady states before the introduction of vaccine $\bar{\lambda}_{0}(a)$ and $\bar{S}_{0}(a)$.

Next, we describe the evaluation of vaccination schemes in order to evaluate the risk of CRS. The discussion is based on the possibility that vaccine could result in CRS when administrated in pregnant women.

3.1.1.   Mathematical models and vaccination schemes

The proposition and evaluation of risk of CRS due to vaccination are based on mathematical model considering constant contact rate $\beta$. Based on a mathematical model assuming constant $\beta$, Massad et al. ^[26] discussed vaccination schemes, as did many other authors (see references therein).

Assuming a constant contact rate $\beta$, thought of being an average value of $\beta \left(a, a^{\prime }\right)$ over ages $a$ and $a^{\prime }$, Anderson and May ^[22] obtained the steady state average ages at infection as $\bar{A}_{0} = 1/\bar{\lambda}_{0}$ and $\bar{A} = 1/\bar{\lambda}$, where $\bar{\lambda}_{0}$ and $\bar{\lambda}$ are the forces of infection (constant contact rate implies in age independent force of infection) before and after vaccination. By the fact that $\bar{\lambda}_{0} > \bar{\lambda}$, that is, the immunization always decreases the force of infection, a paradigm was established: any vaccination scheme always increases the average age $\bar{A}$ in comparison with $\bar{A}_{0}$.

One of a major consequence of such modellings is that the steady state risk of rubella $\bar{\rho}\left(a_{1}, a_{2}\right)$ increases in the age interval $\left[a_{1}, a_{2}\right]$ where women can be pregnant, for instance, $\left[18, 45\right]$ years. In another words, vaccination increases the number of cases of CRS, due to the fact that vaccination increases the number of women infected with rubella in comparison with the situation without vaccine.

Due to the paradigm $\bar{A} > \bar{A}_{0}$, many vaccination strategies were proposed, for instance, a unique dose, or two doses. The first vaccination scheme attempts to eradicate the rubella infection, while the second scheme aims to diminish cases of CRS. For this reason, the second dose is administrated among female youngs in the age interval $\left[14, 20\right]$ years, for instance. It is worth stressing the fact that increasing in the average age at infection with vaccination is a ubiquitous result from constant contact rate modelling.

In ^[20], from an age-structured model considering constant contact rate, the time dependent relative average age at infection $\varphi (t)$, defined by Eq (3.9), was obtained for different age intervals $\left[\alpha _{1}, \alpha _{2}\right]$ vaccinated. The relative average age $\varphi (t)$ always increases as a vaccination effort (described by $\nu$) increases.

Could more elaborated mathematical model and accumulated knowledge about CRS caused by vaccine change the vaccination policy?

3.1.2.   Changing the paradigm and vaccination policy

In ^[27,28], age-structured modeling was developed considering an age-structured contact rate $\beta (a, a^{\prime })$ given by

where $\Gamma (x)$ is the gamma function, and $b_{1}$, $b_{2}$, $b_{3}$ and $\beta _{0}$ are model parameters. These parameters were fitted from an age-specific force of infection of rubella, in São Paulo State, Brazil ^[29]. From these fitted parameters, it was shown that the average age at infection $\bar{A}$, given by Eq (3.7), does not always increase with vaccination. Indeed, if the lower bound $\alpha _{1}$ of the age interval $\left[\alpha _{1}, \alpha _{2}\right]$ is lower than the average age $\bar{A}_{0}$ before the introduction of the vaccination, then, this age displaces to higher ages. However, if the lower bound $\alpha _{1}$ is higher than $\bar{A}_{0}$, that is, vaccine is administrated to older children, then the average age displaces to lower values ^[28]. This shows that the paradigm that average age at infection with vaccine always increases is true only in constant contact rate modellings.

By the fact that the paradigm is not always valid, the vaccination scheme of two doses can be improved. The first dose must be applied in early aged children in order to eradicate the disease, and the second dose, among female youngs before pregnancy in order to avoid CRS. The first dose increases the average age at infection, but the second dose must be managed to decrease the risk of CRS. The reason behind this is that the age-structured contact rate modelling shows that the number of cases of CRS decreases if vaccine is administrated in elder individuals.

Additionally to the improved mathematical modellings, several studies followed up inadvertently vaccinated pregnant women to establish if live-attenuate vaccine may or may not cause CRS. In all these studies none case of CRS was found among unknowingly pregnant women vaccinated ^{[2,24,30,31,32,33,34,35,36,37,38,39]}.

Based on recent findings that vaccine does not generate CRS (from field observations), and that the average age does not always increase (from more elaborated modeling), the conclusion is that any vaccination strategy is safe and if the vaccine is administrated to very young children, the disease can be eliminated faster and easier ^[22] as other childhood infections, for instance measles. However, in ^[2,24,32,39], they stated that none case of CRS was reported among pregnant women, but due to the fact that the actual risk may not be zero, women known to be pregnant should not be vaccinated.

3.2.   Dengue vaccination – Current challange

Infection by one of the four DENV serotypes has been shown to confer lasting protection against homotypic re-infection, but only transient protection against a secondary heterotypic infection. The most commonly accepted theory about the occurrence of SD is known as the secondary-infection or immune enhancement hypothesis. This hypothesis implies that patients experiencing a second infection with a heterologous DENV serotype have a significantly higher risk for developing SD. Preexisting heterologous dengue antibody recognizes the infecting virus and forms an antigen-antibody complex, which is then bound to and internalized by immunoglobulin. Because the antibody is heterologous, however, the virus is not neutralized and is free to replicate once inside the macrophage. Thus, it is hypothesized that prior infection enhances the infection and replication of dengue virus in cells of the mononuclear cell lineage, through ADE ^[40,41,42]. Due to these dengue-specific complexities, vaccine development focuses on the generation of a tetravalent vaccine aimed at providing long-term protection against all virus serotypes ^[43].

A recombinant yellow fever-17D-dengue virus, live, attenuated, tetravalent dengue vaccine (Dengvaxia CYD-TDV by Sanofi Pasteur) has undergone extensive safety and immunogenicity assessment, and is currently in late phase of development ^[44]. This dengue vaccine was first registered in Mexico in December, 2015, and has demonstrated its efficacy against symptomatic, virologically-confirmed dengue during the active surveillance phase ^[45,46,47]. However, by analyzing phase Ⅲ clinical trial of Dengvaxia CYD-TDV, Halstead and Russell concluded that the relative risk of hospitalization for children vaccinated when they were $5$ years or younger was $4.9$ ^[48], and suggested that the possibility of vaccine related ADE may have occurred across all age groups. In response, Hadinegoro et al. ^[49] claimed that a higher relative risk may be only temporary. Recently, Sanofi released a press saying that "the analysis confirmed that Dengvaxia provides persistent protective benefit against dengue fever in those who had prior infection. For those not previously infected by dengue virus, however, the analysis found that in the longer term, more cases of severe disease could occur following vaccination upon a subsequent dengue infection" ^[50]. However, no deaths occurred associated to vaccine ^[51].

There are approximately five additional vaccine candidates under evaluation in clinical trials, including other live-attenuated vaccines, as well as subunit, DNA and purified inactivated vaccine candidates. Additional technological approaches, such as virus-vectored and virus-like particles (VLP)-based vaccines, are under evaluation in preclinical studies ^[43].

3.2.1.   Mathematical model – ADE and tetravalent dengue vaccine

We showed that imperfect vaccine does not generate catastrophic phenomenon. However, imperfect tetravalent dengue vaccine applied to naive individuals could cause SD when infected by circulating DENV due to ADE phenomenon.

Currently accepted hypothesis about the occurrence of severe dengue (SD) is the increased susceptibility, severity and infectiousness in the cases of secondary infections mediated by ADE. Based on this hypothesis, mathematical models were proposed to evaluate the occurrence of SD. For instance, Billings et al. ^[52] proposed and analyzed several serotypes of DENV model, but did not considered mosquito population, by assuming fast dynamics for mosquitoes. They obtained extremely complex dynamics, even chaotic dynamics. However, Esteva and Vargas ^[53] dealt with a model considering two different serotypes of dengue disregarding ADE, but included mosquito population (see Eq (2.4) as an example of secondary infection). The inclusion of mosquito population modulated the complex dynamics, and they obtained a unique steady state at which both serotypes coexisted.

Coudeville and Garnett ^[54] proposed a very complex model to describe all serotypes of dengue infections as did Billings et al. ^[52] with respect to increased susceptibility, but they included mosquito population, three doses of vaccination, waning of immunity induced by vaccine, and age. They introduced seasonal effects only on biting and vector birth rates, letting entomological parameters constant (we recall that mortality and extrinsic incubation period are strongly dependent on temperature ^[55]). Another aspect is the introduction of, besides the usual $\gamma _{m}^{-1}$, a second period of incubation (by delayed equation). They performed numerical simulations. The inclusion of seasonality (non-autonomous dynamics) resulted in the appearance of complex dynamics.

3.2.2.   Mathematical model – Role of inoculation

Yang et al. ^[56] showed that the biting activity increases with temperature as well as the longevity of female mosquitoes. The replication of DENV depends on the metabolic environment of tissues in its midgut ^[57], which may increase with increasing temperature. These behaviors may contribute to increase the potential transmission of DENV with temperature (temperature-depending entomological parameters showed increased endemicity in hot seasons ^[58]).

During low temperature seasons, DENV infection produces mild symptoms, and SD is not observed. Differently, SD is frequently observed during hot seasons. Based on this observation, it is possible to hypothesize that the amount of DENV in the midgut of mosquitoes is proportional to the temperature. Hence, by assuming that high temperatures increase hugely the DENV harboured by infected mosquitoes, model (2.5) can be used to evaluate the effect of temperature on the appearance of cases of SD. Under this hypothesis, imperfect tetravalent dengue vaccine does not play any role in originating SD.

3.3.   Discussion

Villar et al. ^[47] determined serotype-specific vaccine efficacy of a tetravalent dengue vaccine in children in Latin America: $50.3\%$ for serotype 1, $42.3\%$ for serotype 2, $74.0\%$ for serotype 3, and $77.7\%$ for serotype 4. According to Eq (3.4), these vaccine efficacies against four DENV serotypes are maybe in the range of occurrence of backward bifurcation. However, the tetravalent Dengvaxia CYD-TDV by Sanofi Pasteur may not induce protection as well as the natural infection, hence it is more likely that Eq (3.5) is not satisfied. Hence, backward bifurcation can not appear. Notwithstanding, the inefficient tetravalent dengue vaccine may cause SD through ADE phenomenon.

Let us consider dengue epidemics occurred at Belo Horizonte City, Minas Gerais State, Brazil, in $2013$ and $2014$, based on weekly recorded cases of dengue, SD and average temperature. In $2013$, occurred a huge epidemics, with $26$ SD in a total of $98,082$ cases of dengue, while a mild dengue epidemics followed in $2014$, with $26$ SD in a total of $3,152$ cases of dengue. The odds of SD is $31.1$ in $2014$ compared to $2013$. Based on the fact that predominantly DENV-1 and DENV-4 circulated in both years ^[59], and that annual average temperatures in $2013$ and $2014$ were $21.85^{o}C$ and $22.27^{o}C$, respectively, it seems that temperature could be a risk factor of SD associated or not with ADE. But, more detailed analysis must be done taking into account SD cases and temperature, also including more regions where dengue is endemic.

Currently, based on ADE phenomenon, many researchers claim that dengue vaccine must protect simultaneously against four serotypes. An open question is: how far is ADE responsible for severe cases of dengue, due to the fact that this is based on epidemiological findings? Notice that the ultimate effect of ADE phenomenon is amplification of initial inoculation, which, in theory, must occur in all periods of year (specially in sub-tropical regions) when two or more serotypes are circulating. In other words, SD must occur proportionally to the number of cases of dengue. Alternatively, is high inoculation of virus by infected mosquitoes another source of SD? In this case, SD must occur preferentially in hot seasons ^[10], that is, SD must occur in elevated temperatures disregarding the level of epidemics. A third hypothesis is ADE and high viral inoculation are acting together: severe cases of dengue is mediated by ADE, but only if the inoculation situates above a certain threshold ^[60]. In another words, high inoculation of virus encounters in ADE the possibility of amplifying the viremia, which could result in SD.

Other hypotheses about SD assume that dengue viruses, like all animal viruses, vary and change genetically as a result of selection pressures as they replicate in humans and/or mosquitoes and that there are some virus strains that have greater epidemic potential. Phenotypic expression of genetic changes in the virus genome may include increased virus replication and viremia, severity of disease (virulence), and epidemic potential ^[40]. Then previous models, described by equations (2.4) and (2.5), could take into account these factors.

Based on results from directly transmitted infections modelling, maybe currently proposed dengue vaccination scheme should be inappropriate. Most of phase Ⅲ clinical trials vaccinated individuals of $9$ years and over ^[45]. There are records about inefficacy of vaccine in children under $8$ years ^[61]. Based on these findings, if vaccine is applied to individuals $9$ years and older, there arises an undesirable effect provided by mathematical model: the average age at infection could be displaced to earlier ages ^[27,28] (this result comes from a directly transmitted infections, but the result should be applied to vector-borne infections). By the fact that SD is primarily a disease of children under the age of $15$ years, although it may also occur in adults ^[40], the displacement of the average age at infection to earlier ages must be taken into account if vaccine is applied to individuals over $9$ years old ^[61].

4.   Conclusion

Public health authorities postulated that pregnant women must not be vaccinated by assuming that vaccine (attenuate virus) may have similar theoretical risk of congenital rubella as wild virus (maximum of $2\%$). Following this recommendation, vaccination strategies were elaborated and evaluated theoretically by mathematical models in order to find an optimal scheme that minimizes the risk of CRS. However, following unknown pregnant women vaccinated inadvertently during decades, none case of CRS was observed. Additionally, the paradigm of the displacement of the first age at infection to elder ages by vaccination is valid only for mathematical models assuming constant contact rate. In summary, by the fact that rubella vaccine did not induce CRS and a suitable vaccination scheme does not increase cases of CRS, the conclusion is that any vaccination policy is admissible and useful, and the one that drive to the eradication must be chosen. Only for the sake of safety, pregnant women are advised to avoid the vaccination.

According to the ADE theory, elevated viremia and high concentrations of proinflammatory and immunomodulatory cytokines, including interleukin $10$, are associated with SD. Patients with dengue fever have lower peak of viremia titres than do those with SD, and patients who survive SD have higher concentration of circulating interferon than do those who die ^[41]. However, the initial inoculation theory, focusing only on the increased viremia, can explain SD due to high level of virus inoculation by infectious mosquitoes. For instance, Ong et al. found that a primary infection with DENV originated in SD ^[62], which shows that ADE is not essential to cause SD. Hence, there is not yet a clear understanding about the source of SD.

Based on what is known about the source of SD, it seems to be reasonable to vaccine with caution an imperfect tetravalent vaccine. Two aspects however must be taken into account:

(ⅰ) Currently, due to field observations regarding efficacy of vaccine, the vaccination is recommended to $9$ years and older ^[45], assuming that vaccine can induce ADE in children under $9$ years. Mathematical modeling of directly transmitted infections considering age-structured contact rate ^[27,28] showed that in a vaccination scheme of old children the average age at infection is displaced to earlier ages in comparison with the average age before the introduction of vaccine. Hence, if ADE indeed causes SD, this vaccination scheme ($9$ years and older) may not be suitable due to increased cases of infection with wild virus in lower aged children by displacement of the average age at infection to earlier ages. In this vaccination scheme, more children will have SD triggered by wild virus via ADE.

(ⅱ) if ADE alone does not lead to SD, further knowledge about biological aspects of SD must be incorporated to design vaccination policies. However, accepting that ADE and high viral inoculation trigger SD, imperfect tetravalent vaccine currently available could be applied during lower temperature seasons when circulating mosquitoes may harbor low virus load. However, suitable mathematical models must be formulated to assess vaccination scheme in seasonally varying environment (as for dengue vector control in ^[63]).

Acknowledgements

Thanks to anonymous reviewers for providing comments and suggestions, which contributed to improving this paper.

Conflict of interest

The authors declare there is no conflict of interest.

A.   Stability of disease free equilibrium $P^{0}$

Let us consider the equations corresponding to exposed and infectious individuals who are harbouring virus from Eq (2.1), which are written in matrix form as

where $x = \left(E, I\right) ^{T}$ is the vector of state variables, with $T$ standing for the transpose of a matrix. There are two distinct possibilities of writing vectors $f$ and $g$.

The first way is by construction vectors $f$ and $g$ as

The partial derivatives of $f$ and $g$ with respect to $E$ and $I$, denoted respectively by $F$ and $G$, evaluated at the disease free equilibrium $P^{0}$ are

where $S^{0}$ and $V^{0}$ are given by Eq (3.1). The inverse of the matrix $G$ and the matrix of next generation $FG^{-1}$ are

The characteristic equation corresponding to the next generation matrix $FG^{-1}$ is

where $R_{\nu }$ is given by Eq (3.2), and the spectral radius (dominant eigenvalue) is given by $\rho \left(FG^{-1}\right) = \sqrt{R_{\nu } }$. If $\rho \left(FG^{-1}\right) < 1$, then all eigenvalues corresponding to matrix $F-G$ have negative real part, and disease free equilibrium is locally asymptotically stable ^[64].

The second way is by construction vectors $f$ and $g$ as

The partial derivatives of $f$ and $g$ with respect to $E$ and $I$ evaluated at the disease free equilibrium $P^{0}$ are

where $S^{0}$ and $V^{0}$ are given by Eq (3.1). The inverse of the matrix $G$ and the matrix of next generation $FG^{-1}$ are

The characteristic equation corresponding to the next generation matrix $FG^{-1}$ is

and the spectral radius is given by $\rho \left(FG^{-1}\right) = R_{\nu }$.

We have two different expressions for the spectral radius. However, by applying the conjecture proved in ^[65], we concluded that $R_{\nu }$ is the reproduction number.