The Effect of Setting a Warning Vaccination Level on a Stochastic SIVS Model with Imperfect Vaccine ^†

Abstract

This paper deals with a stochastic Susceptible-Infective-Vaccinated-Susceptible (SIVS) model with infection reintroduction. Health policies depend on vaccine coverage, $v_0$, that guarantees herd immunity levels in the population. Vaccine failures occur when an organism develops a disease despite of being vaccinated against it. After vaccination, a proportion of healthy individuals unsuccessfully tries to increase antibody levels and, consequently these individuals are not immune to the vaccine preventable disease. When an infectious process is in progress, the initial vaccine coverage drops down and herd immunity will be lost. Our objective was to introduce a warning vaccination level and define random measures quantifying the time until the number of vaccinated descends to a warning vaccination level (i.e., the so-called sleeping period), and the epidemic size. A sensitivity analysis was performed to assess the influence of the model parameters on the variation and robustness of the sleeping period and the number of infections observed within it.

1. Introduction

Infectious diseases have a huge impact on human society across the globe and throughout history. At the moment, people all over the world are concerned about the novel corona-virus disease (COVID-19) and its dramatic clinical and nonclinical consequences. This is an example of new emerging diseases (avian influenza, Ebola, SARS, etc.) but old-time significant diseases are still present today (diphtheria, measles, polio, tuberculosis, etc.). The main purpose for public health is to increase knowledge about the spread of a disease within a community, with the goal to establish control measures and surveillance strategies to minimize or to stop viral transmission.

Mathematical models have been used to investigate the dynamics of infectious disease since the first mathematical model of Bernouilli, proposed in 1760 for smallpox inoculation [1]. In general terms, there are mainly two mathematical approaches to represent the evolution of an epidemic process: the deterministic approach whose models rely on differential equations and the stochastic approach based on models that employ stochastic processes (Markov chains, branching and diffusion processes, etc.). Historically, deterministic models have received more attention. Although the 1926 stochastic model by A.G. Mckendric [2] preceded his well-known deterministic model in the joint work with W.O. Kermack [3], the early model passed practically unnoticed. Both types of models play an important role in epidemiology showing main differences in their asymptotic behavior. According to several authors [4,5,6], stochastic modeling of epidemics is important in small populations or to model phenomenon where the epidemic outcome depends strongly on the variability in demography, disease transmission or environment.

A possible way to formulate a stochastic model is to directly use its deterministic counterpart [7]. It requires defining random variables depending on individual dynamics. In Markov chain models, the transition probabilities depend on the time between events and the process providing the state of the population at any time satisfy the Markov property. Stochastic differential equations (SDE) models are the natural extension of deterministic models, either by adding a noise term to ordinary differential equations or by randomizing in some sense relevant parameters.

Vaccination is a really effective way to prevent infectious diseases. Vaccines train the body’s immune system in a way that it can recognize and fight pathogens to which it has not been exposed before. Typically, vaccines are given to healthy people to keep them healthy. Hence, they are medications with the highest standards of safety. However a vaccine failure can occur and a proportion of vaccinated individuals do not develop immunity against its infectious pathogen [8]. The main reasons for this failure are failure of the vaccine delivery system and failure of the immune response due to inadequacies of the vaccine or host factors. People that are immune to an infectious disease can impede the transmission of a disease to vulnerable people, this is called herd protection. An interesting question is to determine the number of people in a population that should be immunized against a disease in order to practically remove disease transmission in the population. This threshold, which is called the immunity threshold, varies among diseases, depending on the reproductive potential of the pathogen. For instance, the measles virus needs a high and sustained level of vaccine coverage to interrupt disease transmission. Recent outbreaks of vaccine preventable diseases have arisen in communities with low vaccine coverage because they do not have herd protection [9].

There is a large number of publications dealing with compartmental models involving vaccination as a strategy for disease control [10,11,12,13,14]. Susceptible-Infected-Susceptible (SIS) type models are appropriate for diseases showing repeated infections, where recovery does not confer immunity. Hence, a protective vaccine alleviates, when available, many of these infections and the associated health and economic issues.

Assuming imperfect vaccination, the SIS model has been studied in a deterministic framework in [15,16], assuming nonlinear incidence in [17,18] and including a latency period as well as psychological effects in both susceptible and vaccinated individuals in [19]. A Susceptible-Infective-Vaccinated-Susceptible (SIVS) model under stochastic disturbance approach is discussed in [20,21], in a population of varying size [22], or assuming nonlinear incidence in [23].

In this paper, we deal with a stochastic compartmental model, where the propagation of a contagious disease was modeled in terms of a continuous time Markov chain (CTMC) [7,24,25,26,27]. We assumed that an infectious disease outbreak is taking place within an homogeneous moderately-sized community in which individuals are identical in terms of social mixing and disease contact or recovery characteristics. For community protection, a certain number of individuals were vaccinated against the contagious disease prior to the arrival of the outbreak.

Under the same model assumptions, in [27], the transmission potential of a contagious disease was analyzed by studying $R_{e0}$ and $R_p$, two stochastic descriptors equivalent to the basic reproduction number, $R_0$. Like $R_0$, the  above mentioned stochastic measures provide the minimal vaccination coverage needed to interrupt epidemic transmission. As the administered vaccine was not perfect, since the onset of the epidemic, the group of vaccinated people gradually disappears, leaving the whole population unprotected when the number of vaccinated declines to zero.

The purpose of the present paper was to introduce a threshold for the number of vaccinated individuals, say w, that triggers an alarm when the critical vaccination level drops down. We focused on the sleeping period, defined as the interval of time starting at the epidemic onset and ending when the warning level w is reached. More specifically, we analyzed both the length of the sleeping period and disease incidence during this time interval. To measure the strength of the model across parameter changes and to identify which one of them must be estimated accurately, we conducted a sensitivity analysis on performance measures of the random variables defined to study the length of the interval and disease incidence.

The rest of the paper is organized in the following way: Section 2 contains the description of the stochastic SIVS model with imperfect vaccine and infection reintroduction. Section 3 deals with the warning vaccination level w, the wake up time, $T_w$, and the disease incidence, $N_w$, during the sleeping period $[0,T_w]$. Theoretical and algorithmic results regarding the distribution of both random variables are provided in terms of the model parameters. Theoretical results in Section 4 come from the application of matrix calculus techniques to a perturbation analysis for performance measures of the random variables introduced in Section 3. In Section 5, we illustrate both theoretical and algorithmic results and present a local sensitivity analysis applied to human papillomavirus (HPV). Results and main insights are summarized in Section 6.

2. Model Description

The stochastic SIVS model, used to represent the evolution of a contagious disease within a population, was described in [27]. It refers to a compartmental model where individuals are classified in three groups: namely, susceptible (S), vaccinated (V), and infected (I). The model describes the movement of individuals from one compartment to another as it is shown in Figure 1.

The model tracks the changes in a constant-sized, homogeneous, and uniformly-mixed population suffering from the contagious disease. This population was not isolated. Hence, infection was transmitted due to contacts, whether inside or outside the population, with infected individuals.

To prevent disease spread, health authorities launched a vaccination policy and part of the population received a vaccine that confers immunity. Due to several reasons, the vaccine was not fully effective, a proportion of vaccinated individuals failed to increase antibody levels and was not immune to the vaccine-preventable communicable disease. Therefore, a proportion of contacts between vaccinated and infectious individuals produced new cases of infection. When this occurred, the vaccinated individual loses vaccine protection and becomes an infected individual.

After recovery, any infected individual became susceptible (i.e., unprotected) to the disease, regardless of whether or not he was previously vaccinated. There was no immunity after recovery and consequently any unprotected individual can get the infection several times during an outbreak.

Our mathematical description involved a CTMC that provides the evolution of the disease in terms of the number of individuals present at each compartment at any time t. In this framework, at time $t\ge 0$ we record the number of unprotected susceptible, $S\left(t\right)$, vaccinated, $V\left(t\right)$, and infected individuals, $I\left(t\right)$.

On the other hand, the constant population size hypothesis implies that $S\left(t\right)+V\left(t\right)+I\left(t\right)=N$. Assuming that at $t=0$ there are $v_0$ vaccinated individuals, with $0<v_0\le N$, the underlying stochastic model is a bidimensional CTMC

$$X=\left\{\right(V\left(t\right),I\left(t\right)):t\ge 0\},$$

with state space $\mathbf{S}=\{(v,i):0\le v\le v_0,0\le v+i\le N\}$, that has a cardinality of $(v_0+1)(N+1-v_0/2)$ states.

As provided in [27], the structure of the infinitesimal generator of X presents a block bidiagonal structure that is helpful for computational purposes.

Regarding the dynamics and evolution of the disease, we summarize model parameters in

Table 1. Summary of model parameters.

Size of the population	N
Rate for disease internal transmission	$\beta$
Rate for disease external transmission	$\xi$
Probability of vaccine failure	h
Rate for recovery	$\gamma$

and we describe the effective events that produce a change in the current state of the epidemic process. Therefore,

Table 2. Outgoing events and their transition rates.

Effective Outgoing Event	Transition	Rate
Susceptible-Infectious contagion	$(v,i)\to (v,i+1)$	$(\frac{\beta i}{N}+\xi )(N-v-i)$
Vaccinated-Infectious contagion	$(v,i)\to (v-1,i+1)$	$h(\frac{\beta i}{N}+\xi )v$
Recovery and loss of immunity	$(v,i)\to (v,i-1)$	$\gamma i$

enlists events jointly with outer transitions and rates. In more detail, the effective events were the following:

$E_1$. A susceptible individual gets the infection through a direct contact with an infected individual either within the population or from an external source of infection.

$E_2$. Since the vaccine was not fully effective, a vaccinated individual becomes infected through a contact with an infectious individual either within the population or from an external source of infection.

$E_3$. An infectious individual recovers and becomes susceptible.

Times spent at each state $(v,i)\in \mathbf{S}$ are independent and exponentially distributed random variables, with rates $q_{v,i}$, depending on the state and on the model parameters, which are as follows

$$q_{v,i}=(\frac{\beta i}{N}+\xi )(N-v-i)+h(\frac{\beta i}{N}+\xi )v+\gamma i.$$

(1)

In more detail, transitions out of a specific general state $(v,i)\in \mathbf{S}$ are depicted in Figure 2, where the appearing rates are introduced in order to ease the notation in the sequel. More explicitly, we define transition rates as follows

$$\begin{array}{ccc}\hfill {\gamma }_i& =& (1-{\delta }_{i,0})\gamma i,\hfill \\ \hfill {\lambda }_{v,i}& =& (1-{\delta }_{i,N-v})(\frac{\beta i}{N}+\xi )(N-v-i),\hfill \\ \hfill {\eta }_{v,i}& =& (1-{\delta }_{v,0})h(\frac{\beta i}{N}+\xi )v\hfill \end{array}$$

(2)

where the symbol ${\delta }_{i,j}$ represents the well-known Kronecker delta function, which takes the value 1 if $i=j$ and 0 otherwise.

When no vaccination takes place after $t=0$ and $\xi >0$, the long-term behavior of X is given by the stationary distribution, which is concentrated in the absorbing set of states with no vaccinated individuals $\left\{\right(0,i):0\le i\le N\}$ and agrees with the stationary distribution of the number of infectious individuals in a SIS model with external source of infection [27,28]. That is, the protection provided by the vaccine fades away almost surely in finite expected time, leaving the whole population unprotected and vulnerable to the disease.

The aim of this paper was to introduce a warning vaccination level, $w<v_0$, that alerts health authorities not only about the danger of a major outbreak but also the need to urge some individuals to get re-vaccinated.

3. The Warning Vaccination Level

In this section, we analyze the random variables $T_w$ and $N_w$ related to the warning vaccination level, w. We assume that, for an initial vaccine coverage $v_0$, w is an integer such that $0\le w\le v_0$.

At time $t=0$, $v_0$ is the number of vaccine protected individuals and the first case of infection is detected. During the outbreak, the number of vaccinated individuals $V\left(t\right)$ decreases according to the number of infections that are due to vaccine failure. As soon as the number of vaccine protected individuals drops to w an alert is triggered.

Hence, given a warning level, we first focus on the time at which the alert is activated. We refer to the random variable $T_w$ as the wake-up time and we call the random interval $[0,T_w]$ the sleeping period. Secondly, we put attention on $N_w$, the number of cases of infection among unprotected people taking place until the wake-up time, that is, during the sleeping period $[0,T_w]$.

3.1. The Wake-Up Time, $T_w$

The aim of this section is to analyze the period of time that the infectious process needs to hit the warning level for vaccination.

It is interesting to know how fast the warning level will be reached to take an appropriate action like planning a re-vaccination.

Given a warning level for vaccination w, $0\le w\le v_0$, we are interested in the random variable

$$T_w=inf\{t\ge 0:V\left(t\right)=w\},$$

which indicates how long the number of vaccinated individuals takes to reach the warning level w. Additionally, $T_w$ is called first passage time into the set $\left\{\right(w,i):0\le i\le N-w\}$. It can be shown that $T_w$ presents a phase-type distribution [29], that involves the exponential of a matrix, which is a sub-matrix of the infinitesimal generator of X. In this paper we proceed by using a different approach based on the first-step methodology.

Let us denote by $\widehat{W}=\{(v,i)\in \mathbf{S}:w\le v\le v_0,0\le i\le N-v\}$. We are going to study the probabilistic behavior of the random variable $T_w$ conditioned to the current situation, that is the set of random variables $\{\left(T_w\left|V\left(0\right)=v,I\left(0\right)=i\right.\right)$, for $(v,i)\in \widehat{W}\}$. Therefore, we need to identify Laplace transforms and central moments for these conditioned variables. To ease theoretical derivations, we eliminate the warning level from the notation of transforms and central moments. Hence, for $s\in \mathbb{C}$, with $Re\left(s\right)\ge 0$ and $(v,i)\in \widehat{W}$ we define

$$\begin{array}{ccc}\hfill {\psi }_{v,i}\left(s\right)& =& E\left[e^{-s{T}_w}\left|V\left(0\right)=v,I\left(0\right)=i\right.\right],\hfill \\ \hfill M_{v,i}^k& =& E\left[T_w^k\left|V\left(0\right)=v,I\left(0\right)=i\right.\right]={(-1)}^k\frac{{\partial }^k\left[{\psi }_{v,i}\left(s\right)\right]}{\partial s^k}{|}_{s=0}.\hfill \end{array}$$

When we condition to the initial situation of the outbreak, we can state some basic results.

$$\begin{array}{ccc}\hfill P(T_w& =& 0\left|V\left(0\right)=w,I\left(0\right)=i\right.)=1, \mathrm{for} 0\le i\le N-w.\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill P(T_w& <& \infty \left|V\left(0\right)=v,I\left(0\right)=i\right.)=1, \mathrm{for} (v,i)\in \widehat{W}.\hfill \end{array}$$

(4)

Result (3) is trivially true from the definition of $T_w$ by noticing that the warning level agrees with the initial number of vaccinated individuals. Result (4) comes from the long-term behavior of the CTMC X. On one hand, $\widehat{W}$ is a finite set of transient states. On the other hand, the mean time to absorption, in the set $\left\{\right(0,i):0\le i\le N\}$, is finite. Hence, any level of vaccination will be reached in a finite time almost surely.

Moreover, we notice that Equation (4) gives that

$${\psi }_{v,i}\left(0\right)=M_{v,i}^0=1, \mathrm{for}(v,i)\in \widehat{W}.$$

(5)

Additionally, from Equation (3), we get that

$${\psi }_{w,i}\left(s\right)=1, \mathrm{for} 0\le i\le N-w,$$

(6)

at any point $s\in \mathbb{C}$, with $Re\left(s\right)\ge 0$.

In fact, the result in Equation (6) will be the starting point to determine the remaining set of Laplace transforms for a given number of vaccinated individuals v, with $w<v\le v_0$.

Let us assume that the initial situation of the outbreak is $V\left(0\right)=v$, and $I\left(0\right)=i$, with $w+1\le v\le v_0$ and $0\le i\le N-v$. A first step argument conditioning on the exponentially distributed time to the first transition, out of the state $(v,i)$, gives the following relationship

$$\begin{array}{c}\hfill {\psi }_{v,i}\left(s\right)=\frac{{\gamma }_i}{s+q_{v,i}}{\psi }_{v,i-1}\left(s\right)+\frac{{\lambda }_{v,i}}{s+q_{v,i}}{\psi }_{v,i+1}\left(s\right)+\frac{{\eta }_{v,i}}{s+q_{v,i}}{\psi }_{v-1,i+1}\left(s\right).\end{array}$$

That is equivalent to

$$\begin{array}{c}\hfill -{\gamma }_i{\psi }_{v,i-1}\left(s\right)+(s+q_{v,i}){\psi }_{v,i}\left(s\right)-{\lambda }_{v,i}{\psi }_{v,i+1}\left(s\right)={\eta }_{v,i}{\psi }_{v-1,i+1}\left(s\right).\end{array}$$

(7)

Remark 1.

Notice that, at any point s, the right hand side term of Equation (7) just depends on model parameters and on Laplace transforms of one level less of vaccination. Therefore, it is possible to get the numerical value of any Laplace transform ${\psi }_{v,i}\left(s\right)$ in a recursive way starting from the boundary result (6), in the natural order for v. With the help of numerical methods for Laplace transforms inversion (see [30,31]) it is possible to calculate probability distribution functions. Although the numerical inversion is indeed possible, it is many times computationally not feasible. In  our model, the recursive solution of Equation (7) could be particularly useful to get computational results.

Apart from the above observation, Equation (7) is the basis to get the central moments of the wake up time random variable: $T_w$.

Firstly, we need to introduce some notation. For any number of vaccinated, v, with $w\le v\le v_0$ and any integer $k\ge 0$, we denote

$$\begin{array}{ccc}\hfill {\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}}& =& {(M_{v,0}^k,\cdots ,M_{v,N-v}^k)}^{{}^{\prime }},\hfill \\ \hfill {\tilde{\mathbf{M}}}_{\mathbf{v}}^{\mathbf{k}}& =& {(M_{v,1}^k,\cdots ,M_{v,N-v}^k)}^{{}^{\prime }}\hfill \end{array}$$

where notation ${}^{\prime }$ stands for transpose vectors.

In what follows ${\mathbf{1}}_j$ and ${\mathbf{0}}_j$ will represent $j-$dimensional all-ones and all-zeroes vectors, respectively.

The following theorem provides central moments as the solution of a system of linear equations.

Theorem 1.

Given a fixed warning vaccination level w, with $0\le w\le v_0$, and an integer k, $k\ge 0$, the central moments of order k, are the solution of the following system of linear equations:

$$\begin{array}{ccc}\hfill {\mathbf{M}}_{\mathbf{v}}^{\mathbf{0}}& =& {\mathbf{1}}_{N-v+1},forw\le v\le v_0,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\mathbf{M}}_{\mathbf{w}}^{\mathbf{k}}& =& {\mathbf{0}}_{N-w+1},fork\ge 1,\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{v}}{\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}}& =& k{\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{M}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}},fork\ge 1andw<v\le v_0.\hfill \end{array}$$

(10)

where ${\mathbf{R}}_{\mathbf{v}}$ and ${\mathbf{D}}_{\mathbf{v}}$ are $(N-v+1)$ square matrices whose non-null entries are

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{v}}(i,j)& =& \left\{\begin{array}{cc}\hfill -{\gamma }_i,& ifj=i-1and1\le i\le N-v\hfill \\ \hfill q_{v,i},& ifj=iand0\le i\le N-v\hfill \\ \hfill -{\lambda }_{v,i},& ifj=i+1and0\le i\le N-v-1\hfill \end{array}\right.\hfill \end{array}$$

and ${\mathbf{D}}_{\mathbf{v}}(i,j)={\eta }_{v,i}$, for $j=i$ and $0\le i\le N-v$.

Proof. 

First we notice that Equation (8) summarizes the result appearing in Equation (5). By differentiating repeatedly k times Equation (6) with respect to s, and setting $s=0$, we get the result appearing in Equation (9).

For any initial situation $(v,i)$, such that $w+1\le v\le v_0$ and $0\le i\le N-v$, to obtain higher order central moments, we take $k\ge 1$ derivatives on Equation (7) with respect to s and, after evaluation in $s=0$, we get

$$\begin{array}{c}\hfill -{\gamma }_i{M}_{v,i-1}^k+q_{v,i}{M}_{v,i}^k-{\lambda }_{v,i}{M}_{v,i+1}^k=k{M}_{v,i}^{k-1}+{\eta }_{v,i}{M}_{v-1,i+1}^k.\end{array}$$

(11)

This system of equations can be written in matrix form as appears in Equation (10). □

Theorem 1 provides a recursive scheme for the computation of any central moment $M_{v,i}^k$, for $k\ge 0$ and $(v,i)\in \widehat{W}$. It is summarized in the following algorithmic scheme (Algorithm 1).

Algorithm 1$T_w-$ central moments.

For any integer $k,$ $0\le k\le N-w$, central moments of order k for the random variable $T_w$ are determined from the following scheme: Step 1: Set $j=0$.      Step 1a: For $w\le v\le v_0$, set ${\mathbf{M}}_{\mathbf{v}}^{\mathbf{0}}={\mathbf{1}}_{N-v+1}$. Step 2: Set $j=1$.      Step 2a: Set $v=w$ and ${\mathbf{M}}_{\mathbf{v}}^{\mathbf{j}}={\mathbf{0}}_{N-v+1}$. If $w=v_0$ go to Step 3.      Step 2b: Set $v=v+1$. While $v\le v_0$ compute $${\mathbf{M}}_{\mathbf{v}}^{\mathbf{j}}={\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1}(j{\mathbf{M}}_{\mathbf{v}}^{\mathbf{j}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{M}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{j}}).$$ Step 3: Set $j=j+1$. While $j\le k$, go to Step 2a.

3.2. Epidemic Transmission during the Sleeping Period $[0,T_w]$: Number of Infectious Cases

In this section, we investigate the epidemic transmission until the number of individuals in the vaccinated compartment drops to the warning level, w. More precisely, we observe the number of cases of infection that take place in the time interval $[0,T_w]$. We recall that a vaccinated individual, once infected and subsequently recovered becomes susceptible to the infection. Hence, the epidemic final size in the time interval $[0,T_w]$ amounts to $(v_0-w)$ infections due to vaccine failure plus the number of infections related to susceptible individuals (i.e., unprotected people).

We introduce a random variable, $N_w$, that provides the number of infections of unprotected individuals that take place during $[0,T_w]$. As $T_w$ was analyzed conditioned on the population description at $t=0$, we study a conditional version of $N_w$, that is $\left(N_w\right|V\left(0\right)=v,I\left(0\right)=i)$, for $(v,i)\in \widehat{W}$.

First, we notice that when the warning level is chosen equal to the vaccine coverage at $t=0$ no infections of susceptible can occur, that is

$$\begin{array}{c}\hfill P(N_w=0|V\left(0\right)=w,I\left(0\right)=i)=1,\mathrm{for}0\le i\le N-w.\end{array}$$

(12)

Equation (12) shows that, for this particular initial condition, $N_w$ is a random variable of finite support; however, in general, this property is not satisfied when the warning level is fixed below the initial vaccination level. Instead of that there is an unbounded set of outcomes for the number of infections of susceptible individuals taking place in the sleeping period $[0,T_w]$. Nevertheless, from the long-term behavior of the CTMC X and the result shown in Equation (4) we get

$$\begin{array}{ccc}\hfill P(N_w& <& \infty |V\left(0\right)=v,I\left(0\right)=i)=1, \mathrm{for} (v,i)\in \widehat{W}.\hfill \end{array}$$

(13)

We recall that any outbreak of the infectious process is detected as soon as the first infection occurs. Consequently, we assume that at time $t=0$ the population is composed of $v_0$ vaccinated, a single infectious individual and $(N-v_0-1)$ susceptible.

Let us fix a warning vaccination level w, with $0\le w\le v_0$. To study the probabilistic behavior of $\left(N_w\right|V\left(0\right)=v_0,I\left(0\right)=1)$ we introduce a set of auxiliary variables where, in order to ease theoretical derivations, the warning level has been removed from the notation. Hence, for any $(v,i)\in \widehat{W}$, we denote

$$N_{v,i}\equiv \left(N_w\right|V\left(0\right)=v,I\left(0\right)=i).$$

Probabilistic behavior of any random variable $N_{v,i}$, with $(v,i)\in \widehat{W}$, is analyzed from its mass probability and generating functions, and factorial moments. Notation for probabilities, generating functions and factorial moments is as follows:

$$\begin{array}{ccc}\hfill x_{v,i}^k& =& P(N_{v,i}=k), \mathrm{for} k\ge 0,\hfill \\ \hfill {\varphi }_{v,i}\left(z\right)& =& E\left[z^{N_{v,i}}\right]=\sum _{k=0}^{\infty }{z}^{k}P[N_{v,i}=k], \mathrm{for} z\in \mathbb{C},\left|z\right|\le 1,\hfill \\ \hfill m_{v,i}^k& =& \frac{{\partial }^k\left[{\varphi }_{v,i}\left(z\right)\right]}{\partial z^k}{|}_{z=1}, \mathrm{for} k\ge 0.\hfill \end{array}$$

(14)

Observe that, due to the result shown in Equation (13), the factorial moments of order zero are

$$\begin{array}{c}\hfill m_{v,i}^0={\varphi }_{v,i}\left(1\right)=1,\end{array}$$

(15)

and we have that $m_{v,i}^k=E[N_{v,i}(N_{v,i}-1)\cdots (N_{v,i}-k+1)],$ for $k\ge 1$.

As a consequence of the Markovian property and of Equation (12) we obtain some trivial results when the number of vaccinated individuals reaches the warning level.

$$\begin{array}{ccc}\hfill x_{w,i}^k& =& {\delta }_{0,k}, \mathrm{for} 0\le i\le N-w,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\varphi }_{w,i}\left(z\right)& =& 1, \mathrm{for} 0\le i\le N-w \mathrm{and} \left|z\right|\le 1,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill m_{w,i}^k& =& {\delta }_{0,k}, \mathrm{for} 0\le i\le N-w \mathrm{and} k\ge 0.\hfill \end{array}$$

(18)

Let us begin by exploring generating functions. As we did with the analysis of $T_w$, we use a first step argument conditioning on the first effective event leading to an outer transition from a state $(v,i)\in \widehat{W}$ to get the relationship

$${\varphi }_{v,i}\left(z\right)=\frac{{\gamma }_i}{q_{v,i}}{\varphi }_{v,i-1}\left(z\right)+\frac{{\lambda }_{v,i}}{q_{v,i}}z{\varphi }_{v,i+1}+\frac{{\eta }_{v,i}}{q_{v,i}}{\varphi }_{v-1,i+1}\left(z\right).$$

That is equivalent to

$$-{\gamma }_i{\varphi }_{v,i-1}\left(z\right)+q_{v,i}{\varphi }_{v,i}\left(z\right)-{\lambda }_{v,i}z{\varphi }_{v,i+1}\left(z\right)={\eta }_{v,i}{\varphi }_{v-1,i+1}\left(z\right).$$

(19)

Equation (19) is the starting point in the derivation of factorial moments of the random variables $N_{v,i}$, for $(v,i)\in \widehat{W}$. We follow the lines stated in the previous section and, for any integer $k\ge 0$, we obtain factorial moments of order k as the solution of a system of linear equations. First, we introduce some notation for vectors involving moments $m_{v,i}^k$ that will appear in the Theorem 2 that summarizes the complete set of results for computing any factorial moment

$$\begin{array}{ccc}\hfill {\mathbf{m}}_{\mathbf{v}}^{\mathbf{k}}& =& {(m_{v,0}^k,\cdots ,m_{v,N-v}^k)}^{{}^{\prime }},\hfill \\ \hfill {\widehat{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{k}}& =& {(m_{v,1}^k,\cdots ,m_{v,N-v}^k,0)}^{{}^{\prime }},\hfill \\ \hfill {\tilde{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{k}}& =& {(m_{v,1}^k,\cdots ,m_{v,N-v}^k)}^{{}^{\prime }}.\hfill \end{array}$$

Theorem 2.

Given $k\ge 0$, factorial moments of order k in the set $\{m_{v,i}^k:(v,i)\in \widehat{W}\}$ are recursively computed from the following equations:

$$\begin{array}{ccc}\hfill {\mathbf{m}}_{\mathbf{v}}^{\mathbf{0}}& =& {\mathbf{1}}_{N-v+1},forw\le v\le v_0,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\mathbf{m}}_{\mathbf{w}}^{\mathbf{k}}& =& {\mathbf{0}}_{N-w+1},fork\ge 1,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{v}}{\mathbf{m}}_{\mathbf{v}}^{\mathbf{k}}& =& k{\mathbf{L}}_{\mathbf{v}}{\widehat{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{m}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}},fork\ge 1andw<v\le v_0,\hfill \end{array}$$

(22)

where matrices ${\mathbf{R}}_{\mathbf{v}}$ and ${\mathbf{D}}_{\mathbf{v}}$ are defined in Theorem 1 and ${\mathbf{L}}_{\mathbf{v}}$ is a diagonal matrix of dimension $N-v+1$, with non-null elements ${\lambda }_{v,i}$, for $0\le i\le N-v$.

Proof. 

We start from the explicit results shown in Equations (15) and (18), which can be expressed in matrix form as appears in Equations (20) and (21) in the statement of this theorem.

Remaining moments will arise by taking $k\ge 1$ derivatives regarding z, followed by an evaluation at $z=1$, in Equation (19). In that sense, for any initial state $(v,i)$, with $w+1\le v\le v_0$ and $0\le i\le N-v$, we get

$$\begin{array}{c}\hfill -{\gamma }_i{m}_{v,i-1}^k+q_{v,i}{m}_{v,i}^k-{\lambda }_{v,i}{m}_{v,i+1}^k=k{\lambda }_{v,i}{m}_{v,i+1}^{k-1}+{\eta }_{v,i}{m}_{v-1,i+1}^k,\end{array}$$

(23)

which can be expressed in matrix form as it is shown in Equation (22).

Notice that the right hand side of Equation (23) depends on moments of one order less and moments related to one vaccinated individual less. Hence, we can solve the system of equations given in (22) in a recursive way starting from the explicit results that appear at Equations (20) and (21). □

Remark 2.

Another algorithm can be written for the recursive computation of factorial moments. The scheme is similar to the one shown in the Algorithm 1 with the natural changes in notation and substituting computation in Step 2b by

$${\mathbf{m}}_{\mathbf{v}}^{\mathbf{j}}={\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1}(j{\mathbf{L}}_{\mathbf{v}}{\widehat{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{j}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{m}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{j}}).$$

For any initial situation $(v,i)\in \widehat{W}$, the distribution of $N_{v,i}$, the number of infections taking place in the interval $[0,T_w]$, could be computed recursively by inverting generating functions ${\varphi }_{v,i}\left(z\right)$ with the help of the Equations (17) and (19), and applying a fast Fourier transform (FFT) algorithm [32]. Instead of that, in Theorem 3, we derive a system of equations whose solution provides directly the mass probability function of $N_{v,i}$, introduced in expression (14).

As in the preceding sections, first, we set out some adequate notation. Let define a new set of vectors associated to the probability of having $k\ge 0$ infections of unprotected individuals during the sleeping interval $[0,T_w]$

$$\begin{array}{ccc}\hfill {\mathbf{x}}_{\mathbf{v}}^{\mathbf{k}}& =& {(x_{v,0}^k,\cdots ,x_{v,N-v}^k)}^{{}^{\prime }},\hfill \\ \hfill {\widehat{\mathbf{x}}}_{\mathbf{v}}^{\mathbf{k}}& =& {(x_{v,1}^k,\cdots ,x_{v,N-v}^k,0)}^{{}^{\prime }},\hfill \\ \hfill {\tilde{\mathbf{x}}}_{\mathbf{v}}^{\mathbf{k}}& =& {(x_{v,1}^k,\cdots ,x_{v,N-v}^k)}^{{}^{\prime }}.\hfill \end{array}$$

Theorem 3.

For any initial state $(v,i)\in \widehat{W}$, the distribution of the conditional random variables $N_{v,i}$ is given by the probabilities $\{x_{v,i}^k,k\ge 0\}$, which are determined from the following equations:

$$\begin{array}{ccc}\hfill {\mathbf{x}}_{\mathbf{w}}^{\mathbf{k}}& =& {\delta }_{0,k}{\mathbf{1}}_{N-w+1},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{v}}^{*}{\mathbf{x}}_{\mathbf{v}}^{\mathbf{k}}& =& (1-{\delta }_{0,k}){\mathbf{L}}_{\mathbf{v}}{\widehat{\mathbf{x}}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{x}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}},forw<v\le v_0,\hfill \end{array}$$

(25)

where matrices ${\mathbf{L}}_{\mathbf{v}}$ and ${\mathbf{D}}_{\mathbf{v}}$ are defined in Theorem 2 and ${\mathbf{R}}_{\mathbf{v}}^{*}$ is a square $(N-v+1)$ matrix whose non-null entries are

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{v}}^{*}(i,j)& =& \left\{\begin{array}{cc}\hfill -{\gamma }_i,& if1\le i\le N-vandj=i-1\hfill \\ \hfill q_{v,i},& if0\le i\le N-vandj=i.\hfill \end{array}\right.\hfill \end{array}$$

Proof. 

First, we observe that the result in Equation (24) is the matrix version of Equation (16). For a vaccination level v, such that $v\ge w$, we proceed by using the first step methodology. Again, by conditioning on the first possible transition out of the state $(v,i)$, for $0\le i\le N-v$, we are able to derive the following linear equation

$$x_{v,i}^k=\frac{{\gamma }_i}{q_{v,i}}{x}_{v,i-1}^k+(1-{\delta }_{0,k})\frac{{\lambda }_{v,i}}{q_{v,i}}{x}_{v,i+1}^{k-1}+\frac{{\eta }_{v,i}}{q_{v,i}}{x}_{v-1,i+1}^k,$$

which is equivalent to

$$-{\gamma }_i{x}_{v,i-1}+q_{v,i}{x}_{v,i}^k=(1-{\delta }_{0,k}){\lambda }_{v,i}{x}_{v,i+1}^{k-1}+{\eta }_{v,i}{x}_{v-1,i+1}^k.$$

(26)

For $0\le i\le N-v$, Equation (26) can be expressed in matrix form as appears in Equation (25). □

Recalling the result stated in Equation (13), we get that ${\sum }_{k=0}^{\infty }{x}_{v,i}^k=1$, for $(v,i)\in \widehat{W}$. For every number of cases of infection in the susceptible compartment, $k\ge 0$, the system of equations in (24) and (25) are solved recursively with the help of Algorithm 2. To compute the mass distribution function up to the mass point k, a stopping criterion should be provided to avoid too long iterative runs. In particular, numerical results shown in Section 5 come directly from Algorithm 2 until 99% of the probability mass of $N_{v,i}$ is accumulated.

Algorithm 2$N_w-$ distribution.

Given an integer $k\ge 0$, the set of conditional probabilities $\{x_{v,i}^k,(v,i)\in \widehat{W}\}$ can be determined through the following iterative scheme: Step 1: Set $v=w$.     Step 1a: Set $j=0$ and ${\mathbf{x}}_{\mathbf{w}}^{\mathbf{0}}={\mathbf{1}}_{N-w+1}$.     Step 1b: Set $j=j+1$. While $j\le k$, set ${\mathbf{x}}_{\mathbf{w}}^{\mathbf{j}}={\mathbf{0}}_{N-w+1}$. Step 2: Set $v=w+1$. If $v_0=w$, stop.     Step 2a: Set $j=0$ and compute ${\mathbf{x}}_{\mathbf{v}}^{\mathbf{0}}={\left({\mathbf{R}}_{\mathbf{v}}^{*}\right)}^{-1}{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{x}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{0}}$.     Step 2b: Set $j=j+1$. While $j\le k$, compute $${\mathbf{x}}_{\mathbf{v}}^{\mathbf{j}}={\left({\mathbf{R}}_{\mathbf{v}}^{*}\right)}^{-1}({\mathbf{L}}_{\mathbf{v}}{\widehat{\mathbf{x}}}_{\mathbf{v}}^{\mathbf{j}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{x}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{j}}).$$ Step 3: Set $v=v+1$. If $v\le v_0$, go to Step 2a.

4. Local Sensitivity Analysis

Related to the Markov chain that represents the evolution of the disease, in this section we use matrix calculus to conduct a perturbation or local sensitivity analysis of the steady-state results, presented in the Section 3. There is a large body of literature on the perturbation analysis of Markov chains [33], which involves problems related to measuring the difference between two transition matrices [34,35] and also local sensitivity, which is the problem of quantifying the effect of changes in the parameters defining transition probabilities on the behavior of the Markov chain. This second line of problems involves differentiation and it can be solved using the approach suggested by Caswell for absorbing chains [36] or the extension for non-absorbing structured Markov processes as in [37,38].

The study of sensitivities is of interest in modeling natural phenomena where the estimates of the model parameters come from a set of data (see [39] and the references therein). A local sensitivity analysis is a really helpful methodology to quantify the impact of a small variation of model parameters on the performance measures of interest and also to determine which parameter is of most influence. This methodology is especially relevant for models with many parameters, where this type of analysis is useful to disentangle the effect of each rate over these measures.

In this section, the perturbation problem is approached via derivatives and elasticities (i.e., proportional sensitivities). In that sense, the discussion that follows depends on the derivatives of conditional moments $M_{v,i}^k$ and $m_{v,i}^k$, for $k\ge 1$ and $(v,i)\in \widehat{W}$, with respect to the model parameters. For simplicity we denote $\mathit{\theta }={({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4)}^{\prime }={(\beta ,\xi ,\gamma ,h)}^{\prime }$.

Since we do not have closed form expressions for the moments of $T_w$ and $N_w$, we cannot derive analytical expressions for sensitivities. However, we take advantage of the matrix-form results, shown in Theorems 1 and 2, to derive recursive equations involving sensitivity and elasticity of the moments of the wake-up time, $T_w$, and the epidemic incidence during the sleeping period, $N_w$.

This analysis reveals how parameter variations, close to the baseline input values, influence model behavior. For instance, the sensitivity index defined as the partial derivative $\partial M_{v,i}^k/\partial {\theta }_j$ gives the rate of change of the moment $M_{v,i}^k$ in response to a change in the parameter ${\theta }_j$, while the remaining parameters hold constant. The sensitivity analysis provides a prospective study giving outcomes for changes on $M_{v,i}^k$ if the parameter ${\theta }_j$ were to change.

We deal firstly with the moments of the wake-up random variable $T_w$. By conditioning to the initial situation $(v,i)\in \widehat{W}$, a relationship for partial derivatives of $M_{v,i}^k$ can be obtained from expression (11), which is as follows:

$$\begin{array}{c}\hfill -{\gamma }_i\frac{\partial M_{v,i-1}^k}{\partial {\theta }_j}+q_{v,i}\frac{\partial M_{v,i}^k}{\partial {\theta }_j}-{\lambda }_{v,i}\frac{\partial M_{v,i+1}^k}{\partial {\theta }_j}=k\frac{\partial M_{v,i}^{k-1}}{\partial {\theta }_j}+{\eta }_{v,i}\frac{\partial M_{v-1,i+1}^k}{\partial {\theta }_j}\\ \hfill +\frac{\partial {\eta }_{v,i}}{\partial {\theta }_j}{M}_{v-1,i+1}^k+\frac{\partial {\gamma }_i}{\partial {\theta }_j}{M}_{v,i-1}^k-\frac{\partial q_{v,i}}{\partial {\theta }_j}{M}_{v,i}^k+\frac{\partial {\lambda }_{v,i}}{\partial {\theta }_j}{M}_{v,i+1}^k.\end{array}$$

(27)

From relationship (27), we can derive recursive schemes providing numerical values of $\partial M_{v,i}^k/\partial {\theta }_j$, for any parameter ${\theta }_j$ with $1\le j\le 4$. In fact, if we denote by $A_{v,i}^k\left({\theta }_j\right)=\frac{\partial M_{v,i}^k}{\partial {\theta }_j}$ and, for $k\ge 1$ and $w<v\le v_0$, define the vectors of derivatives as

$$\begin{array}{ccc}\hfill {\mathbf{A}}_{\mathbf{v}}^{\mathbf{k}}\left({\theta }_j\right)& =& {(A_{v,0}^k\left({\theta }_j\right),\cdots ,A_{v,N-v}^k\left({\theta }_j\right))}^{{}^{\prime }},\hfill \\ \hfill {\tilde{\mathbf{A}}}_{\mathbf{v}}^k\left({\theta }_j\right)& =& {(A_{v,1}^k\left({\theta }_j\right),\cdots ,A_{v,N-v}^k\left({\theta }_j\right))}^{{}^{\prime }},\hfill \end{array}$$

then the Equation (27) can be expressed in matrix form as follows, with the help of the matrices involved in Theorem 1

$$\begin{array}{ccc}\hfill {\mathbf{R}}_{\mathbf{v}}{\mathbf{A}}_{\mathbf{v}}^{\mathbf{k}}\left({\theta }_j\right)& =& k{\mathbf{A}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}\left({\theta }_j\right)+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{A}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}}\left({\theta }_j\right)+\frac{\partial {\mathbf{D}}_{\mathbf{v}}}{\partial {\theta }_j}{\tilde{M}}_{v-1}^k-\frac{\partial {\mathbf{R}}_{\mathbf{v}}}{\partial {\theta }_j}{M}_v^k,\hfill \end{array}$$

(28)

where $\frac{\partial {\mathbf{D}}_{\mathbf{v}}}{\partial {\theta }_j}$ and $\frac{\partial {\mathbf{R}}_{\mathbf{v}}}{\partial {\theta }_j}$ are the gradient matrices, respect to a scalar ${\theta }_j$, of ${\mathbf{D}}_{\mathbf{v}}$ and ${\mathbf{R}}_{\mathbf{v}}$, respectively.

However, changes in an outcome conditional moment are caused by changes in the vector of parameters $\mathit{\theta }$. Thus, sensitivity analysis in conditional moments requires more than a single derivative $\frac{\partial M_{v,i}^k}{\partial {\theta }_j}$. In fact, we have to solve the corresponding Equation (28) once for each parameter. We need an approach to differentiate matrix valued functions of vector arguments, whose usefulness is more relevant for model with a large number of parameters. Matrix calculus permits us to collect various partial derivatives with respect to a vector into vectors and matrices of derivatives, simplifying solving systems of differential equations. A revision on notation and properties of matrix calculus can be found in the paper by Magnus and Neudecker [40].

Starting from Equation (10) in Theorem 1, by applying some basic rules and properties of calculating differentials, we get that

$$\begin{array}{cc}\frac{d{\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }}\hfill & =-({\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1}{(k{\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{M}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}})}^{\prime }\otimes {\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1})\frac{dvec{\mathbf{R}}_{\mathbf{v}}}{d{\mathit{\theta }}^{\prime }}\hfill \\ & +{\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1}(k\frac{d{\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}}{d{\mathit{\theta }}^{\prime }}+{\mathbf{D}}_{\mathbf{v}}\frac{d{\tilde{\mathbf{M}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }})+({\left({\tilde{\mathbf{M}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}}\right)}^{\prime }\otimes {\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1})\frac{dvec{\mathbf{D}}_{\mathbf{v}}}{d{\mathit{\theta }}^{\prime }}.\hfill \end{array}$$

(29)

where $\frac{d{\mathbf{M}}_{\mathbf{v}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }}$ is the Jacobian matrix with entries given by $\frac{d{M}_{v,i}^k}{d{\theta }_j}$, ⊗ represents Kronecker product and $vecA$ transforms matrix $\mathbf{A}$ into a vector by stacking the columns of the matrix one underneath the other.

We observe that conditional random variables $\left(T_w\right|V\left(0\right)=v,I\left(0\right)=i)$, for $(v,i)\in \widehat{W}$, are defined as the first passage times to the set $\left\{\right(w,j):0\le j\le N-w\}$ which can be seen as an absorbing set of states and, consequently, the above random variables can be regarded as absorption times. Hence, arguments by Caswell (see [36,41]) for sensitivities and elasticities of continuous time absorbing Markov chains can be applied to get the result in Equation (29).

In addition, the structure of the expression appearing in Algorithm 1 can be exploited to derive an algorithm that provides the computation of partial derivatives of the moments $M_{v,i}^k$, with respect to the vector of parameters $\mathit{\theta }={({\theta }_1,\cdots ,{\theta }_4)}^{\prime }$. The starting point are the frontier conditions $\frac{d{\mathbf{M}}_{\mathbf{v}}^{\mathbf{0}}}{d{\mathit{\theta }}^{\prime }}={\mathbf{0}}_{(N-v+1)\times 4}$, for $w\le v\le v_0$, and $\frac{d{\mathbf{M}}_{\mathbf{w}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }}={\mathbf{0}}_{(N-w+1)\times 4}$, for $k\ge 1$.

Derivatives involving factorial moments of the conditional random variables $N_{v,i}=\left(N_w\right|V\left(0\right)=v,I\left(0\right)=i)$, for $(v,i)\in \widehat{W}$, can be obtained following a parallel reasoning. We start by taking derivatives on Equation (23), respect to any parameter ${\theta }_j$:

$$\begin{array}{c}\hfill -{\gamma }_i\frac{\partial m_{v,i-1}^k}{\partial {\theta }_j}+q_{v,i}\frac{\partial m_{v,i}^k}{\partial {\theta }_j}-{\lambda }_{v,i}\frac{\partial m_{v,i+1}^k}{\partial {\theta }_j}=k{\lambda }_{v,i}\frac{\partial m_{v,i+1}^{k-1}}{\partial {\theta }_j}+{\eta }_{v,i}\frac{\partial m_{v-1,i+1}^k}{\partial {\theta }_j}\\ \hfill +k\frac{\partial {\lambda }_{v,i}}{\partial {\theta }_j}{m}_{v,i+1}^{k-1}+\frac{\partial {\eta }_{v,i}}{\partial {\theta }_j}{m}_{v-1,i+1}^k+\frac{\partial {\gamma }_i}{\partial {\theta }_j}{m}_{v,i-1}^k-\frac{\partial q_{v,i}}{\partial {\theta }_j}{m}_{v,i}^k+\frac{\partial {\lambda }_{v,i}}{\partial {\theta }_j}{m}_{v,i+1}^k.\end{array}$$

(30)

From Equation (30), we can derive recursive schemes to get partial derivatives of factorial moments, regarding any single parameter. Moreover, a matrix calculus approach gives the following relationship between parameters and moment outcomes.

$$\begin{array}{cc}\frac{d{\mathbf{m}}_{\mathbf{v}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }}\hfill & =-({\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1}{(k{\mathbf{L}}_{\mathbf{v}}{\widehat{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}+{\mathbf{D}}_{\mathbf{v}}{\tilde{\mathbf{m}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}})}^{\prime }\otimes {\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1})\frac{dvec{\mathbf{R}}_{\mathbf{v}}}{d{\mathit{\theta }}^{\prime }}\hfill \\ & +{\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1}(k{\mathbf{L}}_{\mathbf{v}}\frac{d{\widehat{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}}{d{\mathit{\theta }}^{\prime }}+{\mathbf{D}}_{\mathbf{v}}\frac{d{\tilde{\mathbf{m}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }})\hfill \\ & +({\left(k{\widehat{\mathbf{m}}}_{\mathbf{v}}^{\mathbf{k}-\mathbf{1}}\right)}^{\prime }\otimes {\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1})\frac{dvec{\mathbf{L}}_{\mathbf{v}}}{d{\mathit{\theta }}^{\prime }}+({\left({\tilde{\mathbf{m}}}_{\mathbf{v}-\mathbf{1}}^{\mathbf{k}}\right)}^{\prime }\otimes {\left({\mathbf{R}}_{\mathbf{v}}\right)}^{-1})\frac{dvec{\mathbf{D}}_{\mathbf{v}}}{d{\mathit{\theta }}^{\prime }}.\hfill \end{array}$$

(31)

Numerical values of the partial derivatives regarding the vector of parameters $\mathit{\theta }$ will result in an iterative manner, starting from boundary conditions $\frac{d{\mathbf{m}}_{\mathbf{v}}^{\mathbf{0}}}{d{\mathit{\theta }}^{\prime }}={\mathbf{0}}_{(N-v+1)\times 4}$, for $w\le v\le v_0$, and $\frac{d{\mathbf{m}}_{\mathbf{w}}^{\mathbf{k}}}{d{\mathit{\theta }}^{\prime }}={\mathbf{0}}_{(N-w+1)\times 4}$, for $k\ge 1$.

In Section 5 we quantify the effect of the model parameters on moments of $T_w$ and $N_w$. To interpret sensitivities we need to be aware that model parameters are measured in different units. Failure probability h only takes values between 0 and 1. Contact or recovery rates have not such restrictions. Therefore, the sensitivity of a given moment to changes in an epidemic rate is difficult to compare with the sensitivity of the failure probability. To avoid this difficulty, we introduce elasticity indices that measure relative changes in moments when a parameter changes. Hence, in what follows, the elasticity index of a moment m, depending differentiably on a parameter $\theta$, is defined as $\frac{\partial m}{\partial \theta }\times \frac{\theta }{m}$.

5. Sensitivity Analysis

In this section, we present some numerical work for the sleeping time, $T_w$, and for the epidemic incidence during this period, $N_w$. Firstly, we present global sensitivity results when the parameters of the model vary over their entire range of interest. Secondly, a local sensitivity analysis is performed to check if results are robust with small changes in the parameters. For illustrative purposes we focus on results related to human papillomavirus.

We consider a population of $N=100$ individuals where an outbreak of a contagious disease is in progress. We assume that some individuals of the population, $v_0$, have received a vaccine to be protected against the disease. The outbreak started at the time at which the first infection occurs. Hence, at time $t=0$, the initial state of the CTMC X is $(v_0,1)$. Therefore, in the sequel we ease the notation and represent $\left(T_w\right|V\left(0\right)=v_0,I\left(0\right)=1)=T_w$ and $\left(N_w\right|V\left(0\right)=v_0,I\left(0\right)=1)=N_w$.

5.1. Global Sensitivity Analysis

Next, we focus on the effect of each parameter of the model on the sleeping period and on the incidence of infections during the sleeping period.

Let us deal first with the sleeping period. Our interest is to observe the influence of the parameters on the expected value and on the standard deviation of $T_w$. Numerical results for $E\left[T_w\right]$ and $\sigma \left[T_w\right]$ have been obtained by means of Algorithm 1.

In the first scenario we assume a 90% of vaccine coverage for a population of $N=100$ individuals. The warning level is fixed at $w=70$. Consequently, at the start of the outbreak the population contains 90 vaccinated, a single infected, and 9 susceptible individuals. The alarm will be triggered as soon as the number of vaccinated individuals drops to 70 vaccinated.

Figure 3 represents the expected sleeping time as a function of he internal transmission rate $\beta$ combined with recovery, external transmission rate, and vaccine failure, respectively. Shaded areas have been obtained by considering $E\left[T_w\right]\pm \sigma \left[T_w\right]$.

We notice that, as a function of the internal transmission rate, the expected value shows a decreasing behavior. This fact can be explained by noticing that an increase in the internal transmission rate $\beta$ increases contacts, contagions and, consequently, the time to reach the warning level is shortened.

Let us go deeper into each graph to analyze the behavior of $T_w$ regarding the remaining parameters. Figure 3a presents $E\left[T_w\right]$ and $\sigma \left[T_w\right]$ as a function of $\beta$ in a scenario with external contact rate $\xi =0.01$ and vaccine failure $h=0.1$. The curves correspond to different recovery rates, namely $\gamma =1.0,5.0$ and $10.0$. We observe that curves for $E\left[T_w\right]$ show a monotonic behavior. For a fixed $\beta$ the expected time increases with the recovery rate. We notice that the slope of the curves decreases for increasing recovery rates. All of these facts agree with the intuition. Higher recovery rates guaranty the existence of a group of susceptible individuals that are more vulnerable to the disease than the vaccinated ones. In consequence, long recovery rates contribute to a longer sleeping periods. Curves in Figure 3b represent the expected sleeping time in a scenario with recovery rate $\gamma =1.0$ and vaccine failure $h=0.1$. Each curve corresponds to internal contact rates $\xi =0.25,0.5$ and $1.5$. We can observe that $E\left[T_w\right]$ decreases for increasing external transmission rates, in agreement with the observation made above on the behavior of the expected time regarding the internal rate $\beta$. Figure 3c corresponds to scenarios with external transmission rate $\xi =0.01$ and recovery rate $\gamma =1.0$. Curves correspond to the choice on the probability of a vaccine failure as $h=0.05,0.1$ and $0.25$. We notice that for a fixed internal transmission rate, as could be expected, a less efficient vaccine makes a decrease in the sleeping period. The influence of the vaccine failure on the sleeping period is more noticeable for diseases with smaller transmission rates.

The following set of numerical results focuses on the influence of vaccine coverage and warning level on the sleeping period. We fix the external contact rate as $\xi =0.01$, the recovery rate as $\gamma =1.0$, and we consider a vaccine whose failure probability is $h=0.1$. Figure 4 displays expected values and standard deviation of $T_w$ for different choices on the internal transmission rate, the initial coverage and the warning level.

Figure 4a presents $E\left[T_w\right]$ and $\sigma \left[T_w\right]$ as a function of $\beta$ for three values of the initial coverage, that is, $v_0=50,70$, and 90 individuals. The warning level is fixed as $w=37$. All curves present a decreasing shape as a function of $\beta$. However, for a fixed internal transmission rate we observe that the expected length of the sleeping period increases with the initial coverage. In Figure 4b we assume an initial vaccine coverage of 90 individuals and represent mean and standard deviation of $T_w$ as a function of $\beta$, for several warning levels, namely $w=20,50$ and 80. As could be expected, under the same internal transmission, higher warning levels give shorter sleeping times. Finally in Figure 4c, the expected and standard deviation are represented as functions of the warning level, the internal contact rate is $\beta =1.15$ and we consider an initial number of vaccinated individuals $v_0=50,70$ and 90. Notice that the warning levels are bounded by $v_0$, that is $0\le w\le v_0$. According to the intuition, expected values and standard deviation decrease for increasing warning levels. When we fix an allowable warning level we observe the same behavior on $v_0$ described for Figure 4a.

Following examples correspond to the disease incidence in the group of susceptible individuals during the sleeping time, $N_w$. Numerical results dealing with probabilities come by applying Algorithm 2 and results for moments come from a modified version of Algorithm 1. As in the preceding examples, we consider a population of $N=100$ individuals, with an initial coverage of $v_0=90$, the recovery rate is $\gamma =1.0$, the external transmission rate is $\xi =0.01$, the vaccine failure is $h=0.1$ and we fixed the warning level at $w=70$ vaccinated.

In Figure 5, the distribution of the number of infections of susceptible individuals, taking place during the sleeping time, is plotted for different values of $\beta =1.15,5.0$ and 10.0. We can observe that for increasing rates of the internal transmission rate, the distribution of $N_w$ is displaced to the right.

Other experiments, not reported in this paper, show similar patterns when we consider the distribution of $N_w$ for several values either of the probability of vaccine failure or the warning level. In both cases, larger values for h or for w produce a right displacement of the distribution.

To have a better understanding of the distribution of the total number of infections during the sleeping time, we represent in Figure 6 a boxplot graph of the random variable $N_w+(v_0-w)$. We assume that the initial number of vaccine protected is $v_0=90$ and the outbreak starts with a single infected individual, $(v_0,i_0)=(90,1)$. Our purpose is to compare patterns of the epidemic, when we modify the warning level. The whiskers start at the lowest value of infectious cases, that is $v_0-w$, and extend to the 99th percentile. We observe that the dispersion of the total incidence, measured in terms of the inter-quartile range, decreases as a function of the warning level and the distribution changes from left-skewed to right-skewed as w increases. As expected, the epidemic process involves more individuals when warning is kept at low levels. Numerical experiments not reported in the paper, when the warning level is $w=0$, indicate that there is a 0.99 probability that the number of infections until the complete loss of vaccine protection fluctuates in the interval $(90,6700)$, with median at 3500 cases.

Our knowledge on the $N_w$ distribution can be applied to determine the optimal warning level to guarantee enough resources with a fixed probability. For instance, if we have at hand only 150 treatment units, from Figure 6, we observe that setting the warning level at $w=55$ vaccinated guaranties that the probability of having at least 150 infections, before the alarm was triggered, is 0.25, while an upper warning level, say $w\ge 65$, guaranties that the probability of having more than 150 infections prior to the alarm activation (i.e., shortage of resources) is less than 0.01.

5.2. Local Sensitivity Analysis Applied to the Spread of Human Papillomavirus Infection

In this section, we present a perturbation analysis on $T_w$ and $N_w$, to observe the effect of changing contact and recovery rates, and also vaccine efficacy. To this end, we derive elasticities of the mean and standard deviations of the time and incidence of infectious cases until the activation of the alarm. We particularize our results for an infection of human papillomavirus.

The human papillomavirus infection is a viral infection that is transmitted among individuals through a skin-to-skin contact. HPV is the most common sexually transmitted infection [42], in fact, at least 70% of sexually active people gets the HPV once in their lives. Most people get a genital HPV infection through direct sexual contact but, as HPV is a contact disease, it can occur without intercourse.

Although the incidence of the infection is high, most of the infections do not cause symptoms and go away spontaneously. However, the absence of symptoms makes it easy to transfer the virus among individuals. Due to the fact that HPV infects and resides in a latent state in epidermal cells, recurrences are very common. The contagiousness or transmissibility of the HPV virus can be measured by the basic reproductive number $R_0$, that is usually estimated depending on the underlying mathematical model. The $R_0$ values have been published for different HPV types and scenarios [42,43,44] and run from 1.09 to 5.6.

There are four HPV vaccines that protect against HPV types that cause most of HPV cancer. Clinical trials have shown that HPV vaccination induces lifelong protection against new HPV infections. All HPV vaccines have been found to have high efficacy, over 80%, depending on the vaccine used and also on the HPV type [45].

Some of the clinical characteristics of the HPV justify the use of a SIVS mathematical model to represent the evolution and transmission of this disease. HPV is a vaccine preventable disease with very common recurrences. However, we can not ignore that our model is probably far too simplistic to represent the complex dynamic of the HPV transmission. We include it as an application of the methodology.

As in the preceding section, we assume that outbreaks are detected as soon as the first case of infection occurs and that population size remains constant during the epidemic process. The baseline parameter values are summarized in

Table 3. Baseline parameter values.

Parameter	Value
N	100
$v_0$	70
$\beta$	$3.5$
$\xi$	$0.01$
$\gamma$	$1.25$
h	$0.1$

. The year is the unit of time, the recovery rate $\gamma$ was fixed as to keep expected recovery times around 300 days, the internal rate $\beta$ was chosen in order to have a basic reproduction number $R_0=\beta /\gamma =2.8$, the external rate is $\xi =0.01$, and we assume that vaccine efficacy is 90% (i.e.; $h=0.1$).

Additionally, we assume that vaccine coverage is $v_0=70$. We consider three scenarios regarding the warning level, w. Say low, medium, and high warning level condition. More specifically, a low warning level implies that the alert is activated when the number of protected individuals drops to the 10% of the initial coverage; medium and high warning levels are associated to 50% and 90% of the initial coverage, respectively.

Table 4. Mean and standard deviation of $T_w$ and $N_w$ for different warning level conditions.

	Low	Medium	High
$E\left[T_w\right]$	$43.96852$	$33.96278$	$18.33432$
$\sigma \left(T_w\right)$	$16.36562$	$16.23378$	$12.76035$
$E\left[N_w\right]$	$777.7604$	$226.6682$	$30.08812$
$\sigma \left(N_w\right)$	$130.1556$	$46.00042$	$12.74736$

shows values for mean and standard deviation of $T_w$ and $N_w$, under low, medium, and high warning conditions. Behavior regarding warning level variations agrees with the one stated in the preceding section.

Elasticities evaluated at the baseline parameter values appearing in the rest of this section were determined numerically by applying the theoretical results in Section 4. In

Table 5. $T_w$ and $N_w$ elasticities versus warning condition.

Warning	Variable	Parameter	Mean	Sd
Low Condition	$T_w$	$\beta$	−2.611602	+2.778346
		$\xi$	−0.477366	−0.253497
		h	−1.092157	+1.231468
		$\gamma$	+2.088969	−1.470340
	$N_w$	$\beta$	−1.891158	+15.90355
		$\xi$	−0.029536	+0.539073
		h	−1.006530	−0.939145
		$\gamma$	+0.920694	+0.911401
Medium Condition	$T_w$	$\beta$	−2.834902	+1.500176
		$\xi$	−0.613948	−0.267127
		h	−1.118875	+0.522567
		$\gamma$	+2.448851	−0.888822
	$N_w$	$\beta$	−1.676597	+10.02865
		$\xi$	−0.068708	+0.905739
		h	−1.019953	−0.890020
		$\gamma$	+0.745304	+0.705858
High Condition	$T_w$	$\beta$	−2.028866	−0.070704
		$\xi$	−0.819625	−0.385364
		h	−1.050907	−0.124000
		$\gamma$	+1.848492	+0.222181
	$N_w$	$\beta$	−1.031938	+1.412817
		$\xi$	−0.204952	+0.571329
		h	−0.997156	−0.893880
		$\gamma$	+0.236890	+0.301462

we show elasticities of the sleeping period, $T_w$, and of the number of cases of infection in unprotected individuals, $N_w$.

In interpreting each elasticity index for $E\left[T_w\right]$, in Table 5, we keep all other parameters fixed. We notice that, no matter what the warning conditions are, increasing contact rates (either internal or external) will lead to a decrease in the expected length of the sleeping period. Similarly, reducing vaccine effectivity shortens the sleeping period; however, increasing the recovery rate $\gamma$ increases $E\left[T_w\right]$.

Looking closer to the results for the sleeping period, we appreciate that the more influential parameters are the internal contact rate $\beta$ and the recovery rate $\gamma$. Changes in these parameters have opposite effects. In this sense, either decreasing $\beta$ or increasing $\gamma$ by a 1% increases the expected length of the sleeping period by around a 2%, no matter the warning condition.

The stochastic uncertainty of $T_w$, represented by $\sigma \left(T_w\right)$, is more affected by perturbations in $\beta$ or $\gamma$ than in h or $\xi$, whenever we set low or medium warning conditions.

Regardless of the warning condition, sign patterns for $E\left[N_w\right]$ agree with the observed for $E\left[T_w\right]$. Indicating that any increase in the length of the sleeping period will be accompanied by an increase in the number of cases of infection observed during this interval of time. Elasticities of the expected number of cases of infection show that changes in the internal contact rate $\beta$ or in the vaccine failure probability h are always among the most significant.

Elasticities for $\sigma \left(N_w\right)$ show that standard deviation increases in terms of transmission rates and of recovery rate and decreases for increasing values of the vaccine failure probability. It is noticeable that the volatility of $N_w$ is highly influenced by the internal transmission rate, under low or medium warning conditions. Where an increase of 1% on $\beta$ produces an increase larger than 10% on the stochastic uncertainty of $N_w$.

6. Conclusions and Future Work

We analyzed a stochastic SIVS model with imperfect vaccine and infection reintroduction. Disease transmission was modeled by a continuous time Markov chain recording, at any time t, both the number of vaccinated and infectious individuals.

Assuming that $v_0$ individuals in the population have been vaccinated prior to the onset of the outbreak, we set a level of vaccinated individuals, w, to trigger an alert when the number of vaccinated reaches‚ w. Our research involves two random variables, $T_w$ and $N_w$, associated to the warning threshold. $T_w$, or the wake-up time, gives the epoch, from the onset, at which the alarm is triggered. $N_w$ helps to measure disease incidence until this moment.

For the continuous variable $T_w$, probability distribution is described in terms of moments, while for $N_w$ we are able to give a complete probabilistic description in terms of its mass distribution function. Algorithmic schemes are provided for getting numerical results. Additionally, the influence of parameter variation on performance measures, related to both variables, was studied using both a global and local sensitivity analysis.

Results are related to outbreaks starting from the initial situation $(v_0,1)$ but algorithmic procedures provide results for higher initial infected individuals i, with $0\le i\le v_0$.

Given a warning level, w, performance measures for the length and disease incidence of the sleeping period are easily computable. Indeed, they are related to first passage times of the underlying Markov chain describing the evolution of the epidemics. Therefore, our study can be adapted to more involved compartmental models that include natural immunization [46] or a relationship structure among individuals [47]. On the other hand, since vaccine failure tends to leave the population completely unprotected against the disease, our future research will focus on re-vaccination models allowing transitions from the susceptible to the vaccinated compartments. This extension could be used to plan vaccination strategies where the warning level will play a role in the re-vaccination schedule.