Modeling the transmission dynamics of varicella in Hungary

Abstract

Vaccines against varicella-zoster virus (VZV) are under introduction in Hungary into the routine vaccination schedule, hence it is important to understand the current transmission dynamics and to estimate the key parameters of the disease. Mathematical models can be greatly useful in advising public health policy decision making by comparing predictions for different scenarios. First we consider a simple compartmental model that includes key features of VZV such as latency and reactivation of the virus as zoster, and exogeneous boosting of immunity. After deriving the basic reproduction number $R_{0}$, the model is analysed mathematically and the threshold dynamics is proven: if $R_{0}\leq 1$ then the virus will be eradicated, while if $R_{0}>1$ then an endemic equilibrium exists and the virus uniformly persists in the population. Then we extend the model to include seasonality, and fit it to monthly incidence data from Hungary. It is shown that besides the seasonality, the disease dynamics has intrinsic multi-annual periodicity. We also investigate the sensitivity of the model outputs to the system parameters and the underreporting ratio, and provide estimates for $R_{0}$.

1 Introduction

The varicella-zoster virus is a highly contagious disease that affects a huge proportion of the population, consequently the varicella incidence is of a similar magnitude to the number of births. Although most people contract the disease in their childhood, when the symptoms are generally mild, complications may occur during the infection. Furthermore, at an older age the risk of serious complications is significantly higher.

In many developed countries, varicella vaccination programs are already implemented. Originally one-dose programs were introduced, which have been replaced by multiple-dose vaccinations in some countries by now. In Hungary, vaccination is marketed for non-routine use, and it has been made available free of charge in a few cities in recent years. There are many country-specific studies regarding the effects and cost-effectiveness of the introduction of varicella vaccination, e.g. [1, 3, 4, 11]. However, there are hardly any studies about Hungary ([10] is a retrospective, descriptive study), where the introduction of varicella vaccination into the routine childhood vaccination program is being introduced. Given the actuality and the importance of this issue, here we summarize the challenges of such a modeling work, draw some conclusions from the qualitative and experimental study of simple compartmental models and devise a plan for comprehensive future work.

First we give a summary of the main issues of the modeling, then develop and investigate the qualitative properties of a simple autonomous compartmental model. Having performed data analysis on the varicella incidence reports in Hungary (2010–2017), we introduce seasonality (according to the academic year) in the model, and fit the system parameters to the data set. We also perform parameter sensitivity analysis, and having observed a strong underreporting, we study the sensitivity of system parameters to the underreporting ratio. Finally, we will give a plan of our research of age-structured models.

2 Challenges in modeling

Below we summarize the major challenges for building a comprehensive varicella model, some of them are specific to Hungary.

Latency of the virus and reactivation as zoster. Upon recovery from the varicella infection, VZV remains in the body in latent form. In general, the individual develops lifelong immunity to VZV. This immunity usually prevents the reappearance of varicella, however the immunity can wane over time, hence the virus may reactivate causing zoster (shingles). Zoster infected people are also infectious, but at a lower rate than varicella infected persons. Shingles cannot be passed from one person to another. However, VZV can spread from a person with active shingles to cause chickenpox in someone who had never had chickenpox or received chickenpox vaccine. The length and the efficacy of VZV immunity can show a wide inter-individual variability. Hence the modeling must take into account both the short period of infections (weeks) and long-time consequences (years), as well as the age-structure of the population.

The hypothesis of exogenous boosting. The waning immunity against VZV can be boosted if the individual has a contact with a VZV infected person. Although there is no clear picture concerning the degree of the boosting effect, the existence of the exogenous boosting seems to be valid [8]. Assuming exogenous boosting, it is reasonable that after introducing vaccination, the number of varicella cases decreases and consequently the zoster incidence temporarily increases [1, 11].

Underreporting. In Hungary, monthly reporting of varicella cases to the public health authorities is obligatory. Unfortunately, the varicella incidence appears to be much higher than reported, since the annual birth number is about 2.5-times higher than the reported varicella cases; and according to serological studies, these two values should be nearly equal [11]. Among others, the main reason is that not every child is taken to the pediatrician, as there is no effective medical treatment.

Seasonality. Available data also reflects a seasonal behavior in varicella incidence (see Fig. 2). It can be traced back to the high number of infected children; consequently the school term and vacation play an important role in the spread of VZV. To describe this phenomenon, time-dependent contact rates are needed in the model.

Lack of zoster data. Contrary to varicella cases, it is not compulsory in Hungary to report the zoster cases. Therefore, there is no available data related to zoster. One needs to make assumptions, based on studies from other countries.

Long term dynamics. Since we need predictions for many years ahead, an age-structured model should handle the transitions between age cohorts, which makes it more difficult than in models for single outbreaks, such as influenza with short-term behavior [7]. Demographic changes also need to be taken into consideration.

Vaccines are already present. Parents have the opportunity to buy the vaccine on the market in Hungary. Some cities have made the vaccine available for free for local children. Thus, a fraction of children have already been immunized. This certainly has some effect on the model, but in this paper we do not take it into account.

Age structure. Since the virus is highly contagious and appears mainly among young children, the transmission dynamics of the virus is largely age specific. Furthermore, varicella has more severe symptoms and higher risk of complications at an older age, and reactivation in the form of zoster also occurs at older age. Hence, age-structured models are necessary to capture these phenomena.

Vaccination efficacy and waning. Since varicella vaccination was licensed in the mid-80s in some European countries, the vaccine parameters are fairly reliable. In case of MMRV vaccine, 65% of the vaccinated population acquires full protection after one dose and 95% after the second dose. The vaccine-induced protection wanes in 15–20 years after one dose; while the two-dose vaccination provides lifelong immunity [11].

Cost-benefit calculations. In 2017, [10] gave a comprehensive study on the economic burden of varicella in Hungary using descriptive statistical methods. There are many uncertainties related to the introduction of VZV-vaccination. Hence, detailed dynamic model-based studies of the economic effects can be extremely useful.

In what follows we consider mathematical models that tackle some of these challenges, such as virus reactivation and zoster, exogeneous boosting, strong seasonality in the data, and the high underreporting.

3 Insights from a simple compartmental model

Based on the known models in the literature [1, 11], we use a simple compartmental system in our studies with the compartments representing the varicella disease states: Susceptible, Exposed, Infectious, Recovered, Susceptible to Zoster, Zoster, Zoster Immune. Maternal immunity is not taken into account in our model. Although the real situation is different (see Sect. 4), for the sake of simplicity we assume that the birth and death rates are equal $(d)$. More precisely, we assume that individuals in each compartment die with rate d, and also produce newborns with rate d (but all of these newborns appear in the susceptible compartment). Under this assumption, the total population is constant and a proportional model can be used where the total population is normalized to $n=s+e+i+r+s_{z}+z+r_{z}=1$, thus the model variables represent the fraction of the population being in the given epidemiological state at time t. We can illustrate the model with the compartmental diagram Fig. 1.

Figure 1

The compartmental system illustrating the VZV model

The model is as follows:

$$ \begin{aligned}& s' = d - \lambda s - d s, \qquad s_{z}' = -\sigma \lambda s_{z} + \zeta r - \eta s_{z} - d s_{z}, \\ &e' = \lambda s - \varepsilon e - d e,\qquad i_{z}' = \eta s_{z} - \kappa i_{z} - d i_{z}, \\ &i' = \varepsilon e - \gamma i - d i,\qquad r_{z}' = \kappa i_{z} - d r_{z}, \\ &r' = \gamma i + \sigma \lambda s_{z} - \zeta r - d r, \end{aligned} $$

(1)

where the variables $s,e,i,r,s_{z},z,r_{z}$ are functions of the time t measured in years, and $(\cdot)'$ represents time derivative. The force of infection, meaning the rate at which susceptibles acquire the infection is denoted by $\lambda = \lambda (t)=\beta (i(t) + \nu i_{z}(t) )$. Note that $\lambda =\lambda (t)$ is not a constant, as in our dynamical model the force of infection is changing in time, yet our system (1) is not a non-autonomous system as $\lambda (t)$ is expressed by other model variables. Newborns are assumed to be susceptible (hence the inflow into the s compartment is $d\cdot n=d$), then, one can become infected upon adequate contact with a varicella or zoster infectious person, with transmission rates β and νβ, respectively. Having been infected, individuals go through a non-infectious latent period, and then they will be infectious. Following the recovery, individuals acquire immunity to VZV. Immunity may wane, and then individuals become susceptible to zoster. One can either be boosted through exposure to VZV and regain immunity with efficiency σ or become zoster infectious through reactivation of VZV with the rate η. Zoster recovered individuals have lifelong immunity to VZV. The average length of time spent in the exposed, infectious, temporary immunity, and zoster states are $\varepsilon ^{-1}$, $\gamma ^{-1}$, $\zeta ^{-1}$ and $\kappa ^{-1}$, respectively.

3.1 The basic reproduction number $R_{0}$

The basic reproduction number $R_{0}$ is a key parameter regarding the level of contagiousness of the disease. It measures the transmission potential of the pathogen, and serves as a threshold determining if the disease can invade the population or not. It is the average number of secondary infections produced by a single infected individual introduced in a fully susceptible population. Now, let us give the basic reproduction number intuitively. The idea of this formulation lays on [12].

In case of VZV not only the primary varicella infections play a role in spreading the disease, but also the zoster patients can transmit the virus to susceptibles. Thus, the basic reproduction number can be calculated as the sum of the average number of new infections produced by a varicella infected individual, plus the average number of new infections produced by this individual upon virus reactivation as a zoster infected individual, factoring in the probability that the patient with the primary infection survives until VZV reactivates as zoster. We denote these two values by $R_{0,v}$ and $R_{0,z}$, respectively. Hence,

$$ R_{0}=R_{0,v}+R_{0,z}. $$

Consider a person being in the Infected state, who transmits the disease to susceptibles. The new infected individuals have to survive first the Exposed period and then get into Infected state. We assume that the waiting times in each compartment are exponentially distributed. Considering both the epidemiological and demographical movements, the average time individuals spend in the Exposed state is $\frac{1}{\epsilon +d}$. Thus, the average fraction surviving this period is $\frac{\epsilon }{\epsilon +d}$. Similarly, the average length of the time when individuals can infect others in Infected state is $\frac{1}{\gamma +d}$. The contact rate is β. Hence, $R_{0,v}$ can be intuitively calculated as the product of these values, i.e.,

$$ R_{0,v}=\beta \frac{\epsilon }{\epsilon +d} \frac{1}{\gamma +d}. $$

A person being in the Zoster state also transmits the disease. Infected individual have to survive not only the Exposed state, but also the Infectious, the Recovered and the Susceptible to Zoster periods. Same as in the previous case, the average fraction to survive these periods are $\frac{\gamma }{\gamma +d}$, $\frac{\zeta }{\zeta +d}$ and $\frac{\eta }{\eta +d}$, respectively. The average length of the time when individuals can infect others in this state is $\frac{1}{\kappa +d}$. It is assumed, that the infectiousness of zoster infected individuals is proportional to the infectiousness of the varicella infected individuals with a scaling factor ν. So $R_{0,z}$ can be calculated as the product of these values, the fraction ν and the contact rate β. Thus,

$$ R_{0,z}=\nu \beta \frac{\epsilon \gamma \zeta \eta }{(\epsilon +d)(\gamma +d)(\zeta +d)(\eta +d)} \frac{1}{\kappa +d}. $$

Finally, we obtain

$$ R_{0} = \frac{\beta \varepsilon }{(\gamma +d)(\varepsilon +d)} + \frac{\nu \beta }{(\kappa +d)} \cdot \frac{\varepsilon \gamma \zeta \eta }{(\varepsilon +d)(\gamma +d)(\zeta +d)(\eta +d)}, $$

(2)

and the very same expression can be derived also from the next generation matrix approach. To make a clear-cut in the formulas, we use the following notations: $\tilde{\epsilon }=\epsilon +d$, $\tilde{\gamma }=\gamma +d$, $\tilde{\zeta }=\zeta +d$, $\tilde{\eta }=\eta +d$ and $\tilde{\kappa }=\kappa +d$, and then

$$\begin{aligned} R_{0}=\frac{\beta \epsilon }{\tilde{\epsilon }\tilde{\gamma }}+ \frac{\nu \beta }{\tilde{\kappa }} \frac{\epsilon \gamma \zeta \eta }{\tilde{\epsilon }\tilde{\gamma }\tilde{\zeta }\tilde{\eta }}. \end{aligned}$$

(3)

3.2 Qualitative analysis of equilibra

To analyze system (1), we consider only the biologically feasible region

$$ \Delta =\bigl\{ (s,e,i,r,s_{z},i_{z},r_{z}) \in \mathbb{R}_{+}^{7}: s+e+i+r+s_{z}+i_{z}+r_{z} \equiv 1\bigr\} . $$

By standard arguments, one can see that solutions of (1) with non-negative initial data remain non-negative for all future time. Combining with $s'+e'+i'+r'+s_{z}'+i_{z}'+r_{z}' \equiv 0$, we obtain that the region Δ is positive invariant with respect to (1). Note also that the solutions with given initial values are unique. Hence, the model is well-defined both mathematically and epidemiologically.

3.2.1 The disease-free equilibrium

It is easy to see, that the point $P_{0}=(1,0,0,0,0,0,0)$ is an equilibrium—called disease-free equilibrium—of the system (1). The biological meaning of this equilibrium is, that there are only susceptibles and no infected individuals. The main result is as follows.

Theorem 1

The disease-free equilibrium is globally asymptotically stable if$R_{0} \leq 1$and unstable if$R_{0}>1$.

Proof

First note that the reproduction number can be written in the following form:

$$ R_{0}=\frac{\beta \epsilon }{\tilde{\gamma } \tilde{\epsilon }} \frac{\tilde{\kappa } \tilde{\zeta }\tilde{\eta }+\nu \gamma \zeta \eta }{\tilde{\kappa }\tilde{\zeta }\tilde{\eta }}. $$

Then the inequality $R_{0}<1$ is equivalent to

$$ \beta \epsilon \tilde{\kappa }\tilde{\zeta }\tilde{\eta }+\beta \epsilon \nu \gamma \zeta \eta < \tilde{\kappa }\tilde{\gamma } \tilde{\epsilon }\tilde{\zeta } \tilde{\eta }. $$

System (1) can be written in the form $x'=f(x)$, where $x=(s,e,i,r,s_{z},i_{z},r_{z})$ and $f(x)$ is the right side of the system.

To prove the stability of the disease-free equilibrium, following the method in [12] we construct a Lyapunov function in several steps. Let the first function be a linear combination of e and i, namely

$$ L_{1}=\epsilon e+ \tilde{\epsilon }i. $$

Then

$$\begin{aligned} \dot{L}_{1}&= \biggl\langle \frac{\partial L_{1}}{\partial x},f(x) \biggr\rangle = \bigl\langle (0,\epsilon,\tilde{\epsilon },0,0,0,0), f(x) \bigr\rangle \\ &=\epsilon e'+\tilde{\epsilon } i' = \epsilon (\lambda s - \tilde{\epsilon } e) + \tilde{\epsilon } (\epsilon e-\tilde{\gamma }i) = \epsilon \lambda s - \tilde{\epsilon } \tilde{\gamma } i \\ & = \epsilon \beta (i+\nu i_{z})s - \tilde{\epsilon } \tilde{\gamma } i \end{aligned}$$

Due to the structure of $L_{1}$ the term $\epsilon \tilde{\epsilon }e$ is eliminated but the term $\tilde{\epsilon }\tilde{\gamma }i$ remains in $\dot{L}_{1}$. Next, we eliminate this term. Thus, we introduce $L_{2}=\gamma L_{1}+\tilde{\epsilon }\tilde{\gamma }r$. Then

$$\begin{aligned} \dot{L}_{2} & = \gamma \dot{L}_{1} + \tilde{\epsilon } \tilde{\gamma } r' \\ & = \gamma \epsilon \beta (i+\nu i_{z})s - \gamma \tilde{\epsilon } \tilde{\gamma }i + \tilde{\epsilon } \tilde{\gamma } \gamma i+ \tilde{\epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} - \tilde{ \epsilon } \tilde{\gamma } \tilde{\zeta }r \\ & = \gamma \epsilon \beta (i+\nu i_{z})s + \tilde{\epsilon } \tilde{ \gamma }\sigma \beta (i+\nu i_{z}) s_{z} - \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta }r. \end{aligned}$$

Now, to eliminate the term $\tilde{\epsilon } \tilde{\gamma }\tilde{\zeta }r$, we take $L_{3} = \zeta L_{2} + \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta }s_{z}$. Its derivative is

$$\begin{aligned} \dot{L}_{3} & = \zeta \dot{L}_{2} + \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } s_{z}' \\ & = \zeta \gamma \epsilon \beta (i+\nu i_{z}) s + \zeta \tilde{ \epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} - \zeta \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } r - \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta }\sigma \beta (i+\nu i_{z})s_{z} + \tilde{\epsilon } \tilde{\gamma }\tilde{\zeta } \zeta r - \tilde{\epsilon } \tilde{\gamma }\tilde{\zeta } \tilde{\eta }s_{z} \\ &=\zeta \gamma \epsilon \beta (i+\nu i_{z}) s - d \tilde{\epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} - \tilde{ \epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{\eta } s_{z}. \end{aligned}$$

To eliminate the last term, we take $L_{4} = \eta L_{3} + \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{\eta } i_{z}$, thus we get

$$\begin{aligned} \dot{L}_{4} & = \eta \dot{L}_{3} + \tilde{\epsilon } \tilde{ \gamma } \tilde{\zeta } \tilde{\eta } i_{z}' \\ & = \eta \zeta \gamma \epsilon \beta (i + \nu i_{z}) s - \eta d \tilde{\epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} - \eta \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{ \eta } s_{z} + \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{ \eta } \eta s_{z} - \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{\eta } \tilde{\kappa } i_{z} \\ & = \eta \zeta \gamma \epsilon \beta (i + \nu i_{z}) s - \eta d \tilde{\epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} - \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{\eta } \tilde{\kappa } i_{z}. \end{aligned}$$

Now, let $L = \tilde{\zeta } \tilde{\eta } \tilde{\kappa } (\epsilon e + \tilde{\epsilon } i) + \nu L_{4}$ be the Lyapunov function, and we obtain

$$\begin{aligned} \dot{L} & = \tilde{\zeta } \tilde{\eta } \tilde{\kappa } \bigl( \epsilon e' + \tilde{\epsilon } i' \bigr) + \nu \dot{L}_{4} \\ & = \tilde{\zeta } \tilde{\eta } \tilde{\kappa } \bigl( \epsilon \beta (i + \nu i_{z}) s - \tilde{\epsilon } \tilde{\gamma } i \bigr) + \nu \eta \zeta \gamma \epsilon \beta (i + \nu i_{z}) s - \nu \eta d \tilde{\epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} - \nu \tilde{\epsilon } \tilde{\gamma } \tilde{\zeta } \tilde{\eta } \tilde{\kappa } i_{z} \\ & = (\tilde{\zeta } \tilde{\eta } \tilde{\kappa } + \nu \eta \zeta \gamma ) \epsilon \beta (i+\nu i_{z}) s - \tilde{\zeta } \tilde{\eta } \tilde{ \kappa } \tilde{\epsilon } \tilde{\gamma } (i + \nu i_{z}) - \nu \eta d \tilde{\epsilon } \tilde{\gamma } \sigma \beta (i+\nu i_{z}) s_{z} \\ & = (i + \nu i_{z}) \bigl( (\tilde{\zeta } \tilde{\eta } \tilde{ \kappa } + \nu \eta \zeta \gamma ) \epsilon \beta s - \tilde{\epsilon } \tilde{ \gamma } (\tilde{\zeta } \tilde{\eta } \tilde{\kappa } + \nu \eta d \sigma \beta s_{z} ) \bigr) \\ & = (i+\nu i_{z}) \bigl( B (R_{0} s - 1) - \tilde{\epsilon } \tilde{\gamma } \nu \eta d \sigma \beta s_{z} \bigr) \leq B (R_{0} s - 1) (i + \nu i_{z}), \end{aligned}$$

(4)

where $B = \tilde{\gamma } \tilde{\epsilon } \tilde{\kappa } \tilde{\zeta } \tilde{\eta }$. Then $\dot{L} \leq 0$ follows from $R_{0} \leq 1$, where equality holds only in two cases: $R_{0}=1$ and the system is at the disease-free equilibrium ($s=1$); or $i=i_{z}=0$. It can easily be seen that the latter relation holds for all t only if the system is at the disease free equilibrium. Thus, the largest invariant set in L̇= is $\{P_{0}\}$, and the first part of the theorem is proven by LaSalle’s invariance principle.

Concerning the case $R_{0}>1$, we have to observe that the characteristic polynomial of the Jacobian of system (1) at the disease-free equilibrium takes the form

$$ P(u)= \bigl(u^{5}+a_{4} u^{4}+a_{3} u^{3}+a_{2} u^{2}+a_{1} u+a_{0} \bigr) (-d-u)^{2}, $$

where $a_{4}=\tilde{\gamma }+\tilde{\zeta }+\tilde{\eta }+\tilde{\kappa }+ \tilde{\epsilon }$ and $a_{0}= (1-R_{0} ) \tilde{\gamma } \tilde{\zeta } \tilde{\eta } \tilde{\kappa } \tilde{\epsilon }<0$. Since $P(0)<0$ and $P(u)>0$ for u large enough, there must exist a positive zero, i.e., a positive eigenvalue of the Jacobian, that implies the instability of the disease-free equilibrium for $R_{0}>1$. □

3.2.2 The endemic equilibrium

A point $P^{*}=(s^{*},e^{*},i^{*},r^{*},s_{z}^{*},i_{z}^{*},r_{z}^{*}) \in \Delta$ is an equilibrium if and only if the coordinates of $P^{*}$ satisfies the following equations:

$$\begin{aligned} & d - \lambda s - d s = 0, \end{aligned}$$

(5)

$$\begin{aligned} &\lambda s - \varepsilon e - d e = 0, \end{aligned}$$

(6)

$$\begin{aligned} &\varepsilon e - \gamma i - d i = 0, \end{aligned}$$

(7)

$$\begin{aligned} &\gamma i + \sigma \lambda s_{z} - \zeta r - d r = 0, \end{aligned}$$

(8)

$$\begin{aligned} &{ -}\sigma \lambda s_{z} + \zeta r - \eta s_{z} - d s_{z} = 0, \end{aligned}$$

(9)

$$\begin{aligned} &\eta s_{z} - \kappa i_{z} - d i_{z} = 0, \end{aligned}$$

(10)

$$\begin{aligned} & \kappa i_{z} - d r_{z} = 0, \end{aligned}$$

(11)

We saw above, that there exists a unique disease-free equilibrium $P_{0}=(1,0,0,0,0,0,0)$, if $R_{0}\leq 1$. An equilibrium $P^{*}=(s^{*},e^{*},i^{*},r^{*},s_{z}^{*},i_{z}^{*},r_{z}^{*}) \in \Delta$ is called endemic, if not only the first element is positive. We can prove the following theorem:

Theorem 2

If$R_{0}>1$, then system (1) has a unique endemic equilibrium,

$$ P^{*}=\bigl(s^{*},e^{*},i^{*},r^{*},s_{z}^{*},i_{z}^{*},r_{z}^{*} \bigr)\in \Delta. $$

In addition, $s^{*}$is the unique solution in the interval$(1/R_{0},1)$of the equation$Q(s)=0$, where

$$ Q(s) = \frac{d \gamma \zeta \epsilon }{\tilde{\epsilon } \tilde{\gamma }} (1-s) s - \frac{d \tilde{\kappa } \tilde{\zeta } \tilde{\eta }}{\beta \nu \eta } (1 - s) (R_{0} s - 1) - \frac{\tilde{\eta }^{2} \tilde{\zeta }^{2} \tilde{\kappa }}{d\beta \nu \eta \sigma } (R_{0} s-1) s. $$

(12)

The other coordinates of$P^{*}$satisfy the following equations:

$$\begin{aligned} &e^{*}=\frac{d}{\tilde{\epsilon }}\bigl(1-s^{*}\bigr), \end{aligned}$$

(13)

$$\begin{aligned} &i^{*}=\frac{\epsilon d}{\tilde{\gamma }\tilde{\epsilon }}\bigl(1-s^{*}\bigr), \end{aligned}$$

(14)

$$\begin{aligned} &r^{*}=\frac{\gamma }{d}i^{*}- \frac{\tilde{\eta }}{d}s_{z}^{*}, \end{aligned}$$

(15)

$$\begin{aligned} &s_{z}^{*}= \frac{\tilde{\kappa }\tilde{\zeta }\tilde{\eta }}{(\tilde{\zeta }-\zeta )\beta \nu \eta \sigma } \bigl(R_{0} s^{*}-1\bigr), \end{aligned}$$

(16)

$$\begin{aligned} &i_{z}^{*}=\frac{\eta }{\tilde{\kappa }}s_{z}^{*}, \end{aligned}$$

(17)

$$\begin{aligned} &r_{z}^{*}=\frac{\kappa \eta }{\tilde{\kappa }d}s_{z}^{*}. \end{aligned}$$

(18)

Proof

Our aim is to show that there exists a unique positive solution $P^{*}$ of equations (5)-(11) whenever $R_{0}>1$. The idea to show that lays on [12]. Using (5) and (7) we obtain $d-ds^{*}-\tilde{\epsilon }e^{*} = 0$, which implies $e^{*}=\frac{d}{\tilde{\epsilon }}(1-s^{*})$. From equation (7), we have $\epsilon e^{*}-\tilde{\gamma }i^{*}=0$, so $i^{*}=\frac{\epsilon }{\tilde{\gamma }}e^{*}= \frac{\epsilon }{\tilde{\gamma }}\frac{d}{\tilde{\epsilon }}(1-s^{*})$. Hence, $e^{*}$ and $i^{*}$ can be calculated provided that $0< s^{*}<1$ is given.

From equation (10) we gain $\eta s_{z}^{*} - \tilde{\kappa } i_{z}^{*} = 0$ thus $i_{z}^{*} = \frac{\eta }{\tilde{\kappa }} s_{z}^{*}$. Combining this with equation (11), we get $r_{z}^{*} = \frac{\kappa }{d}i_{z}^{*} = \frac{\kappa \eta }{\tilde{\kappa }d} s_{z}^{*}$. Equations (9) and (8) imply $\gamma i^{*}-d r^{*} - \tilde{\eta }s_{z}^{*} = 0$, thus $r^{*} = \frac{\gamma }{d} i^{*} - \frac{\tilde{\eta }}{d} s_{z}^{*}$. Hence, if $s_{z}^{*}$ is given, $i_{z}^{*}, r_{z}^{*}$ and $r^{*}$ can be calculated.

Now, we only need to determine $s^{*}$ and $s_{z}^{*}$. Since we are looking for an equilibrium point in Δ, i.e., all derivatives are zero, the coordinates of $P^{*}$ make the value of L̇ zero, where L̇ is given in (4) (see (4) in the proof of Theorem 1). We obtain

$$ \bigl(i^{*} + \nu i_{z}^{*}\bigr) \bigl( \tilde{ \gamma } \tilde{\epsilon } \tilde{\kappa } \tilde{\zeta } \tilde{\eta } \bigl(R_{0} s^{*} - 1\bigr) - \tilde{\epsilon } \tilde{\gamma } \nu \eta d \sigma \beta s_{z}^{*} \bigr) = 0 $$

We are looking for an interior point of Δ, consequently, $i^{*},i_{z}^{*} > 0$. We gain

$$ \tilde{\kappa }\tilde{\zeta }\tilde{\eta } \bigl(R_{0} s^{*} - 1\bigr) - \nu \eta d \sigma \beta s_{z}^{*}=0, $$

and then

$$ s_{z}^{*} = \frac{\tilde{\kappa } \tilde{\zeta } \tilde{\eta }}{\nu \eta d \sigma \beta } \bigl(R_{0} s^{*} - 1\bigr). $$

(19)

It means, if we can find $0< s^{*}<1$, we can calculate $s_{z}^{*}$, too. For determining these two values, we are looking for further relations between $s^{*}$ and $s_{z}^{*}$. The idea is, $(i^{*}+\nu i_{z}^{*})$ appears in several equations, so let solve these equations for $(i^{*}+\nu i_{z}^{*})$ in terms of $s^{*}$ and $s_{z}^{*}$. From equation (5), we obtain

$$ i^{*}+\nu i_{z}^{*} = \frac{d}{\beta } \cdot \frac{1-s^{*}}{s^{*}}. $$

(20)

Multiplying equations (8) and (9) by ζ and ζ̃, respectively, and then adding them, we gain

$$ \gamma \zeta i^{*} - d \sigma \beta \bigl(i^{*}+\nu i_{z}^{*}\bigr) s_{z}^{*} - \tilde{\eta } \tilde{\zeta } s_{z}^{*} = 0. $$

(21)

Similarly, multiplying equations (7) and (6) by ϵ and ϵ̃, respectively, and then adding them, we obtain

$$ \epsilon \beta \bigl(i^{*}+\nu i_{z}^{*}\bigr)- \tilde{\epsilon }\tilde{\gamma }i^{*}=0, $$

and hence

$$ i^{*} = \frac{\beta \epsilon }{\tilde{\epsilon } \tilde{\gamma }}\bigl(i^{*}+ \nu i_{z}^{*}\bigr)s^{*}. $$

(22)

Substituting (22) into (21), we obtain

$$ \bigl(i^{*}+\nu i_{z}^{*}\bigr) d \sigma \beta s_{z}^{*} = \tilde{\eta } \tilde{\zeta }s_{z}^{*} - \gamma \zeta \frac{\beta \epsilon }{\tilde{\epsilon }\tilde{\gamma }}\bigl(i^{*}+\nu i_{z}^{*} \bigr)s^{*}, $$

and so

$$ i^{*} + \nu i_{z}^{*} = \frac{\tilde{\eta } \tilde{\zeta }s_{z}^{*}}{\frac{\gamma \zeta \beta \epsilon }{\tilde{\epsilon } \tilde{\gamma }} s^{*}- d \sigma \beta s_{z}^{*}}. $$

(23)

We found a second solution for $i^{*}+\nu i_{z}^{*}$. Note that, for all $0< s^{*},s_{z}^{*}<1$ the denominator of the fraction on the right hand side is not zero because of (19).

Now, let us combine (20) and (23)

$$ \frac{d}{\beta }\bigl(1-s^{*}\bigr) \biggl( \frac{\gamma \zeta \beta \epsilon }{\tilde{\epsilon } \tilde{\gamma }}s^{*}- d \sigma \beta s_{z}^{*} \biggr) = \tilde{\eta } \tilde{\zeta }s_{z}^{*}s^{*}. $$

(24)

Using the expression of $s_{z}^{*}$ from (19), we can write (24) in the following form:

$$ \frac{d}{\beta } \bigl( 1-s^{*} \bigr) \biggl( \frac{\gamma \zeta \beta \epsilon }{\tilde{\epsilon } \tilde{\gamma }} s^{*} - d \beta \sigma \frac{\tilde{\kappa } \tilde{\zeta }\tilde{\eta }}{d \sigma \beta \nu \eta } \bigl(R_{0} s^{*}-1\bigr) \biggr) - \tilde{\eta } \tilde{\zeta } \frac{\tilde{\kappa } \tilde{\zeta } \tilde{\eta }}{d\beta \nu \eta \sigma } \bigl(R_{0} s^{*}-1\bigr) s^{*} = 0, $$

which depends only on $s^{*}$. Now, we can define the quadratic form Q:

$$ Q(s) = \frac{d \gamma \zeta \epsilon }{\tilde{\epsilon } \tilde{\gamma }} (1-s) s - \frac{d \tilde{\kappa } \tilde{\zeta } \tilde{\eta }}{\beta \nu \eta } (1 - s) (R_{0} s - 1) - \frac{\tilde{\eta }^{2} \tilde{\zeta }^{2} \tilde{\kappa }}{d\beta \nu \eta \sigma } (R_{0} s-1) s. $$

Clearly,

$$ Q \biggl( \frac{1}{R_{0}} \biggr) = \frac{d \gamma \zeta \epsilon }{\tilde{\epsilon } \tilde{\gamma }} \biggl( 1- \frac{1}{R_{0}} \biggr) \frac{1}{R_{0}} > 0. $$

Furthermore,

$$ Q(1)=-(R_{0}-1) \biggl( \frac{\tilde{\kappa }\tilde{\zeta }\tilde{\eta }}{\nu \eta }+ \frac{(\tilde{\eta }\tilde{\zeta })^{2}\tilde{\kappa }}{d\beta \nu \eta \sigma } \biggr)< 0. $$

The intermediate value theorem guarantees that there exists $s^{*} \in (\frac{1}{R_{0}},1 )$ for which $Q(s^{*})=0$. Moreover, this zero position is unique on the interval $(\frac{1}{R_{0}},1 )$ due to the concavity of $Q(s)$.

Finally, it is easy to verify that $s^{*}=1-e^{*}-i^{*}-r^{*}-s_{z}^{*}-i_{z}^{*}-r_{z}^{*}$. The theorem is proved. □

3.2.3 Disease persistence

We introduce the basic notions of persistence [14]. Let X be a nonempty set and $\rho \colon X\to \mathbb{R}_{+}$. A semiflow $\varPhi \colon \mathbb{R}_{+}\times X\to X$ is called uniformly weaklyρ-persistent, if there exists some $\varepsilon >0$ such that $\limsup_{t\to \infty }\rho (\varPhi (t,x))>\varepsilon \text{for all} x\in X, \rho (x)>0$. Φ is called uniformly (strongly)ρ-persistent if there exists some $\varepsilon >0$ such that $\liminf_{t\to \infty }\rho (\varPhi (t,x))>\varepsilon \text{for all} x\in X, \rho (x)>0$. A set $M\subseteq X$ is called weaklyρ-repelling if there is no $x\in X$ such that $\rho (x)>0$ and $\varPhi (t,x)\to M$ as $t\to \infty$. System (1) generates a continuous flow on the phase space Δ.

Theorem 3

The susceptible population is always uniformly persistent. If$R_{0}>1$, then the disease uniformly persists in the population.

Proof

From the inequality $s'\geq d-(\beta +\beta \nu +d)s$ we find $\liminf_{t\to \infty }s(t)>d/(\beta +\beta \nu +d)$, hence the susceptible population is persistent.

Next, consider the persistence function $\rho (x)=L(x)$, $x \in \Delta$, which is a linear combination of the infected compartments $e,i,r,s_{z},i_{z}$. Consider the invariant extinction space of ρ, a collection of $x \in \Delta$ with $\rho (x)=0$, which is given by the subset of Δ with $s+r_{z}=1$ and other components are zero. On this extinction space, solutions converge to $P_{0}$. Recall that

$$\begin{aligned} \dot{L} & =(i+\nu i_{z}) \bigl( B (R_{0} s - 1) - \tilde{\epsilon } \tilde{\gamma } \nu \eta d \sigma \beta s_{z} \bigr), \end{aligned}$$

from which it follows that $P_{0}$ is weakly ρ-repelling whenever $R_{0}>1$, and by Theorem 8.17 of [14], we obtain weak persistence, and strong uniform persistence follows from Theorem 4.5 of [14], since Δ is a compact set.

If $e(t)\to 0$ as $t \to \infty$, then from $(i+r+s_{z}+i_{z})'=\epsilon e -d (i+r+s_{z}+i_{z})$ it follows that all the compartments $i,r,s_{z},i_{z}$ converge to zero as well, contradicting to the persistence of $\rho =L$. Therefore e is uniformly weakly persistent, then similarly as above, e is uniformly strongly persistent as well. Then from $\liminf_{t\to \infty }i(t)\geq (\liminf_{t\to \infty }e(t)/( \gamma +d)$, the uniform persistence is inherited to the i compartment and respectively to $r,s_{z},i_{z}$ as well. We found that all infected compartments are strongly uniformly persistent. □

4 Data analysis and model fitting

Annual Varicella incidence data for 20 years and monthly data since 2010 in Hungary were available to us (red curves in Fig. 2 show the incidence corrected by the fitted underreporting ratio $q=0.4$).

Figure 2

Varicella incidences: data (red) and fitted model (blue)

Due to the strong seasonality of varicella, we replaced the constant β in the system by a periodic function $\hat{\beta }(t)=\beta (1+b \cos (2\pi t -c))$ with $c=0.5$ chosen by a separate fitting process according to the academic year. The seasonal system with parameters d birth-death rate, β (transmission rate), b magnitude of seasonality was fitted to the monthly data. In addition, based on our former arguments in Sect. 2, the underreporting ratio (q) has to be also included into the fitting process.

Note that in our former studies [6], the parameters d and b were fixed, but these parameters are quite uncertain and it turned out that the model is very sensitive to these parameters. Hence here they are also estimated by fitting.

Our methodology is the following. First we assign reasonable values to all the other parameters in the model, based on epidemiological evidence. Due to the lack of available zoster incidence data from Hungary, related parameters were either adapted from the literature [1, 4, 11, 13] or assumed by expert opinion. Having most parameters fixed, we performed the fitting for the key parameters $q,d,b,\beta$. After that, we let all the parameters (be previously fixed or fitted) vary in reasonable intervals to assess their impacts in a sensitivity analysis.

It is known from epidemiological practice that the average length of time spent in the exposed, infectious, and infectious zoster states are $14, 7, 9$ days, respectively. Since these lengths (in years) are $\varepsilon ^{-1}$, $\gamma ^{-1}$, $\kappa ^{-1}$, respectively, we may assume $\varepsilon = 26$, $\gamma = 52$, and $\kappa = 40$.

The length of temporary immunity is estimated as 20 years, and hence $\zeta =0.05$. Similarly, the average length of being zoster susceptible is assumed 30–35 years, hence $\eta = 0.003$. These assumptions agree with the fact that zoster appears at elder ages.

Based on consultations with epidemiologists, the force of infection of zoster is much less then that of varicella, hence $\nu = 0.07$ is assumed. Similarly, we may assume that strength of boosting the immunity by zoster is somewhat weaker than by varicella, hence $\sigma = 0.7$. Our parameters are similar to those typically used in the literature [1, 4, 11, 13].

Finally, values of $(s,e,i,r,s_{z},z,r_{z})$ at any time are not known, but we may assume that varicella incidence is close to the equilibrium, hence initial values of the solutions of the seasonal system were taken equal to the endemic equilibrium according to the values of parameters.

Now, the fitting model is simple: the cumulative growth of $i(t)$ is measured by $\hat{i}(t)$ with $\hat{i}'(t)=\varepsilon e(t)$, and hence the monthly and annual incidences are modeled by $MM(t)=q( \hat{i}(t+1/12) -\hat{i}(t))$ and $AM(t)=q(\hat{i}(t+1) -\hat{i}(t))$, respectively. Obviously, $i(t), \hat{i}(t), MM(t)$ depend also on the parameters $q,d,b,\beta$. The model $MM(t)$ is fitted to the varicella incidence data given monthly in the period 2010–2016 (79 months, 5 months are missing).

Fitting was performed by the sophisticated and well-tested command NonlinearModelFit in Wolfram Mathematica 12.0, which can be applied to implicitly defined models such as parametric numerical solutions of differential equations (ParametricNDSolve), and it can measure the goodness of the fit. Default options and $\mathit{ConfidenceLevel} \rightarrow 0.95$ were used.

The goodness of the fit was measured by the adjusted $R^{2}=0.929$. The fitted values are shown in Table 1 (at significance level: 95%).

Table 1 Results of the parameter fitting

The asymptotic parameter correlation matrix is shown in Table 2.

Table 2 The asymptotic parameter correlation matrix

The result can be seen on Fig. 2. The monthly incidence data and fitted model $MM(t)$ can be found on the left side, while the right one contains the annual data and the fitted model $AM(t)$ as well as the corresponding autonomous model with the same parameters. Finally, we emphasize that although the seasonality is very strong, both the monthly and annual incidence models show a multi-annual periodicity. The yearly peaks have maxima approximately at every four years. This phenomenon is known in the epidemiology of varicella and the value agrees the practice. The same period can be obtained by the autonomous model.

5 Sensitivity analysis

5.1 Sensitivity to system parameters

Now, we present the results of our sensitivity of both the autonomous system (1) and its seasonal version to the system parameters. Although most parameters were fixed before fitting, here we study the sensitivity to all parameters.

Global sensitivity (all parameters may change simultaneously) of the system is investigated by the effective and widely accepted LHS/PRCC method (see [5]). We have also investigated the sensitivity using other methods, for example using the differentiable dependence of solutions of differential equation with respect to the parameters. These other studies gave very similar results to the LHS/PRCC analysis.

Concerning both the original autonomous and the seasonal versions of (1), we investigate the sensitivity of $\hat{i}(t+1) -\hat{i}(t)$ at $t=0,0.5,\ldots4.5,5$. The random sampling is taken from the intervals:

$$\begin{aligned} & d\in [0.0093,0.011]; \qquad \beta \in [700,790]; \qquad \epsilon \in [24,28]; \qquad \gamma \in [46,58]; \\ & \kappa \in [35,45];\qquad \zeta \in [0.04,0.06]; \qquad \sigma \in [0.6,0.8]; \qquad \nu \in [0.06,0.08]; \\ & \eta \in [0.0025,0.0035]; \qquad b \in [0.2,0.3]. \end{aligned}$$

The number of parameter samples is $N=500$. The initial value of the solutions at every parameter combination is the equilibrium of the autonomous system (1) with the fitted parameters. Note that the partial derivatives of them with respect to the parameters are sign-keeping on the given parameter-intervals.

PRCC coefficients are determined at $t=0,0.5,1,\ldots,5$. The dependence of most sensitive parameters with respect to the time at the seasonal system can be seen on Fig. 3. Coefficient with largest absolute values at every parameter are shown on Fig. 4. The values explain the choice of parameters $d,\beta,b$ in the fitting. Sensitivity is also high to ε and γ, but they could be determined from the epidemiological literature. For reference, we present the PRCC values of $\hat{i}(t+1) -\hat{i}(t)$ of the original autonomous system (1).

Figure 3

Dependence of PRCC values on time of new varicella patients at the seasonal version of system (1)

Figure 4

Maximal PRCC values $(t\in [0,5])$ of new varicella patients at the seasonal (left) and autonomous (right) systems

Finally, investigate the sensitivity of basic reproduction number $R_{0}$ of the autonomous system (1). The PRCC coefficients for $R_{0}$ are ($N=1000$) can be seen in Fig. 5. It is obvious that $R_{0}$ is highly sensitive to zoster related parameters, but since no zoster report are available, the uncertainties of their estimates explain the uncertainties of $R_{0}$. See the country-wise difference in [9]. We consider another aspect in the next section.

Figure 5

PRCC coefficiens for $R_{0}$ of the autonomous system (1)

5.2 Sensitivity to underreporting ratio

According to the previous section, varicella cases are likely to be seriously underreported in Hungary $(q \approx 0.4)$. The model fitting is coherent with what the serological studies suggest. In this section we investigate, how sensitive our model is to the ratio of the reported and total cases, i.e., we examine dependence of the basic reproduction number $R_{0}$ (see eqn. (2)) on this ratio q at the parameters fitted above.

Assuming that the population is at the endemic equilibrium in Hungary, the following equality holds:

$$ n_{V} / q = \gamma i^{*}, $$

where $n_{V}$ is the annual reported varicella incidence and $i^{*}$ is the endemic equilibrium of i and $1/\gamma$ is the mean length of the varicella infection. Taking the formula for $i^{*}$ from (12) and (14) and the fitted parameter values given in the previous section, we obtain the following relation between q and $R_{0}$

$$ \frac{n_{v}}{q}=-0.01700 \sqrt{ \biggl(\frac{0.2989}{R_{0}}-0.0296 \biggr){}^{2}+\frac{0.03606}{R_{0}}}-\frac{0.00508}{R_{0}}+0.01068. $$

The result is depicted in Fig. 6 with different values of $n_{v}$. We can see, that $R_{0}$ is very sensitive to the change of q at low values, such as around the fitted $q=0.405$. Considering $n_{v}=0.003658$ the mean annual incidence since 2010 to 2017, the corresponding $R_{0}=8.93$ (see the red curve). Estimating $n_{v}$ by the formula $n_{v}=q d p$, where p is the infected fraction of the total population, we obtain the value $p\approx 0.886$ what looks underestimated comparing to the observation from serological studies that 90–95% of the population becomes infected.

Figure 6

Relation between the underreporting ratio q and the basic reproduction number

On the other hand, $R_{0}=14.56$ (obtained from the fitted parameters) provides $n_{v}=0.003837$, that gives $p=0.9302$ of the Hungarian population is infected, agreeing with the above mentioned observation (see the blue curve).

Finally, considering $p\in [0.9,0.95]$ we obtain the estimation $10.13 \leq R_{0} \leq 20.45$ (turquoise bounded gray region).

Note that in the literature a wide variety of different $R_{0}$ values can be found for the VZV. In [9], the highest value is 16.91 (Netherlands) and the lowest is 3.31 (Italy). Knowing the parameters in these countries, estimation on the underreporting can also be given.

6 Conclusion

Motivated by the ongoing introduction of VZV vaccination into the routine immunization schedule, we gave an overview of the of the current varicella transmission dynamics in Hungary. After summarizing the modelling challenges, first we investigated the qualitative properties of the compartmental model (1), which is a variant of the model studied in [12]. Providing a similar analysis as in [12], we derived the basic reproduction number and proved its threshold property: for $R_{0}\leq 1$ we showed the global asymptotic stability of the disease-free equilibrium, while for $R_{0}>1$ a unique endemic equilibrium exists and the disease is uniformly persistent in the population. Since we had monthly data available from Hungary, showing strong seasonality in disease prevalence, we extended our model with seasonal (time periodic with period one year) disease transmission. We fitted the seasonal version of system (1) to the available data, and found that the strong seasonality of varicella infections and the underreporting are essential factors to be considered in varicella modelling. The model also reproduced the multi-annual cycles observed in varicella dynamics. We performed a sensitivity analysis to see which parameters are the most important in the output of the model. We found strong sensitivity of $R_{0}$ with respect to the underreporting ratio, which can be an explanation to the wide range of basic reproduction numbers estimated in different studies.

The long term aim of our research is to forecast and evaluate the impact of vaccination in Hungary. Introducing VZV vaccination is expected to significantly reduce the number of varicella infections, but it has some negative effects such as the average age at infection is increasing thus conferring higher risks of complications, and a temporarily increase in the number of zoster cases is also expected due to the reduced boosting effect if less VZV is in circulation. WHO recommendations that countries which decide to introduce the varicella vaccination assess the disease burden caused by varicella and continue surveillance after the introduction of vaccination [15], thus we plan to continue our work to this directions. Based on our simple model the global effects and strategic goals can be already visible. To build a realistic model which can be used to evaluate the impact of vaccination policies, our model should be significantly extended by an other modelling arm representing vaccinated compartments, the parameters and the specifics of the given vaccination strategy, seasonal effects, and detailed age structure with age specific parameters and contact patterns.