# Disease Control in Age Structure Population

Authors: Etienne Kouokam; Jean-Daniel Zucker; Franklin Fondjo; Marc Choisy
Journal: ISRN Epidemiology
 (2013)

---

## Abstract

We combine the Leslie model and its derivatives with the classical compartmental SIRS models to build a model of transmission of infected diseases, in a population of hosts, whether opened or closed systems. We calculate the basic reproductive rateR0. Under certain conditions, when R0<1, there is a disease-free equilibrium that is locally asymptotically stable. In contrast, when R0>1, this equilibrium is unstable. Then, through an example, we show how we can define public health strategies to tackle an endemic. Finally we carry a global sensitivity analysis based on this basic reproduction rate to exhibit the most influential parameters of our model that are applied to influenza.


---

## Body

## 1. Introduction

During a certain period it was believed to have vanquished a large number of infectious diseases, until they redo their appearance, sometimes more dramatically (e.g., foot and mouth disease in the United Kingdom in 2001, highly pathogenic avian influenza pathogen in Europe in 2006) [1]. The same observation was made in Africa where, in the 70s, it was believed to have neutralized human African trypanosomiasis, long fought by the WHO until it revives later and often in the same historic foci [2]. This shows how much the prevention and control of transmission of infection diseases in the community requires constant vigilance. Many scientific disciplines have addressed this problem. Among these, the contribution of mathematical models has been very helpful. Specifically, mention may be made of introduction of the concept of R0 and its consequences. Generally, R0, the basic reproductive ratio, defined as the expected number of secondary infections that occur when one infective is introduced into a completely susceptible host population, is used to characterize the nature of the disease [3–9]. So, it is a threshold used to evaluate the conditions for starting an epidemic in a given area and for a specific disease. In general, when this threshold is less than 1, the disease will disappear. By cons, when it is greater than 1, the disease spreads permanently. Knowledge of the R0 was used to evaluate the critical vaccination coverage, that is the minimum proportion of the population that must be immunized in order for the infection to die out in the population (pc=1-1/R0). As an example, when R0=3/2, the vaccination of the third of population is enough to stop the spread of the disease [10].This threshold is even more used in most epidemiological models through compartmental representations [11–13]. In such a representation, compartments are discretized, each corresponding to all individuals in the population with specific epidemiological status (susceptible, latent, infected, removed, etc.). When the dynamics depend on time only, the balance of individuals who leave a compartment and return in different other compartments can be used in a deterministic way to model the phenomenon by a systems of ordinary differential equations (ODE) associated [12]. The spread of epidemics can also be studied in space or another parameter using a system of partial differentiation. Indeed, one can use partial differential equations (PDE) to include forms of population structure other than clinical status. These systems are studied analytically or numerically to determine the conditions under which the disease spreads in the population.However, the proposed methods have limitations because very often they do not take into account the dynamics of the disease itself or the structure of the components of the study environment, such as the age structure. Indeed the structure may be very important for understanding the disease and its control; ignoring age-structured contacts is likely to result in misinterpretation of epidemiological data and potentially costly policy missteps [14–19]. For example, in a human population, the reproduction and the survival rates depend on the structure that we have chosen: a 20-year-old woman is more likely to procreate or to survive than a 60-year-old woman. For this reason, some mathematical tools, such as Leslie matrix (in discrete time) [20, 21] and its derivatives (in continuous time) have emerged [22]. These demographic models are still widely used in population dynamics and ecology because they are well suited to the methods of counting or census of individuals and their use involves relatively simple mathematical tools. An illustration of such a model is shown in Figure 1.Figure 1
Age groups with birth.A
G
i can be considered as individuals who are i years old, assuming that nobody can survive after a maximum of n years. But in general, when people's behaviors are substantially the same depending on the problem studied, several age groups will be grouped within the same single class, taking care to have disjoint classes. In this way, we have less age groups and the problem is more simplified. In this context, fi is the rate of sexual maturation from age group AGi to age group AGi+1 and bi is the reproduction rate, assuming that there is no possible reproduction before a certain age.In this paper, we are interested in the modeling of influenza. It is a contagious disease that is caused by the influenza virus. It attacks the respiratory tract in humans (nose, throat, and lungs). The influenza virus is transmitted in most cases by droplets through the coughing and sneezing of infected persons, but it can be transmitted as well by direct contact. Once recovered, an individual acquires temporary term protection and becomes susceptible again. For these reasons, we will use the SIRS compartmental structure to model the transmission. Moreover, since influenza can be highly contagious, particularly among persons without preexisting antibodies against influenza, such as young children during the interpandemic phase influenza and anyone during a pandemic, we will use a discretized age-class structure.The population is divided into three groups: susceptible (S), infected (I), and removed (R). Each group is made of two subgroups: those who can procreate and those who cannot yet procreate. We assume that the host population can vary but birth and mortality rates are independent of age within the age classes. For this reason, we use ODEs in place of PDEs, according to [23]. We further assume that inside the same age group, the population is perfectly homogeneous in terms of transmission. This transmission can vary across age group. We compute R0 using the proposed van den Driessche method [7]. We show how to use R0 to define some public health strategies. Finally, taking into account the very high number of parameters needed to calculate the basic reproductive rate, we do a global sensitivity analysis that highlights the most important parameters.


## 2. The Model

### 2.1. Description

Our population is subdivided into three disjoint clinical statuses. Each clinical status is then subdivided into two other classes, each class corresponding to one age group. The first age group is people who can not procreate while the second is those who may procreate. The number of hosts in each disjoint class or clinical status can vary with time. A schematic illustration of such the model is described in Figure2.Figure 2
A 2 SIRS age groups model: Dynamics due to infection is represented by continuous lines, while that related to demography is shown in dotted lines.Si, Ii, and Ri correspond to susceptible, infected and removed states, respectively. Here, we have BS=(1-π)b′I2+bS2+b′′R2 and BI=πb′I2.At any time an individual can be in one of the states Susceptible, infected, or removed. Individuals are said to be susceptible if they are not currently infected but can become infected [13]. They are infected if they can transmit the disease to others. As in the earlier models [11, 15], we use the structure of “Who Acquired Infection From Whom” WAIFW square matrix β in which βij is the infection rate of an individual of age class i by another of the age class j. γ is used to represent the recovery rate. So, the infection state lasts an average period of period 1/γ units of time. When an individual recovered from his illness, he enters a state of immunity and is said removed. This immunity lasts a certain period for the hosts who becomes susceptible again and the cycle repeats. The number of hosts is denoted by Si, Ii, and Ri representing the number of susceptible, infected, and removed hosts within a clinical status of age group i. i=1 and i=2 are the values assigned to the first and second age groups, respectively. Taking into account the fact that the infection status of an individual may influence the birth rate, according to [24], we assume that birth rate varies with the epidemiological status: b for susceptible hosts, b′ for those who are infected, and b′′ for removed people. Similarly, a transition rate is used to describe the transition from age group 1 to age group 2: for sake of generality, we consider f for susceptible, f′ for infected, and f′′ for removed hosts; indeed, in animal species, worms, and snail trematode parasites of crustaceans barnacles infect both gonads and can therefore affect the rate of maturation and the birth rate [4, 5]. In addition, mortality can be applied differently depending on the status of infection and age group membership: μ1 (resp., μ2), μ1′ (resp., μ2′), and μ1′′ (resp., μ2′′) for susceptible, infected and removed hosts in age group 1 (resp. 2), respectively. Finally, we assume that it enters or leaves a constant number of susceptible hosts in the compartments S1 and S2 that we denote M1 and M2, respectively. It can be considered as a number representing a balance between input and output of the system. This model leads to the system of ordinary differential equations as follows:
(1)S˙1=-β11S1I1-β12S1I2+bS2+δR1+b′′R2+(1-π)b′I2-(f+μ1)S1+M1,S˙2=-β21S2I1-β22S2I2+fS1+δR2-μ2S2+M2,I˙1=β11S1I1+β12S1I2+πb′I2-(γ+f′+μ1′)I1,I˙2=β21S2I1+β22S2I2+f′I1-(γ+μ2′)I2,R˙1=γI1-(δ+μ1′′+f′′)R1,R˙2=γI2-(δ+μ2′′)R2+f′′R1,
where (S1(0),S2(0),I1(0),I2(0))≥0.In order to meet fairly rigorous models, taking into account the assumption that birth can be contagious, a parent inI2 can give birth to a susceptible newborn. This latter returns to I1 and should instantly join the compartment S1 that is (1-b′)I2. Other births with proportion π been considered as infected for the newborn lead to the quantity πb′I2. Births from removed people are considered noncontagious. The other quantities are given following the principles of the traditional SIRS Model [12, 13].Thus,β11S1I1+β12S1I2 (resp., β21S2I1+β22S2I2) corresponds to the susceptible of age group 1 (resp., age group 2) who become infected and (f+μ1)S1 (resp., (f+μ2)S2) to the output flows from this compartment. Similarly, (γ+f′+μ1′)I1 and (γ+μ2′)I2 correspond to output variations from I1 and I2, respectively.


### 2.2. The Analytical Expression ofR0

Our goal in this subsection is to give an analytical expression forR0. The method of calculating the R0 we use is that definition of van den Driessche and Watmough [7], which asks to reorder the system of ODEs (1) in order to show the first compartments of infected before the others. We have m=2 infected compartments: I1 and I2. A disease-free equilibrium (DFE) leads to the following system of equations:
(2)-(f+μ1)S1+bS2=-M1,fS1-μ2S2=-M2.
The determinant of this system is given by (3) and we can note that if D equals zero, there is an infinite number of DFE, which is irrelevant to us. For this reason, we assume throughout the rest of this paper that the parameters are defined so that D is not equal to zero
(3)D=μ2(f+μ1)-bf.So, a DFE has the formx0=(0,0,S01*,S02*)t, where
(4)S01*=μ2M1+bM2μ2(f+μ1)-bf,S02*=(f+μ1)M2+fM1μ2(f+μ1)-bf.
A condition for this DFE to have a reasonable epidemiological interpretation is that sign(μ2M1+bM2)=sign(μ2(f+μ1)-bf)=sign((f+μ1)M2+fM1). Assuming it is the case, the Jacobian of the DFE is J(x0)=F-V where F is the part due to new infections while V is the part coming from other transfers. We have
(5)F=(β11S01*β12S01*+πb'β21S02*β22S02*);V=(γ+f′+μ1′0-f′γ+μ2′).
If we set detV=(γ+f′+μ1′)(γ+μ2′), because γ>0, β>0 and the other epidemiological parameters are nonnegative, it follows that
(6)V-1=1detV(γ+μ2′0f′γ+f′+μ1′).
For sake of simplicity, we set
(7)u=γ+μ2′v=γ+f′+μ1′.
Then it follows that
(8)FV-1=1detV(β11uS01*+β12f′S01*+πf′b′β12vS01*+πvb′β21uS02*+β22f′S02*β22vS02*)=1detVM,
where
(9)M=(β11uS01*+β12f′S01*+πf′b′β12vS01*+πvb′β21uS02*+β22f′S02*β22vS02*).After calculating, the characteristic polynomial of the matrixM is
(10)Px(M)=x2-tr(M)x+det(M)=x2-(β11uS01*+β12f'S01*+πf′b′+β22vS02*)x+det(M),
where det(M)=(β11β22-β12β21)uvS01*S02*-β21πuvb′S02*. The discriminant for Px(M)=0 is
(11)Δ=(β11uS01*+β12f′S01*+πf′b′+β22vS02*)2-4(β11β22-β12β21)uvS01*S02*+4β21πuvb′S02*.
It is easily shown that since all the parameters are all nonnegative, Δ≥0. Therefore, the spectral radius of M is
(12)ρ(M)=12(β11uS01*+β12f′S01*+πf′b′+β22vS02*+Δ).
Finally, with u an v defined at (7), we obtain
(13)R0=12uv(β11uS01*+β12f′S01*+πf′b′+β22vS02*+Δ).
From (7), we know that 1/u and 1/v are the average duration of infection of I2 and I1, respectively. Consequently, the first part of R0 that is 1/2uv, is the geometric average of infection. The second part, β11uS01*+β12f′S01*+πf′b′+β22vS02*+Δ can be considered as the total number of potential contacts. Seen this way, we can notice that there is an apparent similarity between the formulation of the R0 obtained here and its standard formulation.As shown in [7], if R0<1 then the DFE is locally asymptotically stable, and if R0>1 then the DFE is unstable.


## 2.1. Description

Our population is subdivided into three disjoint clinical statuses. Each clinical status is then subdivided into two other classes, each class corresponding to one age group. The first age group is people who can not procreate while the second is those who may procreate. The number of hosts in each disjoint class or clinical status can vary with time. A schematic illustration of such the model is described in Figure2.Figure 2
A 2 SIRS age groups model: Dynamics due to infection is represented by continuous lines, while that related to demography is shown in dotted lines.Si, Ii, and Ri correspond to susceptible, infected and removed states, respectively. Here, we have BS=(1-π)b′I2+bS2+b′′R2 and BI=πb′I2.At any time an individual can be in one of the states Susceptible, infected, or removed. Individuals are said to be susceptible if they are not currently infected but can become infected [13]. They are infected if they can transmit the disease to others. As in the earlier models [11, 15], we use the structure of “Who Acquired Infection From Whom” WAIFW square matrix β in which βij is the infection rate of an individual of age class i by another of the age class j. γ is used to represent the recovery rate. So, the infection state lasts an average period of period 1/γ units of time. When an individual recovered from his illness, he enters a state of immunity and is said removed. This immunity lasts a certain period for the hosts who becomes susceptible again and the cycle repeats. The number of hosts is denoted by Si, Ii, and Ri representing the number of susceptible, infected, and removed hosts within a clinical status of age group i. i=1 and i=2 are the values assigned to the first and second age groups, respectively. Taking into account the fact that the infection status of an individual may influence the birth rate, according to [24], we assume that birth rate varies with the epidemiological status: b for susceptible hosts, b′ for those who are infected, and b′′ for removed people. Similarly, a transition rate is used to describe the transition from age group 1 to age group 2: for sake of generality, we consider f for susceptible, f′ for infected, and f′′ for removed hosts; indeed, in animal species, worms, and snail trematode parasites of crustaceans barnacles infect both gonads and can therefore affect the rate of maturation and the birth rate [4, 5]. In addition, mortality can be applied differently depending on the status of infection and age group membership: μ1 (resp., μ2), μ1′ (resp., μ2′), and μ1′′ (resp., μ2′′) for susceptible, infected and removed hosts in age group 1 (resp. 2), respectively. Finally, we assume that it enters or leaves a constant number of susceptible hosts in the compartments S1 and S2 that we denote M1 and M2, respectively. It can be considered as a number representing a balance between input and output of the system. This model leads to the system of ordinary differential equations as follows:
(1)S˙1=-β11S1I1-β12S1I2+bS2+δR1+b′′R2+(1-π)b′I2-(f+μ1)S1+M1,S˙2=-β21S2I1-β22S2I2+fS1+δR2-μ2S2+M2,I˙1=β11S1I1+β12S1I2+πb′I2-(γ+f′+μ1′)I1,I˙2=β21S2I1+β22S2I2+f′I1-(γ+μ2′)I2,R˙1=γI1-(δ+μ1′′+f′′)R1,R˙2=γI2-(δ+μ2′′)R2+f′′R1,
where (S1(0),S2(0),I1(0),I2(0))≥0.In order to meet fairly rigorous models, taking into account the assumption that birth can be contagious, a parent inI2 can give birth to a susceptible newborn. This latter returns to I1 and should instantly join the compartment S1 that is (1-b′)I2. Other births with proportion π been considered as infected for the newborn lead to the quantity πb′I2. Births from removed people are considered noncontagious. The other quantities are given following the principles of the traditional SIRS Model [12, 13].Thus,β11S1I1+β12S1I2 (resp., β21S2I1+β22S2I2) corresponds to the susceptible of age group 1 (resp., age group 2) who become infected and (f+μ1)S1 (resp., (f+μ2)S2) to the output flows from this compartment. Similarly, (γ+f′+μ1′)I1 and (γ+μ2′)I2 correspond to output variations from I1 and I2, respectively.


## 2.2. The Analytical Expression ofR0

Our goal in this subsection is to give an analytical expression forR0. The method of calculating the R0 we use is that definition of van den Driessche and Watmough [7], which asks to reorder the system of ODEs (1) in order to show the first compartments of infected before the others. We have m=2 infected compartments: I1 and I2. A disease-free equilibrium (DFE) leads to the following system of equations:
(2)-(f+μ1)S1+bS2=-M1,fS1-μ2S2=-M2.
The determinant of this system is given by (3) and we can note that if D equals zero, there is an infinite number of DFE, which is irrelevant to us. For this reason, we assume throughout the rest of this paper that the parameters are defined so that D is not equal to zero
(3)D=μ2(f+μ1)-bf.So, a DFE has the formx0=(0,0,S01*,S02*)t, where
(4)S01*=μ2M1+bM2μ2(f+μ1)-bf,S02*=(f+μ1)M2+fM1μ2(f+μ1)-bf.
A condition for this DFE to have a reasonable epidemiological interpretation is that sign(μ2M1+bM2)=sign(μ2(f+μ1)-bf)=sign((f+μ1)M2+fM1). Assuming it is the case, the Jacobian of the DFE is J(x0)=F-V where F is the part due to new infections while V is the part coming from other transfers. We have
(5)F=(β11S01*β12S01*+πb'β21S02*β22S02*);V=(γ+f′+μ1′0-f′γ+μ2′).
If we set detV=(γ+f′+μ1′)(γ+μ2′), because γ>0, β>0 and the other epidemiological parameters are nonnegative, it follows that
(6)V-1=1detV(γ+μ2′0f′γ+f′+μ1′).
For sake of simplicity, we set
(7)u=γ+μ2′v=γ+f′+μ1′.
Then it follows that
(8)FV-1=1detV(β11uS01*+β12f′S01*+πf′b′β12vS01*+πvb′β21uS02*+β22f′S02*β22vS02*)=1detVM,
where
(9)M=(β11uS01*+β12f′S01*+πf′b′β12vS01*+πvb′β21uS02*+β22f′S02*β22vS02*).After calculating, the characteristic polynomial of the matrixM is
(10)Px(M)=x2-tr(M)x+det(M)=x2-(β11uS01*+β12f'S01*+πf′b′+β22vS02*)x+det(M),
where det(M)=(β11β22-β12β21)uvS01*S02*-β21πuvb′S02*. The discriminant for Px(M)=0 is
(11)Δ=(β11uS01*+β12f′S01*+πf′b′+β22vS02*)2-4(β11β22-β12β21)uvS01*S02*+4β21πuvb′S02*.
It is easily shown that since all the parameters are all nonnegative, Δ≥0. Therefore, the spectral radius of M is
(12)ρ(M)=12(β11uS01*+β12f′S01*+πf′b′+β22vS02*+Δ).
Finally, with u an v defined at (7), we obtain
(13)R0=12uv(β11uS01*+β12f′S01*+πf′b′+β22vS02*+Δ).
From (7), we know that 1/u and 1/v are the average duration of infection of I2 and I1, respectively. Consequently, the first part of R0 that is 1/2uv, is the geometric average of infection. The second part, β11uS01*+β12f′S01*+πf′b′+β22vS02*+Δ can be considered as the total number of potential contacts. Seen this way, we can notice that there is an apparent similarity between the formulation of the R0 obtained here and its standard formulation.As shown in [7], if R0<1 then the DFE is locally asymptotically stable, and if R0>1 then the DFE is unstable.


## 3. Numerical Solutions

The results presented hereafter correspond to the data resumed in Table1 below. The demographic data are obtained from [25, 26]. Our time unit is the day. So, the infected (resp., removed) period lasts 21 (resp., 14) days; birthrate is considered the same for susceptible, infected, and removed hosts. We are interested here in a population where women are likely to deliver early, from about 15 years; for this reason, the age group I is made of people under 15 years old. The life expectancy is 75 years at birth, so the mortality at birth is set to 1/(75*365) and is considered the same for all those in group age 1. Taking into consideration the fact that the disease can affect the life, the mortality for infected hosts is set to doubled that of susceptible. At each time step, the balance between those who leave the system and those who return is set to 10 for susceptible in age group 1 and 50 for those in age group 2. We have S01=23*80000, S02=77*80000, I01=1, I02=0 while R01=0 and R02=0.Table 1
The parameters of the simulations.
γ
δ
b
=
b
′
=
b
′
′
f
=
f
′
=
f
′
′
π
μ
1
=
μ
1
′
′
μ
1
′
μ
2
=
μ
2
′
′
μ
2
′
M
1
M
2
1
/
21
1
/
14
(
1
/
365
)
*
16.86
%
1
/
(
15
*
365
)
0
1
/
(
75
*
365
)
2
/
(
75
*
365
)
1
/
(
60
*
365
)
2
/
(
60
*
365
)
10
50


### 3.1. CaseR0>1

Withβ11=(4/3), β12=2, β21=4, and β22=0.00000008, it follows that R0=7.0708. Figure 3 presents the results we obtained. In all the simulations we have done with R0>1, there was always an asymptotic equilibrium. We have noticed that even for only one susceptible host in the numerical simulations, the disease spreads out.Figure 3
Results of numerical solutions whenR0>1.


### 3.2. CaseR0<1

In this case, We keep the same parameters as above, except eachβij is divided by 10. It follows that R0=0.70708. Asymptotically, this case always leads to a disease-free equilibrium for the numerical model, as presented in Figure 4.Figure 4
Results of numerical solutions whenR0<1.


### 3.3. Public Health Policy

The formula (13) can further be used to explore some public health interventions. For example, with the parameters above, when R0>1, the disease spreads out as illustrated in Figure 3. If the demographic parameters are modified in an endemic situation so as to allow a positive balance of only 2 susceptible in age group 1 and 5 others in age group 2 instead of 10 and 50, respectively, R0 will switch from 7.0708 to 0.81547 that is the disease will disappear after a while. This case is presented in Figure 5.Figure 5
Numerical simulations when applying health strategies.


## 3.1. CaseR0>1

Withβ11=(4/3), β12=2, β21=4, and β22=0.00000008, it follows that R0=7.0708. Figure 3 presents the results we obtained. In all the simulations we have done with R0>1, there was always an asymptotic equilibrium. We have noticed that even for only one susceptible host in the numerical simulations, the disease spreads out.Figure 3
Results of numerical solutions whenR0>1.


## 3.2. CaseR0<1

In this case, We keep the same parameters as above, except eachβij is divided by 10. It follows that R0=0.70708. Asymptotically, this case always leads to a disease-free equilibrium for the numerical model, as presented in Figure 4.Figure 4
Results of numerical solutions whenR0<1.


## 3.3. Public Health Policy

The formula (13) can further be used to explore some public health interventions. For example, with the parameters above, when R0>1, the disease spreads out as illustrated in Figure 3. If the demographic parameters are modified in an endemic situation so as to allow a positive balance of only 2 susceptible in age group 1 and 5 others in age group 2 instead of 10 and 50, respectively, R0 will switch from 7.0708 to 0.81547 that is the disease will disappear after a while. This case is presented in Figure 5.Figure 5
Numerical simulations when applying health strategies.


## 4. Sensitivity Analysis

In this section, we carry a sensitivity analysis based onR0 in order to determine the most important paremeters.


### 4.1. Theory

We consider a generic model with independent variablesXi(14)Y=f(X1,X2,…,Xn).V is the variance of Y, Vi=V(E[Y∣Xi]) the resulting variance of Y, taken over Xi, and Vij=V(E[Y∣Xi,Xj]) the average of Vi over all possible point Xi.Sensitivity analysis studies how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input [27, 28]. The methods of sensitivity analysis can be grouped into three classes: screening methods [29], which consist of a qualitative analysis of the sensitivity of the output variable to the input variables; the methods of local analysis, which quantitatively assess the impact of a small variation around a set of input values data; methods of global sensitivity analysis [30], interested in the variability of the model output in its entire range of variation. In this paper we opt for methods of global sensitivity analysis not only because they are interested in both qualitative and quantitative aspects but they save the effort spent to effectively carry out a series of local sensitivity analyzes. Indeed, they determine the entries responsible for this variability by using sensitivity indices known as Sobol indices [31]. Three of them are widely used: the first-order, the second-order, and the total indices.The first-order index generally denoted bySi quantifies the sensitivity of the output Y to the input variable Xi that is the proportion of variance in Y due to the variable Xi, while the second-order index Sij expresses the sensitivity of the variance of Y to the interaction of 2 variables Xi and Xj. When the number p of input variables is too large, the number of sensitivity indices explodes (≈2p-1). The estimation and interpretation of all these indices become almost impossible. In that case, it is convenient to use the total sensitivity index denoted by STi. The mathematical formulation of these indices is given in (15) where Ki represents the set of indices containing the index i:
(15)Si=ViV,Sij=VijV,STi=∑k∈Ki‍Sk.
The higher the index will be (close to 1), the more important the variable will be.


## 4.1. Theory

We consider a generic model with independent variablesXi(14)Y=f(X1,X2,…,Xn).V is the variance of Y, Vi=V(E[Y∣Xi]) the resulting variance of Y, taken over Xi, and Vij=V(E[Y∣Xi,Xj]) the average of Vi over all possible point Xi.Sensitivity analysis studies how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input [27, 28]. The methods of sensitivity analysis can be grouped into three classes: screening methods [29], which consist of a qualitative analysis of the sensitivity of the output variable to the input variables; the methods of local analysis, which quantitatively assess the impact of a small variation around a set of input values data; methods of global sensitivity analysis [30], interested in the variability of the model output in its entire range of variation. In this paper we opt for methods of global sensitivity analysis not only because they are interested in both qualitative and quantitative aspects but they save the effort spent to effectively carry out a series of local sensitivity analyzes. Indeed, they determine the entries responsible for this variability by using sensitivity indices known as Sobol indices [31]. Three of them are widely used: the first-order, the second-order, and the total indices.The first-order index generally denoted bySi quantifies the sensitivity of the output Y to the input variable Xi that is the proportion of variance in Y due to the variable Xi, while the second-order index Sij expresses the sensitivity of the variance of Y to the interaction of 2 variables Xi and Xj. When the number p of input variables is too large, the number of sensitivity indices explodes (≈2p-1). The estimation and interpretation of all these indices become almost impossible. In that case, it is convenient to use the total sensitivity index denoted by STi. The mathematical formulation of these indices is given in (15) where Ki represents the set of indices containing the index i:
(15)Si=ViV,Sij=VijV,STi=∑k∈Ki‍Sk.
The higher the index will be (close to 1), the more important the variable will be.


## 5. Application

Formula (13) shows that the model output can be represented only by a simple function given by R0. Therefore, the model output sensitivity is reduced to the sensitivity of this function relative to variations of 16 model parameters among 21. Indeed, 5 parameters f′′, b′′, μ1′′, μ2′′, and δ which are the output of compartments R1 and R2, are not involved in the formula R0. This situation is analogous to the classical SIRS model, where the output parameters of the compartment R which are not input to S, are not involved in the formulation of the basic reproductive ratio R0. For this reason, we can write:
(16)R0=f(β11,β12,β21,β22,γ,f,f′,b,b′,μ1,μ2,μ1′,μ2′,π,M1,M2).
To carry out the sensitivity analysis, we use the Monte-Carlo based numerical procedure for computing the full set of first-order and total-effects indices for a model of 16 factors. The method that is used here is due to Saltelli [32]. Each parameter has a standard uniform distribution U(0,1) except M1 and M2 that follow the uniform distribution U(-50,50). For the Saltelli model with 16 parameters, 100 000 perturbations have been generated.The first-order sensitivity indices (Si) for all 16 parameters are calculated numerically. It is obvious from Figures 6 and 7 that the single indices are very small as comparing to the total indices. Indeed, the sum of first-order effects is approximately 2.37×10-3, while the sum of the total indices is 4.55; as these two sums are both different from 1, this implies high importance of parameter interactions. Moreover, given that all the parameters have total indices greater than their first orders, we conclude that all of them are taking part in interactions.Figure 6
Single sensitivity indices forR0.Figure 7
Total sensitivity indices forR0.The pictures given in Figure8 show the values of R0 across the 100 000 perturbations that have been generated. In these figures, an attempt at clustering shows that the infection spreads in approximately 99% of cases while only 1% corresponds to cases where it recedes. This may explain the fact that flu rife in an endemic scale. There are almost 42% that correspond to the prohibited cases in the model above, that is, the cases where R0<0. In addition, the total indices give more information about model sensitivity and factors determining it. From Figure 7, we can distinguish three groups of parameters. As the maturation rate (f), birth rate of susceptible (b), and the mortality rates of susceptible individuals (μ1 and μ1′) retain the leading role with the total index of about 1, they are the most important parameters. The total indices for each of these parameters belongs to [0.95,1.03]. It is important to note that these four parameters are the only ones involved in formula (3) of the determinant at the DFE. Another group is made of the migration parameters (M1 and M2) which are the second most influential parameters with total indices up to 0.75. However, all the 10 remaining parameters are less important but still influential, with total indices close less than 0.22 for each of them.The variability ofR0.
(a)
(b)


## 6. Discussion

In this paper, we have coupled models that are well known and widely used not only in population dynamics but also in the transmission of infected diseases, namely, the Leslie model and its derivatives with the SIRS compartmental model. Applying our model to influenza, under certain assumptions, we have deduced an analytical expression of the basic reproductive rateR0 for which the DFE is locally asymptotically stable when it is below 1, but unstable when it is greater than 1. From an example, we have shown how to define public health policies to influence or even eliminate a disease that spreads out. This last result has already been highlighted in [33] in the context of a vector-born transmission disease. In addition, we have made a global sensitivity analysis, allowing us to accurately determine the parameters most influencing the asymptotic dynamics represented by R0. Finally, we were able to note that the numerical simulations confirm the analytical results we have exhibited in this paper. Another track will be to complete our model in order to extent it to more age groups and other compartmental models, such as SEIR or SEIRS [13, 34, 35]. Finally the effect of spatialization in such models can be investigated as it has been done in [33, 34].



---

*Source: 703230-2012-11-26.xml*