Dynamical analysis, optimal control and spatial pattern in an influenza model with adaptive immunity in two stratified population

1.   Introduction

Disease epidemics are the outcomes that occur when microorganism (bacterium, germ, plague etc.) interact with healthy organism and give rise to infection by drastically altering normal activity of bodies or encouraging immune system to build defensive comeback that results in disease symptoms. Epidemics are major considerable causes of death globally. Infections are hard to wipe out and some times leads into pandemic. Influenza is a profoundly irresistible respiratory infection that can cause serious health issues. Influenza is classified into four genera (type A, B, C and D), among which Influenza A(IAV) can infect wide range of animals, birds and humans. These situations can occur in those places which are overpopulated such as specific slum areas or those areas where the environmental conditions act as a catalyst in spreading the infection. Influenza A is one of the serious epidemic infecting which has given rise to several pandemics like swine flu. It specially affects people of different ages and the symptoms remain in a body for 6–8 days. As mentioned earlier also, growth of this disease is not very tough and there are several reasons responsible for the spread of it. For example travelling, coming in contact with the infected person, air, coughing in front of someone, etc. Breathing of respiratory droplets (which contain virus particles) from infected human usually commence influenza infection ^[1]. As per the recent data ^[2], urban area people are more affected by this disease in comparison to rural areas. One of the reasons behind it is the immunity of people staying in these regions due to which the spread of infection in both the regions follow a different pattern ^[3]. Also, this uncommon spread is due to demographic regions, populace portability, financial levels and clinical administrations between these metropolitan and provincial regions which might impact flu scourges, seriousness furthermore, transmission. According to Naba Kumar Goswami ^[4], influenza has infected around 5 million individuals and approx. 2–5 million deaths are guessed globally. Four pandemics have occurred in the last century as a result of the development of a novel influenza strain for which human population has little or no immunity. These pandemics are (i) Spanish flu resulted $40-100$ million deaths in $1918-1920$ due to H1N1 strain (ii) Asian flu resulted $1-2$ million deaths in $1957-1958$ due to H2N2 strain (iii) Hong Kong flu $0.5-2$ million deaths in $1968-1970$ due to H3N2 strain (iv) Swine Flu resulted upto $575,000$ million deaths in $2009-2010$ due to H1N1 strain ^[5,6,7]. Several researchers have proposed immunological mathematical models for influenza transmission dynamics ^[8,9,10,11]. As per ^[12] during contamination, innate immunity resistance gives principal line of protection and triggers provocative reactions. Adaptive immunity additionally assumes a basic part in the leeway of viral microbes during the later phases of contamination and it is different in different age group ^[13]. The modelling of immune in influenza is further essential because it can help in prediction of the fact that during the immune response to viral infection how antiviral therapy would be more effective. In ^{[14,15,16,17,18]}, it has been shown that how the antiviral therapy is effective to attack pathogen in presence of immune response to Influenza. Their model advised that prolonged pathogen restricts immune cells production. So, it is suggested to give treatment within two days of infection. Vaccination is also a fundamental technique to thwart contamination and reduce death rate ^[19]. The World Health Organization delivered a proposal for the arrangement of occasional flu immunization every year based on strains that are anticipated to course ^[20,21]. Despite the fact that vaccination brings defensive safe reactions against the surface antigens of flu, constant hereditary transformations permit the infection to ultimately sidestep antibody initiated assurance. The optimal control hypothesis is a fantastic scientific tool for making decisions in complex dynamical contexts. It's also been linked to infectious disease issues and it's a good strategy for figuring out how well we can control a disease ^{[22,23,24,25]}. Another significant active research which has shown very less light is the study of spatial models via spatial patterns. Because of the localised transmission or another form of interaction, many critical epidemiological and immunological phenomena are strongly influenced by space. As a result, mathematical models with time and space help to investigate the process of disease spreading. In 1952, Alan Turing was the first researcher who proposed the idea of pattern formation ^[26]. According to him, a couple of reaction-diffusion equations can be used to describe spatial patterns. In the system, a system of an ordinary differential equation is asymptotically stable without diffusion, but with the addition of diffusion, the system loses its stability and becomes unstable, which gives rise to the condition of Turing instability. This Turing instability is responsible for spatial patterns, and this idea was firstly investigated in morphogenesis and then used in other branches of science. These spatial patterns can identify the exact distribution in both space and time of the population. The spatial patterns are of two types Turing and Non-Turing. The stationary solution of the reaction-diffusion system gives rise to Turing patterns such as spots, stripes, labyrinthine, and mixed. Non-Turing patterns occur when temporal oscillation produced by Hopf-bifurcation collide with spatio-temporal patterns. The non-Turing patterns contain travelling waves, spiral, and target patterns. Segal and Jackson firstly applied the Turing instability idea to the prey-predator system for nonhomogeneous prey-predator interaction ^{[27,28,29,30,31]}. Hence, in reference to the above literature, we have developed two models in context to rural and urban areas focusing exclusively on the immune response of people in these stratified population. We have assumed that people staying in urban area have a weak immune response ^[3] so that the transmission rate of infection follow Holling type-I functional response in comparison to rural areas the infection rate follow Holling type-II functional response. We will try to conduct the comparative analysis for both the models which would include the dynamical analysis, optimal control and spatial dynamics.

1.1.   Biological process

The attachment of a virion to a target cell (a susceptible cell) is the first step in viral replication for influenza. Upon passage, the viral envelope wires with the endosomal film, empowering the viral ribonucleoproteins to enter the core of the cell, synthesise viral RNAs, and the synthesised RNA is translated into proteins. At long last, the proteins are collected into new virions at the plasma film and the new virions are let out of the cell. Innate signalling and adaptive immune responses both play important roles in protecting against influenza virus infections and achieving viral clearance. Innate immunity assumes a basic part in proficient and fast constraint of viral contamination and for adaptive immunity initiation. T-cells and antibodies (B-cells) are two adaptive immune responses that are certain to all invaders. T cells play a key role for reduction of virion. When someone is infected with Influenza A, $CD8^+$ cells are activated to attack the infected cells. It differentiates into CTLs that target the virus by producing cytotoxic granules that are responsible for restricting virus replication. $CD4^+$ cells contribute towards B cells activation and antibody production. The morbidity is reduced by B cells and it also helps in recovery upon infection. The antibodies generated by B cells make it easier to eliminate infection and also speed up the expansion of memory $CD8^+$ cells after infection. Antibody can inhibit the pathogen causing Influenza and is more important to inhibit transmission ^[32].

1.2.   Paper structure

The paper is organised as follows : In section $2$ we present our mathematical models with adaptive immune response (T cell and B cell) in context to urban and rural areas and check positivity and boundedness of the system. In section $3$, the summarized results of reproduction number, equilibrium points and local stability is done and detailed part is shown in Appendix. In subsection $3.1$ global stability is done by graph theoretic approach. Section $4$ covers the optimal control problem using two controls. Finally, section $5$ encapsulate the development and analysis of a spatial linear and non-linear model to visualize the turing effect. Further, in the numerical section $6$, all the solutions are validated and spatial patterns are obtained in validation to section $5$ followed by a detailed summarizing of our results in the conclusion discussed in section $7$.

2.   Formulation of model

We have formulated two immune response models. In model $1$, we have assumed Holling type-I functional response for transmission rate of infection whereas Holling type-II functional response is assumed for model $2$. As mentioned earlier, the idea behind suggesting two models with different functional response is due to the fact that the rate of healthy cells getting infected in rural area is less in comparison to urban area due to the fact of direct correlation between immunity and infection rate of population ^[33], further immune response of ruralite is better then metropolis ^[2,3], which is also evident from Figures 1 and 2 ^[2].

(Model 1: Urban Area) The model along with Holling type-I functional response takes the following form:

(Model 2: Rural Area) The proposed model is defined into following five classes that are density of distinct populations in time $t$.

The description of these five classes and the parameters used are mentioned in Table 1.

The flow chart showing the interaction between each of the compartment is shown in the Figure 3.

2.1.   Positivity & boundedness of the system

In this section we would be discussing the positive invariance and boundedness of the system (2.1) $\forall t > 0$.

Theorem 2.1. The system (2.1) is positively invariant.

Proof. Let $B(0) = B_0\in \mathbb R^5_+$, where $B = (X, Y, V, W, Z) \in \mathbb R^5_+$ than the system (2.1) can be written in matrix form as $\dot B = E(B)$ given as

where, $E : D_+ \rightarrow \mathbb R^5$ and $E \in D^\infty(\mathbb R^5)$.

Now, we observed that $E_i(B)|_{B_{i} = 0} = \lambda > 0$ for $(i = 1, 2, 3, 4, 5)$ where $B(0)\in \mathbb R^5_+$ as a result $B_{i} = \lambda$. Thus, the system (2.1) is positively invariant in $\mathbb R^5_+$ $\forall t > 0$ ^[34].

Theorem 2.2. $E = ((X, Y, V, W, Z) \in \mathbb R^{5}_+:0 < X+Y < \frac{\lambda}{\mu}, 0 < V < k\frac{\lambda}{\mu^2}, 0 < W+Z < \frac{\theta_1}{\theta_2})$ is the bounded region for the system.

Proof. From first two equations of system (2.1), we get

which implies $\frac{d(X+Y)}{dt}\leq \frac{\lambda}{\mu}$ by theorem of inequality ^[35].

On putting the upper bound in third equation we get $\frac{dV}{dt}\leq k \frac{\lambda}{\mu}-\mu V$, resulting in $\frac{dV}{dt} \leq \frac{k \lambda}{\mu^2}$ by theorem of inequality ^[35].

Let us assume $\theta_1 = max(gk\frac{\lambda}{\mu^2}, c\frac{\lambda}{\mu}), \theta_2 = min(h, b)$, than

which implies $\frac{d(W+Z)}{dt} \leq \frac{\theta_1}{\theta_2}$ by theorem of inequality ^[35].

This proves the boundedness of (2.1). In similar manner, the boundedness and positivity for non-linear model can be proved.

3.   Dynamical analysis

The analytic results which involves basic reproduction number, local stability for both the systems is summarized in Table 2 and the detailed calculations are mentioned in the Appendix.

3.1.   Global stability

This section would focus on the global stability of the non-trivial equilibrium point for both the Models $1$ and $2$ using graph-theoretic method as in ^[36,37] to construct the Lyapunov function we shall use a directed graph. A directed graph has a set of ordered pair say $(i, j)$ and vertices where $(i, j)$ is known as arc to terminal vertex $j$ from initial vertex $i$. For the terminal vertex $j$, $d^{-}(j)$ is the in-degree of $j$ which denotes the number of arcs in the digraph and for initial vertex $i$, $d^{+}(i)$ is the out-degree of vertex $i$ which denotes the number of arcs in the digraph. Let us consider a weighted directed graph say $\chi(N)$ over a $m\times n$ weighted matrix $N$, where the weights ($a_{ij}$) of each arc if they exist are $a_{ij} > 0$ otherwise $a_{ij} = 0$. We consider $c_{i}$ as the co-factor of $l_{ij}$ of the Laplacian of $\chi(P)$ which is given by:

If there is a strongly connected path i.e., directed to and fro path for the arcs in $\chi(P)$ then $c_i > 0 \forall i = 1, 2..., q$. We shall also use Theorem 3.3 and Theorem 3.4 from ^[36,38], which will help us in the construction of Lyapunov function. The theorems states:

Theorem 3.1.    ● If $a_{ij} > 0$ and $d^{-}(i) = 1$, for some $i, j$, then

● If $a_{ij} > 0$ and $d^{+}(j) = 1$, for some $i, j$, then

We shall also use the below theorem of ^[36].

Theorem 3.2. Let us consider an open set $L\subset R^{m}$ and a function $f: L\rightarrow R^{m}$ for a system

and assuming:

a) $\exists$ $A_{i}:L\rightarrow R$, $G_{ij}:L\rightarrow R$ and $a_{ij}\geq 0$ such that

b) For $N = [a_{ij}]$, of $(H, N)$ each directed cycle $D_{c}$ satisfies:

where $\epsilon(D_{c})$ is set of arcs in $D_{c}$

Then for $c_{i}\geq0, i = 1, \dots, m$ the function is:

satisfies $A'\vert_{3.1}\leq0$, that is, $A(z)$ is a Lyapunov function for 3.1.

Proof. Construction of Lyapunov function : Let us assume that $A_{1} = \frac{(X-X^{*})^2}{2}$, $A_{2} = Y-Y^{*}-Y^{*}\ln\frac{Y}{Y^{*}}$, $A_{3} = \frac{(V-V^{*})^2}{2}$, $A_{4} = \frac{(W-W^{*})^2}{2}$, and $A_{5} = Z-Z^{*}-Z^{*}\ln\frac{Z}{Z^{*}}$. Now by differentiation, we get $A'_{1} = (X-X^{*})X'\leq (\lambda+\mu X^*)(X+Y)+ \beta X V X^* = a_{12}H_{12}+a_{13}H_{13}$, where $a_{12} = (\lambda+\mu X^*), a_{13} = X^*$. $A'_{2} = \frac{Y-Y^*}{Y}Y' \leq \beta X V V^{*}+(\mu+pZ)(1+Y)Y^* = a_{31}H_{31}+a_{25}H_{25}$, where $a_{31} = \beta V^*, a_{23} = \beta Y^*$. $A'_{3} = (V-V^{*})V'\leq kVY+(\mu+qW)V V^* = a_{32}H_{32}+a_{34}H_{34}$, where $a_{32} = k, a_{34} = V^*$. $A'_{4} = (W-W^{*})W'\leq g(W-W^*)(V-V^*) = a_{43}H_{43}$, where $a_{43} = g$ and $A'_{5} = (\frac{Z-Z^*}{Z}) Z' \leq c(Z-Z^*)(Y-Y^*) = a_{52}H_{52}$, where $a_{52} = c$. Thus, we get an associated weighted directed graph as shown in Figure 4. Then by Theorem 3.5 ^[36] $\exists c_{i}'s, 1\leq i \leq 5$ such that $A = \sum^{q}_{i = 1}c_{i}A_{i}$ is a Lyapunov function. Using Theorem 3.3 and 3.4 we get the relation between $c_{i}$. For $a_{14} > 0$ and $d^{+}(4) = 1$, we get $c_{4}a_{43} = c_{3}a_{34}$ and for $a_{25} > 0$, and $d^{-}(2) = 1$, we get $c_{2}a_{25} = c_{1}a_{12}+c_{3}a_{32}+c_5a_{52}$. Hence, $c_{1} = c_{3} = c_{5} = 1$, $c_{4} = \frac{V^*}{g}$ and $c_{2} = \frac{\lambda+\mu X^*+k+c}{Y^*}$. Thus, the Lyapunov function is $A = A_{1}+\frac{\lambda+\mu X^*+k+c}{Y^*}+A_{3}+\frac{V^*}{g}+A_5$ and for $A'$ : $A' = (X-X^{*})X'+\frac{\lambda+\mu X^*+k+c}{Y^*}(\frac{Y-Y^{*}}{Y})Y'+(V -V^{*})V'+\frac{V^*}{g}(W -W^{*})W'+(Z -Z^{*})Z'$. Global stability for non-linear model follow the same line.

4.   Optimal control

This section would focus on the optimal control problem. Our main objective is to control the growth of infection, virus and interaction with new infected cells. The objective is to limit the mass of clinically unhealthy cells over a finite time interval $[0, T]$ at a minimal cost of efforts during influenza epidemic outbreak. For both the models, we have put two controls $T_1, T_2$, which gives the controlled equations :

$T_1$ = Efficacy of Vaccination to restrict new infection, $T_2$ = The effectiveness of antidote treatment to prevent production of virus. $(1-T_1)$ is the infection rate during vaccination and $(1-T_2)$ is the virion production rate under antidote treatment.

The aim of optimization is to maximize the following objective function

where, ${t_e}$ is the time period of treatment and the positive constants $B_1$ and $B_2$ are the benefits and costs of the introduced treatments. The controls $T_1(t)$ and $T_2(t)$ are assumed as bounded and Lebesgue integrable.

controls set is defined below as

Now we move further for Existence Theorem. For performing existence of control system we use result from Lukes^[39,40].

Theorem 4.1. There exist a function used for control pair as $(T_1^*, T_2^*)$ such that

Proof. For the proof following properties are needed ^[41]:

(1) Control set and state variable set correlating with it are not empty.

(2) $T$ is convex and closed.

(3) RHS of state system is bounded by linear function in the state and control variable.

(4) Integrand of objective functional is concave on $T$.

(5) There exist constants $m_1, m_2 > 0$ and $f > 1$ such that the integrand $I(X, W, Z, T_1, T_2)$ of the objective function satisfies

To prove these requirements, we have used Lukes'findings ^[39] to provide the existence of solutions of system (4.1), which produces condition $(1)$. By definition, the control set is convex and closed, resulting in condition $(2)$. As our system is bilinear in $T_1$ and $T_2$, the RHS of the system is also bilinear. Using the boundedness of the solutions, (4.1) satisfies condition $(3)$. Our objective functional's integrand is concave, so we also have the condition that is now required

where $m_2$ is determined by the upper bound on $X, W$ and $Z$ because $B_1 > 0$ and $B_2 > 0$. We conclude that there is an optimal control pair $(T_1^*, T_2^*)\in T$ such that $J(T_1^*, T_2^*) = max(J(T_1, T_2)), (T_1, T_2)\in T.$

The maximum principle proposed by Pontryagin ^[42] establishes the necessary conditions for an optimal control problem.

The Hamiltonian is defined as

where $f_1 = \lambda-\mu X-(1-T_1)\frac{\beta XV}{X+V}$, $f_2 = (1-T_1)\frac{\beta XV}{X+V} -\mu Y-pYZ$, $f_3 = (1-T_2)kY - \mu V - qVW$, $f_4 = gVW-hW$, $f_5 = cYZ-bZ$.

For further solution Pontryagin's maximum principle is applied and we have the following theorem.

Theorem 4.2. With the optimal condition $T = (T_1, T_2)$, solutions $G = (X^*, Y^*, V^*, W^*, Z^*)$ for corresponding state system (4.1), there exist adjoint variable $\chi_{1}, \chi_{2}, \chi_{3}, \chi_{4}, \chi_{5}$, which satisfies

with transversality conditions $\chi_{j} = 0, \; j = 1, ...., 5.$ Further more, the best control for optimality is given as $T_1^* = min(1, max(0, \frac{1}{B_1}(\chi_{2}-\chi_{1})\frac{\beta X^*V^*}{V^*+X^*})), T_2^* = min(1, max(0, \frac{1}{B_2}(kY^*\chi_{3})))$

Proof. With the help of Pontryagin's Maximum Principal the above equations can be acquired as

For solving $T_1^{*}$, $T_2^*$, optimality conditions used are $\frac{\partial H}{\partial T_1} = 0, \frac{\partial H}{\partial T_2} = 0$.

Thus, $\frac{\partial H}{\partial T_1} = B_1T_1+(\chi_{1}-\chi_{2})\frac{\beta X^*V^*}{V^*+X^*} = 0$ and $\frac{\partial H}{\partial T_2} = B_2T_2- kY^*\chi_{3} = 0$. Hence by $T$ bound fact $T_1, T_2$ can be obtained. Further substitution of $T_1, T_2$ in defined control system (4.1) following optimality system is acquired.

$\begin{cases} \frac{dX^*}{dt} = \lambda- \mu X^* - (1-T_1^*)\frac{\beta X^*V^*}{V^*+X^*}\\ \frac{dY^*}{dt} = (1-T_1^*)\frac{\beta X^*V^*}{V^*+X^*} - \mu Y^* - pY^*Z^*\\ \frac{dV^*}{dt} = (1-T_2^*) kY^* - \mu V^* - qV^*W^*\\ \frac{dW^*}{dt} = gV^*W^*-hW^*\\ \frac{dZ^*}{dt} = cY^*Z^*-bZ^*\\ < /p > < p > \chi'_{1} = 1+\chi_1 \mu+(\chi_1-\chi_2)(1-T_1^*)\frac{\beta V^2*}{(V^*+X^*)^2}\\ < /p > < p > \chi'_{2} = (\mu+pZ)\chi_2-\chi_3k(1-T_2^*)-\chi_5cZ^*\\ \chi'_{3} = (1-u_1)(\chi_1-\chi_2)\frac{\beta X^2*}{(V^*+X^*)^2}+\chi_3(\mu+qW^*)-\chi_4gW^*\\ \chi'_{4} = 1+\chi_3 qV^*-\chi_4 (gV^*-h)\\ < /p > < p > \chi'_{5} = 1+pY^*\chi_2-\chi_5 (cY^*-b)\\ T_1^* = min(1, max(0, \frac{1}{B_1}(\chi_2-\chi_1)\frac{\beta X^*V^*}{V^*+X^*})))\\ < /p > < p > T_2^* = min(1, max(0, \frac{1}{B_2}(kY^* \chi_3)))\\ < /p > < p > \chi_j(t_e) = 0, j = 1, ....., 5 \end{cases}$

The outcomes of both the controls that include vaccination & antiviral drug therapy are also represented graphically. Optimal control for linear model can be calculated in the same manner and is given by $T_1^* = min(1, max(0, \frac{1}{A_1}(\chi_2-\chi_1)\beta X^*V^*))$ and $T_2^* = min(1, max(0, \frac{1}{A_2}(\chi_3 kY^*)))$, where $A_1$ and $A_2$ are the benefits and expenses of the applied treatments.

5.   Spatial model for the models

This section would explicitly focus on the formulation of spatial model and its linear analysis. The model is formulated with Holling type-I functional response between healthy cells and virus cells. For, spatial dynamics of the models, we add diffusion to the system

Linear model

Non-linear model

with the following initial conditions and zero-flux boundary conditions:

where $r$ represents the position in $(x, y)$, $(x, y) \in \Omega = [0, L] \times [0, L]$ and $\frac{\partial}{\partial \eta}$ is the normal derivative along unit outward normal vector to $\partial \Omega$. $d_{1}$, $d_{2}$, $d_{3}$, $d_{4}$, $d_{5}$ and $d_{1}'$, $d_{2}'$, $d_{3}'$, $d_{4}'$, $d_{5}'$ are the non-negative diffusion coefficient for linear and non-linear models corresponding to $X(t)$, $Y(t)$, $V(t)$, $W(t)$, and $Z(t)$ respectively. Here $\nabla^{2}$ is the Laplacian operator in two-dimensional space. The description of the rest parameters used in model (2.1) are already described in Table 1.

5.1.   Linear stability analysis and Turing instability

In this subsection, we have studied the local stability of system (5.1). For getting the equilibrium points and their stability, substituting R.H.S of Eq (2.1) is equal to zero, we get

Solving Eq (5.4), gives five set of equilibrium points i.e., $(X_{1}, 0, 0, 0, 0)$, $(X_{2}, Y_{2}, V_{2}, 0, 0)$, $(X_{3}, Y_{3}, V_{3}, 0, Z_{3})$, $(X_{4}, Y_{4}, V_{4}, W_{4}, 0)$, and $(X^{*}, Y^{*}, V^{*}, W^{*}, Z^{*})$. Since we are only interested in coexistence equilibrium point where $X(t)$, $Y(t)$, $V(t)$, $W(t)$, and $Z(t)$ are presents. Hence we have used $E^{*} = (X^{*}, Y^{*}, V^{*}, W^{*}, Z^{*})$ in our analysis.

The variational matrix of the system (2.1) about $E^{*}$ is rendered as:

where

Then, the characteristic equation of (5.5):

where

For the stability of system (2.1), the eigenvalues of Eq (5.6) must have negative real part. As the local stability at $E^{*}$ are very complex to show analytically. So we check the stability by numerical calculation in Numerical section.

Now, we study the effect of diffusion for the spatially homogeneous equilibrium point $E^{*}$ for system (5.1).

Linearizing the system (5.1) about $E^{*}$, gives

where

and $V_{E^{*}}$ is the variational matrix defined in Eq (5.5).

Next, we assume that (5.7) has the solution of the form

where $\sigma$ and $k$ are the growth rate of perturbation and the wave number along $r$ in time $t$, and $c.c.$ represents the complex conjugate. Substituting Eq (5.1) into Eq (5.7), and obtaining characteristic equation:

Solving (5.10), yields the characteristics equation

where

Since, the system (2.1) is stable, therefore the system (5.1) is also stable. For Turing patterns system (5.1) breaks its stability and becomes unstable for some values of diffusion coefficients. Now, we discuss the Turing instability condition of the diffusive model (5.1). Mathematically, Turing instability for (5.1) occurs if at least one root of Eq (5.11) has a positive real part for some $k\ge0$. Therefore, diffusion-driven instability occurs when $\Gamma_{5}(k^{2})\le0$ holds for at least one $k$. Hence, the condition for diffusion instability is given by

where $\Psi_{1} > 0$ and $\Psi_{6} > 0$.

The minimum value of $\Gamma_{6}(k^{2})$ occurs at $k = k_{T}$, where $k_{T}$ is the critical wave number.

Solving the above Eq (5.13), gives the positive value of $k_{T}$. From all the above analytical study, we conclude that for obtaining the Turing patterns $\Gamma_{5}(k^{2})\le0$ for some positive $k$.

After the spatial study of linear model (Urban Area), the same analysis can be carried out for non-linear model (Rural Area) and the results have been discussed in numerical section.

6.   Numerical section

In this section, we would compare the results of both the models numerically and visualize the interpretations. The hypothetical values for parameters used are given in Table 3. Since the density of population in urban area is high as people are now migrating in these areas. Hence, infection rate $\beta$ is assumed more in urban areas in comparison to rural areas assuming same viral load in both the areas. Further, since rural area people have good immune response against the disease due to healthy diet, less population and pollution. Therefore, decay rate of CTl cells ($b$) is assumed to be lower in rural areas in comparison to urban areas.

Comparative study in Urban and Rural area

From Figure 5 and Table 4, we have observed that the greater part of population is getting infected, viral load is increasing and healthy cells are less in the people of urban areas in long run on account of diversely living population and more physical contacts while in rural area, the number of infective are much lesser and healthy cells are increasing. This may be due to less or scattered population in rural areas. Figures 1 and 2 also support our theory. Also Figure 5 and Table 4 depict that the antibodies due to natural body mechanism is higher in urban areas as due to infection the natural body mechanism get activated and CTL cells start protecting the body against infected cells and in case of rural areas since the infection spread is less, the natural body mechanism is not that much activated ^[43]. Figure 6 shows that because of vaccine and anti viral drug administration, the viral load in urban areas decreases resulting in better condition and hence reducing the infective and gradually increasing healthy cells. The same is happening in rural area also but with a very less rate. Also, we have observed that with control the CTL cells and antibody response in rural areas have improved but it is still less than that of urban areas. CTL cells & antibodies play a key role in diminishing viral load and infective and increasing healthy cells. This may be due to no or poor access of medical facilities in rural areas and social stigma is also a reason for lesser impact of control variables ^[44]. From Figures 7 and 8, it is clear that the $\beta$ is directly proportional to $R_0$ as when the value of $\beta$ increases $R_0$ approaches to brightest colour Yellow. Also it is clearly visible that reproduction number in urban area is greater as compare to rural area.

6.1.   Spatial dynamics of linear model

In this section, we perform the numerical simulations to validate our findings for spatial models for the same set of parameters as mentioned in Table 3. For the fixed values of parameters, we plotted a disperse curve between Turing function $(\Gamma_{5}(k^2))$ Vs wave number $k$ for different values of $d_{1}$, $d_{2}$ and $d_{3}$ (see Figure 9) for linear model, which clearly shows that after spatial perturbation the spatio temporal dynamics remains stable. Hence Turing instability condition does not hold here. The same analysis has been done for non-linear system there is also no emergence of Turing instability. So there is no possibility of occurring the Turing patterns for both linear and non-linear model.

Next we see the spatial distribution of system (5.1) in two dimensional plane on solving systems (5.1) and (5.2) by Finite Difference Method and Euler method approximation with time step ($\Delta t = 0.01$) and space step ($\Delta x = 0.5 = \Delta y$). For fixed value of parameters with fixed $d_{1} = d_{2} = d_{3} = d_{4} = d_{5} = 0.01$, the system (5.1) and system (5.2) gives interesting non-Turing patterns (see Figures 10 and 11) for linear model and Figures 12 and 13 for non-linear model. The chosen initial conditions for spatial patterns are:

For Figures 10 and 13

For Figures 11 and 14

For Figures 12 and 15

For linear system

In Figures 12, the spatial distribution of species for fixed value of diffusion and variation in the value of $h$ (decay rate of B-cells) has been studied. As the value of $h$ decreases $0.1 \rightarrow 0.05 \rightarrow 0.01$, population of healthy cells are increasing (Figure 12: $A1\rightarrow B1\rightarrow C1)$, virus load is decreasing (Figure 12: $A3\rightarrow B3\rightarrow C3)$, antibodies are increasing (Figure 12: $A4\rightarrow B4\rightarrow C4)$ and $CTl$ cells are also increasing (Figure 12: $A5\rightarrow B5\rightarrow C5)$. In the spatial distribution the red colour represents high population density and blue colour represents low population density in circular region.

For non-linear system

In Figure 15, the spatial distribution of species for fixed value of diffusion and variation in the value of $b$ (decay rate of $CTL$ -cells) has been studied. As as the value of $b$ decreases $0.2 \rightarrow 0.11 \rightarrow 0.02$, healthy cells are increasing (Figure 15: $A1\rightarrow B1\rightarrow C1)$ where as viral load is decreasing (Figure 15: $A3\rightarrow B3\rightarrow C3)$ and no change take's place in antibodies and $CTL$ cells.

7.   Conclusions and future work

In this paper, we have done a comparative analysis for model $1$ (urban area) and $2$ (rural area) which involves the dynamical analysis of temporal as well as spatial model. We have observed that the rate of infection ($\beta$) between virus cell and uninfected host, decay rate of $CTL$ cells ($b$) are the most sensitive parameters in case of temporal models as they have shown significant impact on influenza in both the areas. From the spatial analysis we have concluded that in urban area the parameter $h$ plays a vital role to cure the viral infections whereas the parameter $b$ plays an important role in rural area. Hence, it can be concluded that in urban area, since the immunity is less, the decay rate of $B$ cells ($h$) and $CTL$ cells ($b$) are the most sensitive parameters in spatial as well as temporal models. In addition, our simulation results also show that availability of treatments like antiviral drugs as well as vaccines play a very important role in curbing influenza in both the areas. Antiviral drug diffuses infected cells in the body, thereby bringing some relief and making the body immune from the infection for a very short span. However, this may not be suitable for massively spread epidemic or pandemic situation. On the other hand, vaccine develops antibodies helps in breaking chain of infection. If, administrated to a population, a vaccine has a tendency to create group immunity or herd immunity leading to major control of an infection. Although, it continues to be difficult to eradicate influenza entirely from a community due to the virus's rapid evolution and ability to shift from one season to the next. Therefore, along with both control mechanisms, it is also important to keep in consideration the areas which are affected majorly due to the infection so that the availability of the antidots and vaccines should be sufficient enough for the people dwelling there.

We would enhance the model by considering different age groups in specific stratified population.

Conflict of interest

The authors declare no conflict of interest.