Numerical study of an influenza epidemic dynamical model with diffusion

In this paper, a deterministic model is formulated in the aim of performing a thorough investigation of the transmission dynamics of influenza. The main advantage of our model compared to existing models is that it takes into account the effects of hospitalization as well as the diffusion. The proposed model consisting of a dynamical system of partial differential equations with diffusion terms is numerically solved using fast and accurate numerical techniques for partial differential equations. Furthermore, the basic reproduction number that guarantees the local stability of disease-free steady state without diffusion term is calculated. Various numerical simulation for different values of the model input parameters are finally presented in order to show the effect of the effective contact rate on the steady state of the different population compartments.

1. Introduction. Influenza (commonly known also as flu) is a viral infection, transmittable from one person to another through respiratory droplets produced as a consequence of coughing or sneezing by an infected person, or by contaminated hands and surfaces [6,24,35]. Its common symptoms include fever, cough, headache, muscle and joint pain, and reddened eyes [24,35]. In some severe case, the infection may lead to death, especially among young children, old persons, pregnant women and patients with weakened immune system or special medical conditions [2,35]. In fact, since influenza virus often changes, human immune system usually fails to recognize it quickly. Hence, a person can get infected by influenza repeatedly [25]. Despite all the new advancements in health technologies, influenza is still considered a serious threat to public health [1,6,19,31,32,35]. According to estimates, about 36, 000 deaths and 114, 000 hospitalizations in the United States alone, take place every year as a consequence of influenza [5]. According to the World Health Organization, any wave of influenza epidemic can be categorized as seasonal influenza, pandemic influenza, and zoonotic or variant influenza [36]. The most common type of influenza is the seasonal influenza, which infects humans each year [2,36,37]. Influenza pandemics (such as 1918 Spanish flu) are caused by viruses to which humans do not have any immunity, but appear among humans during the pandemic [19,36]. Zoonotic or variant influenza viruses are mostly prevalent in animals, and are transmissible to humans through direct contact with the infected animal [36].
In this paper, we will formulate a reaction-diffusion model which will depict the transmission dynamics of seasonal influenza. In fact, reaction-diffusion equations are being widely used as models for spatial effects in epidemiology [20,4]. One main characteristic of reaction-diffusion models is their ability to depict the striking behavior of a population when its growth is regulated by density-dependent mortality (not exponential) [33]. This feature allows rendering the stochastic assumptions about local movement into deterministic descriptions of global densities. Furthermore, reaction-diffusion models are spatially explicit and treat space and time as continuous. This makes them more suitable in epidemiology than some other types of spatial models such as interacting particle systems and meta-population models [9]. However, solving nonlinear reaction-diffusion equations analytically is usually not an easy task. Frequently, it is even impossible to find an analytical solution for a system of nonlinear partial differential equations. Hence, in recent years, many research efforts have been devoted to finding numerical solutions for nonlinear PDEs. In our study, the proposed dynamical model will be numerically solved using an explicit finite difference method, which is a fast and efficient numerical technique.
In the past, various deterministic epidemic models studying the transmission dynamics of influenza have been proposed [2,3,7,8,11,12,16,24,28,30]. Samsuzzoha et al. [24] studied the role of diffusion on the dynamics of influenza by incorporating diffusion in a SEIR (Suceptile, Exposed, Infected, and Recovered) model. González-Parra et al. [11] studied the potential mechanisms behind 2009 epidemic waves of A−H1N1. Whereas, Alexander et al. [2] investigated the effects of vaccination on the dynamics of influenza. Moreover, A.B. Gumel [12] studied the global dynamics of a two-strain avian influenza model. Arino et al. [3] investigated the role of vaccination and antiviral treatment at the same time. On the other hand, Sharomi et al. [28] proposed a model which incorporated an imperfect vaccine and assumed antiviral dosage to infected individuals. In addition, Chong et al. [7] studied the effect of half-saturated incidence on the dynamics of influenza. Chowell et al. [8] investigated the effects of hypothetical interventions during the Spanish flu epidemic of 1918. Stilianakis et al. [30] studied the effects of drug-resistant influenza viruses by analyzing a mathematical model. Moreover, Ruan and Wang [22] investigated the global dynamics of an epidemic model with vital dynamics and nonlinear incidence rate of saturated mass action. Nuño et al. [21] proposed a mathematical model for the evaluation of the pandemic flu preparedness plans of the United States, the United Kingdom, and the Netherlands aiming to assess the role of basic control measures, antivirals, and vaccine in curtailing pandemic influenza. Furthermore, Samsuzzoha et al. [25,26,27] studied various diffusive models aiming to describe the transmission of the influenza as an epidemic and to gain basic understanding of the influenza virus behavior, especially the H 1 N 1 virus. In the aforementioned works, the splitting method has been used to numerically solve the dynamical systems. Furthermore, we mention the work of Imran et al. [14,17] who designed a deterministic model considering the administration of antiviral both as a preventive and as a therapeutic agent. One of the main purposes of the aforementioned study was to analyze the transmission dynamics and the impact of antiviral drugs in controlling the spread of the 2009 swine influenza pandemic. Finally, a deterministic SEIR model investigating the spread of the Zika fever has been developed in [15], which takes into account the effects of both horizontal and vertical disease transmissions.
In the present manuscript, we will formulate a diffusion model which will depict the transmission dynamics of seasonal influenza. In fact, diffusion equations are being widely used as models for spatial effects in epidemiology. One main characteristic of reaction-diffusion models is their ability to depict the striking behavior of a population when its growth is regulated by density-dependent mortality (not exponential) [33]. These features allow rendering the stochastic assumptions about local movement into deterministic descriptions of global densities. Furthermore, diffusion models are spatially explicit and treat space and time as continuous. This makes them more suitable in epidemiology than some other types of spatial models such as interacting particle systems and meta-population models [9].
The purpose of our study is to investigate the dynamics of influenza by mathematically modeling a system which, in particular, accounts for the diffusion as well as the hospitalization of infected individuals. Hence, we propose a nonlinear diffusion model consisting basically of a deterministic dynamical system of partial differential equations with diffusion terms, for which we will carry out some basic analysis like steady-state analysis and local stability. After performing the stability analysis, we aim to assess the transmissibility of the infection by estimating the basic reproduction number. However, solving a nonlinear system of diffusion equations analytically is usually not an easy task; and in many cases, it is even impossible to find an analytical solution. Hence, in recent years, many research efforts have been devoted to finding numerical solutions for nonlinear PDEs and for nonlinear systems of PDEs. In our study, we opt to numerically solve our nonlinear diffusion model using an explicit finite difference method. In fact, explicit finite difference schemes are easy to implement and can lead to accurate solutions if the discretization time step and the discretization space step are chosen carefully [18]. In our case, we perform a centered difference approximation of the second partial derivatives in space, and a forward discretization of the partial derivatives in time. Such discretization leads to a consistent explicit scheme which yields an accurate numerical solution of our dynamical system, under some stability assumptions. Hence, we establish the stability condition relating the time step to the space step for our model, guaranteeing the convergence of the scheme to an accurate solution of the nonlinear diffusion system. The proposed method is simple, fast, and efficiently solve the considered model. Furthermore, numerical simulation will be performed and conclusions will be drawn about the spread of influenza with different contact rates. The effect of other model parameters such as the inhibition factors will be also studied.
This paper is organized as follows. In Section 2, we present the deterministic dynamical model proposed in this study. In Section 3, we detail a stability analysis of our model and determine the eventual stability conditions. In Section 4, we present an efficient scheme to numerically solve the dynamical system of partial differential equations. We exhibit various numerical simulations in Section 5 for different values of the model parameters and propose their interpretations, then we finally conclude with a few remarks in Section 6.

2.1.
The proposed model. The core of our study is a deterministic PDE influenza epidemic model in which the total population, denoted by N (x, t) at time t and for infection distance x around an individual, is divided into five mutually exclusive compartments. Thus, it is the sum of individual populations in each compartment which includes susceptible S(x, t), exposed E(x, t), infected I(x, t), hospitalized H(x, t), and recovered R(x, t) such that, Since there is no vertical transmission of the infection, we assume that all newborns are susceptible. The susceptible population increases at a recruitment rate Π (natural birth rate) and acquires infection following an effective contact with the infected individuals at a rate β (effective contact rate). Furthermore, the susceptible population also decreases at the natural death rate µ. We also introduce the inhibition factors a and b so that the terms 1 1 + aI and 1 1 + bH measure the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infective and hospitalized individuals respectively. In fact, the inhibition factors a and b reduce the possibility of a disease spread because the endemic equilibria disappear when these two parameters are increased [22], while η represents the fraction of inhibition effect from the hospitalized individuals. Hence, the differential equation governing the dynamics of susceptible population is given by where d 1 is a diffusions constant. The exposed population is generated after susceptible individuals acquire infection at a rate β. Population in this class diminishes when individuals get infected at the rate σ. Exposed individuals also die at natural death rate µ. Thus, the PDE governing the dynamics of susceptible population is given by The population of infected individuals, generated at rate σ, decreases when these individuals are hospitalized at rate τ , or recovered (without hospitalization) at rate γ. It also decreases at the natural death rate µ and the disease-induced death rate α. Then, the PDE governing this compartment population is given by The population of hospitalized individuals is generated, when infected individuals are hospitalized at rates τ . It diminishes when individuals recover at a rate θ, and when they die naturally or due to infection at rates µ, and δ respectively. Hence, the dynamics of hospitalized population is governed by Finally, the population of recovered individuals is generated when infected and hospitalized individuals recover at rates γ and θ respectively. It decreases when recovered individuals die naturally at rate µ. Since we are investigating the transmission dynamics of only one epidemic season, we assume that every recovered individual gains immunity for the rest of the season. The differential equation for recovered individuals is therefore Taking all the equations into consideration, the model depicting the transmission dynamics of influenza with diffusion is a system of partial differential equations written as follows In all the aforementioned equations the constants d i , i = 1, . . . , 5, are the diffusivity constants. A description of all the variables and the parameters used in the aforementioned model is presented in Tables 1 and 2.

Variable
Description Total Population of all individuals S(x, t) Population of susceptible individuals E(x, t) Population of exposed individuals Population of recovered individuals Table 1. Description of the variables of the deterministic model (3) Parameter Description β Effective contact rate Π Recruitment rate µ Natural death rate for susceptible individuals α Disease-induced death rate for infected individuals δ Disease-induced death rate for hospitalized individuals σ Infection rate for exposed individuals τ Hospitalization rate for infected individuals γ Recovery rate for infected individuals θ Recovery rate for hospitalized individuals a and b Half saturation constants (inhibition factors) η fraction of inhibition effect from the hospitalized individuals Table 2. Description of the parameters of the deterministic model (3) A special case of the model (3) is the case when d i = 0, i = 1, . . . , 5, which corresponds to the following model without diffusion The analysis and computation are carried out over the space domain [−c, c] (with c > 0), for time t ≥ 0, with the following boundary conditions at x = −c and x = c, respectively: and Furthermore, similarly to [25], we use the following three different cases of initial conditions: Well-posedness of the model. The well-posedness of the model (4) will now be investigated through the following Lemma. Proof. Adding all the equations of the model (4) gives Since Hence, the model (4) is well-posed. Thus, we will now investigate its stability.
3. Existence and stability analysis of steady states. In this section, we will perform an in-depth stability analysis of the model (4). Specifically, we will study the local stability of the disease-free equilibrium along with the existence of endemic equilibrium of our model. The disease-free equilibrium state is when there is absence of infection. In that case, the entire population will consist of susceptible individuals only; all the other population compartments will be zero.

3.1.
Local stability of disease-free equilibrium (DFE). The model (4) has a DFE given by: The local stability of the DFE is explored using the next generator operator method proposed in [10]. The non-negative matrix, F , of the new infection terms, and the M -matrix, V , of the transition terms associated with the model (2.1), are respectively given by The basic reproduction number (R 0 ) is defined as the expected number of secondary cases produced by a single infection in a completely susceptible population, and is denoted by R 0 = ρ(F V −1 ), where ρ is the spectral radius [10]. The basic reproduction number of our model is given by Hence, the following result is established [10]: The DFE of our model is locally-asymptotically stable if R 0 < 1, and unstable if R 0 > 1.
The above lemma implies that the disease can be eliminated from the population whenever R 0 can be brought to (and maintained at) a value below unity, however, when it is greater than one, the disease persists as an epidemic. Mathematically speaking, the result implies that the disease can be eliminated from the population (R 0 < 1) if the initial sizes of the subpopulations of the model (4) Denote by s(A) the spectral bound of matrix A. Let ρ : Theorem 1. If R 0 > 1 then the disease is strongly uniformly ρ-persistent: whenever ρ(x(0)) > 0. Where x(t) = (S(t), E(t), I(t), H(t), R(t)) be a solution of model (4). In this case there exists an endemic steady state.
Note that X, as well as R 5 + \ X, are positively invariant. Also, all solutions originating in X converge to E 0 as t → ∞. E 0 is asymptotically stable in X.
Hence E 0 is isolated in X.
Then, from Theorem 8.17 in [29] we have that the semiflow generated by (2.4) is uniformly weakly ρ-persistent.
From the positive invariance of D, we have that (2.4) is point dissipative. Then, according to Theorem 2.28 in [29], there exists a compact attractor of points for the model (4). This, together with uniformly weakly ρ-persistent imply (15) (see "Theorem 5.2" from [29] for more details). In this case there exists an endemic steady state [29]. 4. Numerical method for solving the dynamical model.

4.1.
Finite difference discretization. In order to numerically solve the dynamical system (3), let us consider a uniform mesh {x i , i = 0, . . . , m} over the space domain[−c, c] with m elements. Let us also consider for T > 0, a uniform mesh {t j , j = 0, . . . , n} over the time domain[0, T ] with n elements. We refer to the spatial mesh size as h x and to the time mesh size as h t . We also refer to the numerical solutions with a bar notation (i.e. the numerical solution for the susceptibles S is denoted asS, etc...). Using the standard centered difference approximation of the second partial derivative with respect to x we can write, for 1 ≤ i ≤ m − 1 and 0 ≤ j ≤ n, We also use the standard forward difference approximation of the first partial derivative with respect to t we can write, for 0 ≤ i ≤ m and 0 ≤ j ≤ n − 1, Equation (2a) will be written, for j = 0, . . . , n − 1, as while for i = 0 and i = m, applying the boundary conditions (6) using a forward approximation of (6a) and a backward approximation for (6b), we havē Using similar notations and similar approximations for the partial derivatives of E, I, H, and R, equations (2b)-(2e) will be respectively written as I0,j+1 =Ī1,j+1, andĪm,j+1 =Īm−1,j+1.

4.4.
Convergence of the scheme. It is well known that explicit finite difference numerical schemes for diffusion partial differential equations of the type C ∂ 2 u ∂x 2 + f (u) = ∂u ∂t , with C > 0, are stable under the condition C h x h 2 t < 1 2 . Similarly, the proposed discretization scheme, for the system of reaction-diffusion equations, is stable provided that [13,34] where d max is denoting the maximum value of the diffusion coefficients, i.e.
5. Numerical simulation and interpretation. In this section, we present various numerical simulations for diverse values of the model parameters. In Table 3 5 show the spatial variation of the population proportions S, E, I, H, and R using the different initial conditions (i), (ii), and (iii) in both cases of presence and absence of diffusion. The constant value of the effective contact rate β = 0.514 is considered for all the aforementioned cases, similarly to [25]. The values of η and the inhibition factors a and b are assumed to be all equal to 1, unless mentioned otherwise. Results of the simulation are graphed at different times t = 0, 5, 10, 15, and 20. Additionally, three dimensional displays combining space and time variations of all the population proportions S, E, I, H, and R are exhibited in Figures 6 and 7 for the both cases: without and with diffusion, respectively.
In Figure 1, since the initial condition (i) presents uniform values along the domain, then there the population proportions S, E, I, H, and R have also uniform values along the domain for all t. In Figure 2, using initial condition (ii), we observe a slight variation throughout the domain, however the peaks of S, E, I, H, and R are all at x = 0. Similar observation is seen in Figure 3 with initial condition (iii). In both Figures 2 and 3, since there is no diffusion, the population proportions are all centered around x = 0. Nevertheless, in Figures 4 and 5, since the diffusion case is considered, we observe that the population proportions are more spread along the x axis, though, the peak values are still at x = 0. The three dimensional surface graphs of all the population proportions S, E, I, H, and R exhibited in Figures 6  and 7 give a better understanding of the spread of the disease in both space and time. Comparing both figures, it is well observed that the peak value of the infected population is at x = 0 and t = 0, then infected population decreases with time. Furthermore, with diffusion the various population compartments distributions are more spread along the x-axis. Figures 8 and 9 show the time variation of the population proportions S, E, I, H, and R for different values of the effective contact rate β, the cases without diffusion and with diffusion, respectively. In Figure 10 the time variation of S, I, H, and R is again exhibited in a different way. To avoid redundancy, only the case without diffusion is considered. All the aforementioned results are computed using the initial condition (iii). Obviously, for small values of the effective contact rate β, the proportions of infected and hospitalized populations are small. The steady state ends up with a relatively small proportion of the population in the "Recovered" compartment R (individuals who have been infected and/or hospitalized then recovered), while the main proportion of the population remains in the "Susceptible" compartment S (did not catch the disease). However, for relatively moderate and higher values of β, an important proportion of the population ends up (at the steady state) in the "Recovered" compartment R (i.e. most of the individuals of the population caught the disease and have been infected and/or hospitalized then they have recovered). In this case, just a relatively small proportion of the population remains in the "Susceptible" compartment S. It is also observed that the exposed, infected, and hospitalized populations vanish after a reasonable amount of time, while each one of the susceptible and the recovered populations reach a non zero constant steady state value.
We also show the effect of the variation of η, as well as the inhibition factors a and b. Results are displayed in Figures 11 and 12, in addition to Figure 13. It is remarkable from Figures 11 and 12 that when η increases, the peak value of the proportion of infected population as well as the peak value of the proportion of hospitalized population is higher. This conclusion applies to both cases (without diffusion and with diffusion). Nevertheless, the effect of changing the inhibition factors a and b is not major, though, Figure 13 confirms that increasing a and b would results in reducing the possibility of a disease spread (exposed, infected, and hospitalized populations proportions decrease). In fact, increasing a and b implies the disappearing of the endemic equilibria. This result is in agreement with the conclusions reached by [22]. 6. Conclusion. In this paper, we presented an influenza epidemic model consisting of a deterministic dynamical system accounting for the diffusion as well as the hospitalization of infected individuals. A simplified analysis of the model with and without diffusion is performed. The model with diffusion has been solved using an explicit finite difference numerical method. Various numerical simulations for different values of the model parameters have been plotted and discussed.