SS symmetry Article On a SIR Model in a Patchy Environment Under Constant and Feedback Decentralized Controls with Asymmetric Parameterizations Manuel De la Sen 1, * , Asier Ibeas 2 , Santiago Alonso-Quesada 1 and Raul Nistal 1 1 Institute of Research and Development of Processes IIDP, University of the Basque Country, Campus of Leioa, Barrio Sarriena, 48940 Leioa, Bizkaia, Spain; santiago.alonso@ehu.eus (S.A.-Q.); raul.nistal@gmail.com (R.N.) 2 Department of Telecommunications and Systems Engineering, Universitat Autònoma de Barcelona, UAB, 08193 Barcelona, Spain; Asier.Ibeas@uab.cat * Correspondence: manuel.delasen@ehu.eus  Received: 15 February 2019; Accepted: 14 March 2019; Published: 22 March 2019  Abstract: This paper presents a formal description and analysis of an SIR (involving susceptible- infectious-recovered subpopulations) epidemic model in a patchy environment with vaccination controls being constant and proportional to the susceptible subpopulations. The patchy environment is due to the fact that there is a partial interchange of all the subpopulations considered in the model between the various patches what is modelled through the so-called travel matrices. It is assumed that the vaccination controls are administered at each community health centre of a particular patch while either the total information or a partial information of the total subpopulations, including the interchanging ones, is shared by all the set of health centres of the whole environment under study. In the case that not all the information of the subpopulations distributions at other patches are known by the health centre of each particular patch, the feedback vaccination rule would have a decentralized nature. The paper investigates the existence, allocation (depending on the vaccination control gains) and uniqueness of the disease-free equilibrium point as well as the existence of at least a stable endemic equilibrium point. Such a point coincides with the disease-free equilibrium point if the reproduction number is unity. The stability and instability of the disease-free equilibrium point are ensured under the values of the disease reproduction number guaranteeing, respectively, the un-attainability (the reproduction number being less than unity) and stability (the reproduction number being more than unity) of the endemic equilibrium point. The whole set of the potential endemic equilibrium points is characterized and a particular case is also described related to its uniqueness in the case when the patchy model reduces to a unique patch. Vaccination control laws including feedback are proposed which can take into account shared information between the various patches. It is not assumed that there are in the most general case, symmetry-type constrains on the population fluxes between the various patches or in the associated control gains parameterizations. Keywords: epidemic model; irreducible matrix; Metzler matrix; disease transition and transmission matrices; decentralized control; disease-free and endemic equilibrium points; Moore–Penrose pseudoinverse; next generation matrix; patchy environment; vaccination controls 1. Introduction Usually, populations mutually interact through migrations and immigrations to and from other environments. Therefore, the study of more general epidemic models based on interacting subsystems, patches or frame-worked in patchy environments is of a major interest. See, for instance [1–8], and references therein. Then, the implementation of decentralized treatment or vaccination strategies in Symmetry 2019, 11, 430; doi:10.3390/sym11030430 www.mdpi.com/journal/symmetry Symmetry 2019, 11, 430 2 of 42 health centres [9] is of interest, so as to increase their efficiency, by taking into account not only the fixed population assigned to them but also the available information about the fluctuant population associated with migration and punctual travelling. It can be pointed out that the topic of Decentralized Control is very important in a variety of complex problems where control decisions have to be locally taken for the integrated subsystems due to a lack of full information on the coupling dynamics from and to the remaining coupled subsystems taking part of the whole dynamic systems [10–12], the first one concerning with decentralized control while the two last ones are concerned with positivity. In [13], some useful numerical tools are given concerning the non-singularity of perturbed matrices which are used in this paper. Background literature on dynamic systems, including its role on epidemic modelling, is given in [14–19]. In this context, typical situations which need relevant attention when dealing with epidemic models, thinking of their usefulness in their practical implementation in health centers are: (a) The implementation of mixed constant and feedback controls with eventual alternative controller parameterizations and supervisory switching actions between them according to optimization trade-off criteria on the vaccine costs, or their availability, and the infection evolution through time [15,20]. The supervisory scheme chooses online the best appropriate controller parameterization that minimizes the loss function. These considerations could be also of potential applicability interests in the cases of quarantine evaluation on certain parts of the population [17], or occurring transfers from infectious to susceptible individuals [21]. (b) The need for a development of adequate strategies for online either commissioning data [22], or intervention strategies [23], or even the programming of useful strategies for vaccine procurement in due time towards its application to the population [24]. (c) The design of control strategies to fight against the epidemic spreading on multiplex networks which are subject to nonlinear mutual interaction [25], or in cases when the vaccination [16,26–28] is imperfect so that certain amounts of vaccinated susceptible subpopulation are not, in fact, removed from the susceptible subpopulation and transferred to the recovered one. It can be pointed out that patch models have also been used for description of diseases spreading in the real world. In particular, these kind of models have been used to simulate and predict the spatial spreading of infectious diseases. For instance, it is concluded in [29] that the analysis the disease dynamics by considering the effective distances leads to understand complex contagion mechanisms in multiscale networks. The performed analysis showed that network and flux information are sufficient to predict the dynamics and the arrival times. Finally, it was pointed out that the study could be extended to other contagion phenomena, such as activated bio invasion or the spread of rumors. On the other hand, an operational forecast system was developed and verified in [30] that can successfully predict the spatial transmission of influenza in the United States at the state and county levels. On the other hand, we point out that there are other epidemic problems which involve couplings of dynamics between different compartments and subsystems like, for instance, when there are combined diseases and/or the influence of vectors in their propagation. See, for instance [31]. The designed system included processes of surveillance data from multiple locations, forecast accuracy for onset week, peak week, and peak intensity. This paper is focused on the study of the disease-free and endemic equilibrium points as well as the global stability in a patchy environment with multiple patches when there are travelling populations coming into and leaving the various patches. Vaccination strategies are proposed so that each health centre at a particular patch can have and use some certain crossed shared complete or partial information from the remaining patches. It is not assumed, in the most general case, that there are symmetry-type constraints related to the mutual interchanges of populations between pairs of patches or in the control gain parameterizations. The paper is organized as follows. Section 2 describes the proposed SIR epidemic model in a patchy environment of n patches under vaccination control laws which consist of constant and proportional to the susceptible subpopulation actions and which are implemented at each compartment of the patchy structure. Symmetry 2019, 11, 430 3 of 42 The model has travel matrices which take into account the acquisitions and loses of the subpopulations from the other patches due to populations travelling interchanges between each particular patches. The complete model is described in the presence of a feedback vaccination law which contains, in general, constant and feedback linear information on the susceptible subpopulations. It is assumed, in the most general case, that each community health centre can have either a total, a partial, or none information about the susceptible subpopulations of the remaining patches. Such an information can be suitably used, if desired, to generate the whole vaccination control law. Such a law might take into account at each patch not only the subpopulation information of such a concrete patch but, eventually, a total or a partial information of the remaining patches in the whole disposal. These above cases related to the control synthesis rely on the well-known frameworks of centralized control, partially decentralized control, or (fully) decentralized control which are usually invoked in classical Control Theory research [10], especially when the controlled system is complex or distributed in patches which can be physically distributed [10,18,19]. Section 2 also studies the non-negativity of the solutions with initial conditions in the first orthant of the state space and the allocation and uniqueness of the disease-free equilibrium point. Section 3 characterizes the basic reproduction number of the disease by defining the next generation matrix and using its spectral radius as well as the local and global stability and instability properties of the disease-free equilibrium point according to the value of the disease reproduction number compared to unity. The disease-free equilibrium point is calculated as being explicitly dependent on the disease parameters in the model and the control gains. Special particular results are focused on in the cases when some of the relevant travel matrices are irreducible. The endemic equilibrium points are also studied. It is proved that there is at least one endemic equilibrium which is positive and stable (then attainable, that is, allocated within the first orthant of the state space) if the reproduction number equals or exceeds unity. Such an equilibrium point is confluent with the disease- free one if the reproduction number is unity. It is seen, in particular, that if the infectious travel matrix is irreducible, then either all the infectious subpopulation are zero or none of them is zero. This is a very relevant result since with such a kind of conditions, it can be argued that the infectious subpopulations are non-zero at any patches for any endemic equilibrium point. Parallel results are observed in cases when the susceptible travel matrix is irreducible. The characterization of the whole set of endemic equilibrium points is described via the Moore–Penrose pseudoinverse matrix tools [32] by defining a linear algebraic system which contains a partial information of the potential existing set of endemic equilibrium points by neglecting the influence of the quadratic terms associated with the coefficient transmission rates. A complementary nonlinear equation system which is informative about the quadratic terms taking account from the contacts susceptible-infectious in all the patches is then coupled to the above linear system as an extra constraint. If such an algebraic system is compatible indeterminate then there are infinitely many endemic equilibrium solutions including the attainable and un-attainable ones. Section 4 is devoted to the study of the proposed vaccination controls and their implementation in a fully or partly decentralized control context. In particular, the proportional vaccination to the susceptible subpopulation at each patch can be applied only on the susceptible of that patch by taking into account the susceptible subpopulations of those of the other patches which supply it with such an information. The main objective is to distribute the whole set of available vaccines among all the community health centres by sharing such an information. Another potential strategy can be the implementation of vaccination control strategies at each particular health centre of a concrete patch not only on its assigned recorded susceptible but on the travelling susceptible subpopulations coming into it from other patches. Simulated Examples are given and discussed in Section 5. Finally, conclusions end the paper. The proofs of some of the involved results of Section 3 are given in the Appendices A and B. Notation n = {1 , 2 , ... , n}, ei is the i-th unity Euclidean canonical vector of Rn and In is the n-th identity matrix. Symmetry 2019, 11, 430 4 of 42 R+ = R0+ ∪ {0}; R0+ {z ∈ R:z ≥ 0} are the sets of positive and non-negative real numbers, respectively. Z+ = Z0+ ∪ {0}; Z0+ {z ∈ R:z ≥ 0} are the sets of positive and non-negative integer numbers, respectively. A ∈ Rn×n is a Metzler matrix, denoted by A ∈ MEn×n , if all its off-diagonal entries are non-negative.  A  0 (in words, A is non-negative) means that the real matrix A = aij has non-negative entries; A  0 (in words, A is positive) means that aij ≥ 0; ∀i, j ∈ n and there is some (i, j) ∈ n × n such that aij > 0; and A  0 (in words, A is strictly positive) means that all the entries of the real matrix or real vector A are positive. Similar notations are kept for vectors being non-negative (all the components are non-negative), positive (if non-negative with at least one positive component), and strictly positive (all the components are positive). A B, respectively A  B, respectively, A  B means that A − B0, respectively A − B  0, respectively, A − B  0. On the other hand, A ≺ 0 is identical to − A  0, and A ≺ B to B  A. Similar considerations stand “mutatis–mutandis” for the various notations with the symbols “≺”, “≺≺”. ei is the i-th canonical Euclidean vector of the real space Rr whose i-th canonical is unity where the dimension r depends on context. The superscripts T and † stand for transpose and Moore–Penrose pseudoinverses, respectively. If A is a square real non-singular matrix then the transpose of the inverse, identical to inverse of the transpose is denoted by A−T . The symbols ∨ and ∧ stand for logic disjunction and conjunction, respectively. Aij . If A = ( A1 , A2 , .... , An ) T is a real vector, then   If A = Aij is a real matrix | A | = | A| = (| A1 | , | A2 | , .... , | An |)T . If A is a square matrix then ρ( A) is its spectral radius, k Ak2 is the `2 (or spectral) norm and λmax ( A), and respectively, λmin ( A) is its maximum, and respectively, minimum eigenvalue provided that it is real. k Ak1 and k Ak∞ denote, respectively, the `1 and `∞ norms. The time argument in the time-varying variables of differential equations is suppressed for the sake of simplicity when no confusion is expected. We point out that patches could also be referred to as “nodes” (villages, suburbs, towns or regions, each one with a health centre) while “compartment” is each individual subpopulation of susceptible infectious or recovered at each node and “subsystem” is each SIR epidemic mathematical model located at each node in the sense that its describes the self-dynamics at any patch of the whole model including the effects of couplings to other compartments or subsystems. Thus, in our model, the whole system has n subsystems, each one located at one of the n patches, and each subsystem has three compartments, one for each subpopulation. 2. SIR Epidemic Model in a Patchy Environment Under Constant and Proportional Vaccination Controls Consider the following epidemic model in a patchy environment with constant and proportional to the susceptible vaccination controls, which are assumed being monitored in a patchy environment as well: . Si (t) = Λi − β i Si (t) Ii (t) − diS Si (t) + ∑nj(6=i)=1 aij S j (t) − a ji Si (t) − Vi  . (t) I i (t) = β i Si (t) Ii (t) − diI + γi Ii (t) + ∑nj(6=i)=1 bij Ij (t) − b ji Ii (t)   (1) . R n Ri (t) = γi Ii (t) − di Ri (t) + ∑ j(6=i)=1 cij R j (t) − c ji Ri (t) + Vi (t),  ∀i ∈ n, subject to initial conditions Si0 = Si (0) ≥ 0, Ii0 = Ii (0) ≥ 0 and Ri0 = Ri (0) ≥ 0. In the above model, Si (t), Ii (t) and Ri (t) are the susceptible, infectious and recovered (or immune) subpopulations in the i-th patch for i ∈ n, respectively, while β i and γi are, respectively, the disease transmission coefficient rate between susceptible and infectious individuals and the recovery rate of the infectious Symmetry 2019, 11, 430 5 of 42 in the i-th patch. The parameter Λi is the influx of population into the i-th patch. It can be mentioned that in the real word, the influx may also include infectious and immunized subpopulations. However, the influx to infectious and immunized subpopulations is smaller in general than the one to the susceptible subpopulation. In this way, the model only considers the influx affecting the susceptible. The parameters diS , diI and diR are death rates of the susceptible, infectious and recovered, respectively, in the i-th patch. All the parameters of the epidemic model (1) are assumed non-negative and, furthermore, Λi , β i , diS , diI and diR are assumed to be positive for any i ∈ n. The travel matrices    A = aij 0, B = bij 0 and C = cij 0 are not necessarily symmetric and this fact does not affect to the problem formulation. Note that the immigration and outmigration amounts are proportional to the subpopulation values at the various patches. However, the stationary populations never reach zero values at any patch if the respective influx term is nonzero. The description of (1) can be made through the susceptible, infectious and recovered vectors S(t) = (S1 (t) , S2 (t) , ... , Sn (t)) T , I (t) = ( I1 (t) , I2 (t) , ... , In (t)) T and R(t) = ( R1 (t) , R2 (t) , ... , Rn (t)) T , respectively. The vaccination controls are assumed to be monitored via linear feedback information from the susceptible and have the form: n Vi (t) = Vi0 + ∑ Kij Sj (t), i = 1, 2, · · · , n(n ≥ 2) (2) j =1 for given prefixed control gains Kij . The replacement of (2) into (1) yields: . Si (t) = − β i Si (t) Ii (t) − diS Si (t) + ∑nj(6=i)=1aij − Kij S j (t) − a ji Si (t) + Λi − Vi0 − Kii Si (t)   . I i (t) = β i Si (t) Ii (t) − diI + γi Ii (t) + ∑nj(6=i)=1 bij Ij (t) − b ji Ii (t)   (3) . Ri (t) = γi Ii (t) − diR Ri (t) + ∑nj(6=i)=1 cij R j (t) − c ji Ri (t) + Kij S j (t) + Vi0 + Kii Si (t),  ∀i ∈ n. In the sequel, and for the sake of simplicity, the dependence of the variables from time is deleted in the notation when no confusion is expected. The first part of the subsequent result relies on the existence, uniqueness and attainability (or reachability), in the sense that it has no negative component, of the disease-free equilibrium point. The second part of such a result establishes that, for identically zero infection levels through time, the disease-free equilibrium point is globally exponentially stable. The proof is based on the fact that the opposed matrix to an M-matrix is a Metzler matrix and a Metzler matrix is a stability matrix if and only if it is non-singular and its minus inverse is positive: Theorem 1. Define two real vectors P and Λ and a real square matrix D as follows: " # h iT h iT DSS DSR P = ST , RT ∈ R2n ; Λ = ΛST , Λ TR ∈ R2n ; D = ∈ R2n×2n (4) DRS DRR where: S = [S1 , S2 , ... , Sn ] T ; R = [ R1 , R2 , ... , Rn ] T (5) ΛS = Λ − Λ R ; Λ = [Λ1 , Λ2 , · · · , Λn ] T ; Λ R = V0 = [V10 , V20 , · · · , Vn0 ] T (6) d1 + ∑nj=2 a j1 + K11  S K12 − a12 ··· K1n − a1n  S n   K21 − a21 d2 + ∑ j(6=2)=1 a j2 + K22 ··· K2n − a2n   DSS = .. .. .. ..  (7) .    . . .  Kn1 − an1 Kn2 − an2 · · · dSn + ∑nj=−11 a jn + Knn Symmetry 2019, 11, 430 6 of 42 d1R + ∑nj=2 c j1 −c12 ··· −c1n     −c21 d2R + ∑nj(6=2)=1 c j2 ··· −c2n   DRR =  .. .. .. ..  .    . . .  −cn1 −cn2 · · · dnR + ∑nj=−11 c jn   (8) −K11 −K12 · · · −K1n  −K21 −K22 · · · −K2n   ; DSR = 0 ∈ Rn×n  DRS = −K =  .. .. .. .. . . . .     −Kn1 −Kn2 · · · −Knn and assume that the control gains are fixed as follows:  Vi0 ∈ [0 , Λi ] ; Kij ∈ 0 , aij ; ∀i, j(6= i ) ∈ n   (9) Kii > − diS + ∑nj(6=i)=1 aij ; Kii ≥ −∑nj(6=i)=1 Kij ; ∀i, j(6= i ) ∈ n such that Vi0 = Λi for some i ∈ n. Then, the following properties hold: (i) The disease-free equilibrium point of Equation (1), under the vaccination control Equation (2) exists, it is unique and attainable, and given by  T T T T   xd∗ f = xd∗1f , xd∗2f , .... , xd∗nf T ; xd∗if = Sid ∗ ∗ f , 0 , Rid f ; ∀i ∈ n (10) ∗ = e T S∗ , R∗ = e T R∗ ; ∀i ∈ n, where: with Sid f i df id f i df  T Sd∗ f = S1d ∗ ∗ ∗ f , S2d f , · · · , Snd f −1 = DSS ΛS (11)  T     R∗d f = R1d ∗ ∗ ∗ f , R2d f , · · · , Rnd f −1 = DRR Λ R + | DRS | Sd∗ f = DRR −1 −1 | DRS | DSS ΛS + V0 leading to a disease-free equilibrium total population vector:  T   Nd∗f = N1d ∗ f , N ∗ 2d f , · · · , N ∗ nd f = S ∗ df + R ∗ df = In + D −1 RR | D RS | −1 DSS −1 ΛS + DRR V0 (12) and, in the particular case that di = diS = diR ; ∀i ∈ n to the following disease-free equilibrium total population amount: n n Λ NTd∗ f = ∑ id f ∑ dii N ∗ = (13) i =1 i =1 This limit total population is also reached under any existing endemic equilibrium points. Furthermore, the total population N (t) is bounded for any finite initial conditions and all t ≥ 0. (ii) The solution trajectory of the linearized system around the disease-free equilibrium point of the model Equation (3) within the zero-infective (I ≡ 0 ∈ Rn )2n-dimensional subspace of R3n is non-negative for any non-negative initial conditions Si (0), Ri (0); ∀i ∈ n and it is also globally exponentially stable irrespective of the vaccination controls. Proof. Note that the epidemic model (1) is subject to the parametrical constraints that Λi , β i , diS , diI and diR are positive for any i ∈ n, and A = aij 0, B = bij 0 and C = cij 0 under the vaccination    controls (2) subject to (9). Therefore each two terms aii Si and each two terms cii Ri , with opposed signs, become cancelled, respectively, in the first and third equation of Equations (3) for all i ∈ n. Then, Symmetry 2019, 11, 430 7 of 42 one can fix aii = cii = 0 for i ∈ n in Equations (7) and (8) with no loss in generality by keeping the summations from one to n. The disease-free equilibrium point satisfies the constraints: n − diS Si + ∑ aij − Kij S j − a ji S j + Λi − Vi0 = 0   j =1 n − diR Ri + ∑ cij R j − c ji Ri + Kij S j + Vi0 = 0;  j =1 ∀i ∈ n, by fixing aii = cii = 0 for i ∈ n. Note that DRR has non-positive off-diagonal entries with the sum of all the entries per column being positive. Thus, it is a non-singular M-matrix −1 with DRR  0. Also, DSS is has non-positive off-diagonal entries with the sum of all the entries −1 per column being positive from Equation (9). Thus, it is a non-singular M-matrix with DSS  0 [1]. Furthermore, − DRS = | DRS |  0. Therefore, the disease-free equilibrium point is unique and defined by Equations (10) and (11) subject to Equations (4)–(9). The total disease-free equilibrium population Equation (12) follows directly from Equation (11) and the disease-free total population vector is Nd∗f = Sd∗ f + R∗d f . It is attainable in the sense that it has no negative components and it is also nonzero, since DSS and DRR are non-singular from Equation (11), subject to Equations (4)–(9). Equation (13) follows since the total population satisfies the constraint: . n n   n h  i NT = ∑ Ni = ∑ Si + I i + Ri = ∑ Λi − diS Si + diI Ii + diR Ri i =1 i =1 i =1 and, for the disease-free equilibrium point with di = diS = diR ; ∀i ∈ n, . ∗ h i ∗ N Td f = ∑in=1 Λi − di Nid n n f + ∑hi =1 ∑ j=1     aij S j − a ji Si + bij Ij − b ji Ii + cij R j − c ji Ri i = ∑in=1 Λi − di Nid ∗ n f + 0 = ∑ i =1 (0 ) = 0 so that Nd∗f = ∑in=1 Nid∗ = n Λi . It follows that N ( t ) is bounded for any finite initial conditions f ∑ i =1 d i ∗ for all t ≥ 0 and N (t) → NTd f as t → ∞ . Property (i) has been proved. To prove Property (ii), first note that the Jacobian matrix of the linearized system (1), subject to Equation (2), or equivalently Equation (3), about xd∗ f within the manifold I ≡ 0 is Jd∗ f = − D. Since the conditions Equations (9) hold then D is an M-matrix with D −1  0. Thus, Jd∗ f ∈ MEn×n so that the linearized solution trajectory is non-negative for any given set of non-negative initial conditions since a time-invariant linear system has a non-negative solution trajectory irrespective of any given non-negative initial conditions if and only if its matrix of dynamics is a Metzler matrix [11,12]. Furthermore, the Jacobian matrix is invertible −1 satisfying − Jd∗ f = D −1  0. Since a Metzler matrix is a stability matrix if and only if it is non-singular and its minus inverse is positive, one concludes that the linearized system around the disease-free equilibrium point is globally exponentially stable since it is time-invariant so that the asymptotic stability is also exponential.  If, for generality purposes and coherency with the generality of the model, it is supposed in Theorem 1 (i), Equation (13), that, in general, di 6= diR , with di = diR = diS + dei ; ∀i ∈ n in the sense that if the parameters differ from each other, then the mortality of the recovered who already suffered the disease is slightly higher than that of the susceptible since they suffered from the illness. Thus, one gets: n n Λi − dei R∗ n id f Λ NTd f = ∑ Nid f = ∑ ∗ ∗ (<) ≈ ∑ i i =1 i =1 d i d i =1 i Symmetry 2019, 11, 430 8 of 42 Remark 1. Note from Equation (1) and Equation (2) that if Ii (0) = 0; for some i ∈ n then Ii (t) = 0; ∀i ∈ n, t ≥ 0. Under these conditions Theorem 1 (ii) applies. Remark 2. Note from Equations (2), (3), (4) and (9) that, although Kij ≥ 0; ∀i ∈ n in the vaccination law, it is not requested for any particular gain Kii to be positive. The subsequent result relies on some disease-free equilibrium point results based on the positivity and irreducibility of some relevant travel matrices and constraints on the vaccination control describing population fluxes between patches of the model. Theorem 2. The following properties hold:  (i) Assume that B = bij is irreducible. Then, Ii (t) = 0; ∀t ∈ [t1 , t2 ] for some i ∈ n implies that Ij (t) = 0; ∀t ∈ [t1 , t2 ], ∀ j ∈ n irrespectively of the vaccination control law. (ii) Assume that Vi0 = Λi ; ∀i ∈ n and assume also that A − K = aij − Kij is irreducible with A  K.   Then, S j (t) = 0; ∀t ∈ [t1 , t2 ], ∀ j ∈ n if Si (t) = 0; ∀t ∈ [t1 , t2 ] for some i ∈ n. If B = bij  and C = cij are irreducible, K = 0 and Vi0 = 0; ∀i ∈ n then R j (t) = 0; ∀t ∈ [t1 , t2 ], ∀ j ∈ n if Ri (t) = Ii (t) = 0; ∀t ∈ [t1 , t2 ] for some i ∈ n.   (iii) Assume that the conditions of Property (ii) hold and that, furthermore, Kij ∈ 0 , aij ; ∀i, j(6= i ) ∈ n,   Kii > − diS + ∑nj(6=i)=1 aij and Kii ≥ −∑nj(6=i)=1 Kij ; ∀i ∈ n. Then, Nd∗f = R∗d f and NTd ∗ f = ∗ , that is the total population is recovered at the disease-free equilibrium point. ∑in=1 Rid f . Proof. Assume that Ii (t) = 0; ∀t ∈ [t1 , t2 ], then I i (t) = 0; ∀t ∈ (t1 , t2 ) for some i ∈ n, ∀t ∈ (t1 , t2 ) and assume also that there are j(6= i ) ∈ n and t ∈ [t1 , t2 ] such that Ij (t) 6= 0. One concludes . from the second equation of (3), if Ii (t) = 0 for t ∈ [t1 , t2 ], so that  I i (t) = 0 for t ∈ (t1 , t2 ), that ∑ j(6=i)=1 bij Ij (t) = ∑ j=1 bij Ij (t) = BI (t) = 0; ∀t ∈ [t1 , t2 ]. Then, ∑nj=−01 B j I (t) = 0; ∀t ∈ [t1 , t2 ]. n n   But B is irreducible if and only if ∑nj=−01 B j  0, since B  0, and then ∑nj=−01 B j I (t)  0 for any t ∈ [t1 , t2 ] if there is at least one Ij (t) 6= 0 for some j(6= i ) ∈ n and some t ∈ [t1 , t2 ], a contradiction . to ∑nj(6=i)=1 bij Ij (t) = 0; ∀t ∈ [t1 , t2 ]. Then, Ij (t) = 0; ∀t ∈ [t1 , t2 ] so that I j (t) = 0; ∀t ∈ (t1 , t2 ), ∀ j ∈ n. Property (i) has been proved. On the other hand, one concludes from the first equation of (3) if . Si (t) = 0 for t ∈ [t1 , t2 ], so that Si (t) = 0 for t ∈ (t1 , t2 ), that n n ∑ ∑   aij − Kij S j = aij − Kij S j = 0; ∀t ∈ [t1 , t2 ] j(6=i )=1 j =1 . provided that Vi0 = Λi ; ∀i ∈ n and, if Ri (t) = 0 and Ii (t) = 0 for t ∈ [t1 , t2 ], so that Ri (t) = 0 for t ∈ (t1 , t2 ), one concludes that, if in addition Vi0 = 0; ∀i ∈ n then ∑nj(6=i)=1 cij R j − c ji Ri = 0.  The proof of Property (ii) is completed under similar reasoning as that used in the proof of Property (i). Finally, Property (iii) follows directly from Property (ii) and Theorem 1 (i) via Equation (9).  It has to be pointed out that a particular version of Theorem 2 (i) for the case of absence of vaccination controls has been proved in another way in [1]. In the total absence of vaccination parameterized by the vector Ω = 0, the vectors and matrices of Equations (4)–(8) are subject to the following replacements Λ R → 0 , DSS → DSS0 , DRS → 0 ; and DRR and DSR = 0 are kept identical with: Symmetry 2019, 11, 430 9 of 42 d1S + ∑nj=2 a j1 − a12 ··· − a1n     − a21 d2S + ∑nj(6=2)=1 a j2 ··· − a2n   DSS0 = DSS ]Ω=0 = .. .. .. ..  .    . . .  − an1 − an2 · · · dSn + ∑nj=−11 a jn 3. Basic Reproduction Number: Attainability of the Endemic Equilibrium versus Instability of the Disease-Free One Define the following matrices:   ∗ ∗ ∗ F = Diag β 1 S1d f , β 2 S 2d f , · · · , β n S nd f (14) d1I + γ1 + ∑nj=2 b j1 −b12 ··· −b1n     −b21 d2I + γ2 + ∑nj(6=2)=1 b j2 · · · −b2n   U= .. .. .. ..  (15) .    . . .  −bn1 −bn2 · · · dnI + γn + ∑nj=−11 b jn The basic reproduction number is R0 = ρ FU −1 , where (−U ) is the transition matrix, F is the  transmission matrix and FU −1 is the next generation matrix. The following positivity and stability result, proven in Appendix A, holds: Theorem 3. The following properties hold: (i) (−U ) ∈ MEn×n is stability matrix. (ii) If β i = 0; ∀i ∈ n then the disease-free equilibrium point is globally exponentially stable and any solution trajectory is non-negative for all time for any given non-negative initial conditions. (iii) If R0 < 1 then the disease-free equilibrium point xd∗ f is locally asymptotically stable and, if R0 > 1, such an equilibrium point is unstable. (iv) The reproduction number satisfies the subsequent upper-bounding constraint:   −1 −1 R0 ≤ R01 = β max ( β ir )k DSS (Λ − V0 )k2 ρ1/2 U T U (16) 1≤ i ≤ n where β ir = β i /β; ∀i ∈ n, are relative transmission coefficient rates. Assume, in addition, that kU0d k2 < 1/kUd−1 k2 and k DSS0d k2 < 1/k DSSd −1 k2 , where Ud and Uod are the diagonal and off-diagonal parts part of U = Ud + U0d and DSSd and DSSod are the diagonal and off-diagonal parts of DSS . Then, −1 k DSSd k2 kUd−1 k2 R0 ≤ R02 = β max ( β ir ) −1 kΛ − V0 k2 (17) 1≤ i ≤ n 1 − k DSSd k2 k DSS0d k2 1 − kUd−1 k2 kU0d k2 with β ≥ 0 being a prefixed reference value of the coefficient transmission rate. (v) R02 is minimized for any given model parameterization and any given constant vaccination vector V0 if the vaccination control gains for the susceptible are chosen as Kij = aij ; ∀i, j ∈ n\{1}. Such a reproduction number upper-bound is zeroed if each whole influx of population in all patches are vaccinated by constant controls. Remark 3. Note that β can be, in practice, one of the coefficient rates (for instance, its maximum or minimum value). Note that the choice β = 0 is feasible if and only if β i = 0; ∀i ∈ n. The non-negativity of the linearized solution proved in Theorem 1 (ii) also applies to the whole non-linear system under weak conditions as follows. Symmetry 2019, 11, 430 10 of 42 Theorem 4. Assume that the vaccination control constrains Equations (9) hold and that AK. Then, the following properties hold: (i) Any solution trajectory of the whole non-linear system Equation (1) is non-negative and bounded for all time for any given finite non-negative initial conditions (ii) Assume, furthermore, that R0 ≥ 1. Then, there exists at least one endemic equilibrium point. If, in addition, B  0 then any endemic equilibrium point has a positive infective population at any patch. If A − K  0 is irreducible then any endemic equilibrium point has a positive susceptible population at any patch even under a maximum constant vaccination Vi0 = Λi ; ∀i ∈ n. (iii) There is no attainable endemic equilibrium point if R0 < 1 while, if R0 ≤ 1, then the unique disease-free equilibrium point is globally asymptotically stable. If R0 = 1 then such a disease-free equilibrium point coincides with one of the existing attainable endemic equilibrium points. Proof. From Theorem 1 (i), the total population N (t) is bounded for all time. By inspecting Equation (1), one concludes that if any susceptible, infectious or recovered subpopulation at any patch and time instant is zero then its time-derivative cannot be negative since AK, B0 and C 0 and Equation (9) hold. Therefore, min (Si (t) , Ii (t) , Ri (t)) ≥ 0 ⇒ min (Si (t) , Ii (t) , Ri (t)) ≥ 0 ; ∀t ≥ 0. i ∈n i ∈n If, furthermore, max (Si (0) , Ii (0) , Ri (0)) < +∞ then, sup max (Si (t) , Ii (t) , Ri (t)) < +∞ since i ∈n t∈R0+ i ∈n N (t) < +∞; ∀t ≥ 0. Property (i) has been proved. Property (ii) is proved by contradiction for the case R0 > 1. Assume that R0 > 1 and since no endemic equilibrium point exists. Thus, the disease-free equilibrium point is unstable, any state solution trajectory has bounded non-negative components for any time and any finite non-negative initial conditions, and no endemic equilibrium point exists. Thus, it follows from Poincaré’s index that a stable bounded limit cycle should surround the disease-free equilibrium point which is the unique (unstable) equilibrium point which has a unity Poincaré’s index. But this feature contradicts that the state solution trajectory is non-negative for all time and any non-negative initial conditions so that no stable limit cycle can surround the unstable disease-free equilibrium point. Therefore, at least one endemic equilibrium point muss exist if R0 > 1. The first part of Property (ii) has been proved. Now, if, in addition, B is irreducible then any zero infectious subpopulation at any patch implies that the infectious total population is zero from Theorem 2 (i). By its equivalent contra-positive implication logic proposition, since the endemic equilibrium point has a nonzero total infectious population, any endemic equilibrium infectious subpopulation is nonzero at any patch. Thus, the infectious subpopulation is nonzero at any patch at the endemic equilibrium points. It follows in the same way that, if ( A − K )  0 is irreducible, then the endemic susceptible subpopulation has to be nonzero at any patch. Property (ii) has been proved for R0 > 1. Now, assume that R0 = 1. In this case, the disease-free equilibrium point is critically stable so that it has at least either one centre (i.e., a critical point with two imaginary complex eigenvalues in one of the two-dimensional partial Jacobian matrices) or one spurious patch (i.e., a critical point with one zero eigenvalue and the other one real positive in one of the two-dimensional partial Jacobian matrices) in at least a two-dimensional hyperplane of the phase space. This situation is also incompatible with the non-negativity of the solution trajectory so that the conclusion on the existence of an endemic equilibrium point is similar to the former part of the proof of this property. Proposition (ii) has been proved. To prove Property (iii), assume that there is an attainable (i.e., with no negative component) endemic equilibrium point if R0 < 1 and note, from Equations (1), (14) and (15), that ∗ Iiend + eiT ( F − U )−1 (− β 1 S1end ∗ ∗ I1end ∗ , − β 2 S2end ∗ I2end ∗ , . . . , − β n Snend ∗ Inend ) T ; ∀i ∈ n (18) Symmetry 2019, 11, 430 11 of 42 where ( F − U )−1 exists and −( F − U )−1  0 since ( F − U ) ∈ MEn×n is a stability matrix since (−U ) ∈ MEn×n is a stability matrix, so U −1  0, and R0 = ρ FU −1 < 1. Thus, ( F − U )−1 has at least  one positive entry per column and one positive entry per row. Then, the above equation holds for ∗ min β i > 0 with Iiend > 0; ∀i ∈ n if and only if S∗jend < 0 for at least a j ∈ n. Thus, there is no attainable i ∈n endemic equilibrium point if R0 < 1 and min β i > 0. Since an endemic equilibrium point exists for i ∈n R0 = 1 from Property (ii), the fact that Equation (18) also holds for R0 = 1, as a result, and the fact that the subsequent constraint stands for the disease-free equilibrium point if R0 < 1:   ∗ Iidf = eiT (−U )−1 − β 1 S1df ∗ ∗ ∗ ∗ I1df , − β 2 S2df ∗ I2df , . . . , − β n Sndf ∗ Indf = 0; ∀i ∈ n (19) it follows from continuity arguments of the equilibrium points with respect to R0 that one of the endemic equilibrium points necessarily coincide with the disease-free one for R0 = 1. Now since: (a) the disease-free equilibrium point is unique and the unique attainable equilibrium point for R0 < 1 (Theorem 1 (i)); and (b) such a point is furthermore locally asymptotically stable, since its linearized version around it is asymptotically stable (Theorem 3 (iii), one concludes that the disease-free equilibrium point is globally asymptotically stable if R0 ≤ 1. Property (iii) has been proved.  Remark 4. Theorem 4 (ii) establishes that, if the disease-free equilibrium point is unstable or critically stable, then an endemic equilibrium point has to exist. With some extra irreducibility-type conditions on the B-travel matrix and on the ( A − K )-travel matrix, it is proved that the infectious and susceptible endemic equilibrium amounts are nonzero at any patch. It can be argued that the matrix of proportional vaccination gains K can modify the irreducibility or reducibility properties of the travel matrix A related to the respective properties of ( A − K ). This fact can imply that, if in the absence of proportional vaccination to the susceptible subpopulation, the endemic equilibrium point has nonzero susceptible (respectively, zero amounts of susceptible at least at one patch) subpopulations at any patch, then, under some kind of proportional vaccination law even for a constant vaccination constraint Vi0 = Λi ; ∀i ∈ n, the endemic susceptible could be zeroed at least at one patch but not in all ! patches. To visualize the above argument, note that the matrix constraint i ∑in=−01 Ai  −∑in=−01 ∑ij=1 Ai− j (−K ) j guarantees that ( A − K ) is irreducible since j n −1 i n −1 i ! n −1 n −1 i ! i i ∑ ( A − K) = ∑∑ j A i− j (−K ) =j ∑ i A + ∑∑ j Ai− j (−K ) j  0. i =0 i =0 j =0 i =0 i =0 j =1 The characterization of the whole set of endemic equilibrium points is addressed in the following result, which is proved in Appendix B, by using algebraic tools: Theorem 5. Assume that R0 ≥ 1 and define the following matrices:     a11 − a12 · · · − a1n b11 −b12 · · · −b1n  − a21 a22 · · · − a2n   −b21 b22 · · · −b2n  AS =  ; A I =  (20)     ··· ··· ··· ··· ··· ··· ··· ···      − an1 − an2 · · · ann −bn1 −bn2 · · · bnn   c11 −c12 · · · −c1n  −c21 c22 · · · −c2n  A RI = Diag [−γ1 , −γ2 , . . . , −γn ]; A R =  (21)   ··· ··· ··· ···    −cn1 −cn2 · · · cnn Symmetry 2019, 11, 430 12 of 42 n   ∗ ∗ ∑  ∗ ∗ Λi − Vi0 = aii Siend + bii Iiend − aij − Kij Sjend + bij Ijend (22) j(6=i )=1 n ∗ Vi0 = cii Riend ∗ + γi Iiend − ∑ ∗ cij Rjend ∗ + Kij Sjend (23) j(6=i )=1 where n n n aii = diS + Kii + ∑ a ji , bii = diI + γi + ∑ b ji , cii = diR + ∑ c ji ; ∀i ∈ n (24) j(6=i )=1 j(6=i )=1 j(6=i )=1 Then, the following properties hold: (i) The following rank condition holds: rank(b, A) = rank A (25) where the limit total population is N ∗ irrespective of the equilibrium point as time tends to infinity, and N∗     Λ1 − V10      ..  1 1 ··· 1 1   .   (2n+1)×3n A= AS A I 0 ∈R ; b= Λn − Vn0  ∈ R2n+1 (26)       0 A RI A R  V10  ..      .  Vn0 The whole set of endemic equilibrium solutions, including both the attainable and unattainable ones, is given by x(y) = A† b + (I3n − A† A)y (27) subject to the n constraints: h  i diI + γi + ∑nj(6=i)=1 b ji enT+i − ∑nj(6=i)=1 bij enT+ j A† b + I3n − A† A y   βi = T ; ∀i ∈ n (28) eiT [ A† b + ( I3n − A† A)y][ A† b + ( I3n − A† A)y] en+i ∗ (y), S∗ (y), . . . S∗ (y), I ∗ (y), I ∗ (y), . . . I ∗ (y), R∗ (y), R∗ (y), . . . R∗ (y) T  with x(y) = S1end 2end nend 1end 2end nend 1end 2end nend and ei is the Euclidean canonical vector whose its ith component is unity; ∀i ∈ 3n,  −1 T  −1 T An = D T DD T C C C ∈ R3n×(2n+1) is the Moore–Penrose pseudoinverse of A, provided that A of rank p ≤ 2n + 1 is factorized as A = CD with existing matrices C ∈ R2(n+1)× p and D ∈ R p×(2n+1) both or rank p, and y ∈ R3n is arbitrary except that it is subject to fulfill Equation (28) for the given coefficient transmission rates β i for i ∈ n, where A† T = A T [32]. The set of attainable endemic equilibrium points is given by Equation (27) subject to the constraints Equation (28) for any y ∈ Y with Y = {z ∈ R3n :x(z)(∈ R3n ) 0}.  (ii) If B = bij is irreducible then the set of attainable endemic equilibrium points is given by (27), subject to the constraints (28), for any y ∈ Ya with Ya = {z ∈ R3n :(x(z)(∈ R3n ) 0) ∧ (xi (z) > 0; i = n + 1, n + 2, ..., 2n)} ⊂ Y (iii) Vi0 = Λi ; ∀i ∈ n and both B = bij and A − K = aij − Kij , with A  K, are irreducible then the set   of attainable endemic equilibrium points is given by Equation (27), subject to the constraints Equation (28), for any y ∈ Yb with Yb = {z ∈ R3n :(x(z)(∈ R3n ) 0) ∧ ([xi (z) > 0; i = 1,2,..., n] ∨ [xi (z) = 0; i = 1,2,..., n]) ∧ (xi (z) > 0; i = n+1,n+2,..., 2n)} ⊂ Ya Symmetry 2019, 11, 430 13 of 42 (iv) If K = 0, Vi0 = Λi = 0; ∀i ∈ n and B = bij , A − K = aij − Kij , with A  K, and C = cij are    irreducible, then the set of attainable endemic equilibrium points is given by Equation (27) subject to the constraints Equation (28) for any y ∈ Yc ⊂ Yb with Yc = z ∈ R3n : (`1 ∧ `2 ∧ `3 ∧ `4 ) holds    `1 := x (z) ∈ R3n  0 `2 := ( xi (z) > 0; i = 1, 2, ..., n) ∨ ( xi (z) = 0; i = 1, 2, ..., n) `3 := ( xi (z) > 0; i = n + 1, n + 2, ..., 2n + 1) `4 := ( xi (z) > 0; i = 2n + 1, 2n + 2, ..., 3n) ∨ ( xi (z) = 0; i = 2n + 1, 2n + 2, ..., 3n) The conditions for the uniqueness of the existing attainable endemic equilibrium point for R0 ≥ 1 are given in the following result which is a direct conclusion of Theorem 5: Corollary 1. Assume that R0 ≥ 1. Then, the attainable equilibrium point is unique if and only there is a y ∈ R3n such that (1) y + A† (b − Ay)  0, E y + A† (b − Ay)  0(respectively,  0 if B is irreducible), where E = [0n×n 0n×n ] ∈ Rn×3n ,  (2) In×n (3) The n constraints (28) hold. One such a vector y ∈ R3n always exists. The following counterpart result to Theorem 5 and Corollary 1 holds for the case when there is only one patch in the epidemic model so that the transportation matrices are zero. The result, proved in Appendix B, gives a nice physical interpretation of the basic reproduction number and its relation to the stability properties and to the attainability of the endemic equilibrium point. Theorem 6. Assume that there is only one patch (i.e., n = 1) and that Λ > V with V being a constant vaccination effort. Then, there is a unique stable attainable endemic equilibrium point if the coefficient dS (d I +γ) S∗ transmission rate fulfills β ≥ β c = Λ−V , equivalently, if the reproduction number, R0 = S∗d f ≥ 1, end Λ −V where Sd∗ f = dS is the susceptible subpopulation at the disease- free equilibrium point, the immune one at the disease-free equilibrium being R∗d f = V dR . Such an endemic equilibrium point is: ∗ = d I +γ ∗ = β(Λ−V )−dS (d I +γ) β(d I +γ)V +γ[ β(Λ−V )−dS (d I +γ) ] Send β ; Iend ;R∗end = (29) β( d I +γ ) βd R (d I +γ) And the following properties hold: (i) If Λ = V then there is a unique disease-free equilibrium point Sd∗ f = Id∗ f = 0 ; R∗d f = dVR while the endemic one does not exist. (ii) If R0 = 1 then the disease- free and the endemic equilibrium points coincide. (iii) If R0 < 1 then the disease- free equilibrium point is globally asymptotically stable and the endemic one is not attainable. (iv) If V (t) = V0 + KS(t) then Sd∗ f = Λ −V0 dS +K , R∗d f = KΛ +dS V0 , and d R (dS +K ) ! Λ dS − d R dS V0 + KΛ Λ K dS − d R dS − d R V0    Nd∗f = S+ = S 1+ R S + R S d dS d R (dS + K ) d d (d + K ) d (d + K ) Λ In the absence of vaccination, Nd∗f = Sd∗ f = dS and R∗d f = 0. Symmetry 2019, 11, 430 14 of 42 The following result, which is proved in Appendix C, relies on the feature that the reproduction number can be reduced by the vaccination controls. This feature implies that the global asymptotic stability towards the disease-free equilibrium point can be guaranteed under smaller values of the coefficient transmission rates via an appropriate monitoring of such controls. Although the proposed model has an identical transmission matrix U for the vaccination-free and vaccinated models, it is assumed for analysis generality purposes that that associated to the vaccination case Uc can be distinct to that associated to the vaccination-free one Uun . This is the case, for instance, if an additional treatment control is injected on the infectious subpopulation. See, for instance [14,15].   Theorem 7. Define Uc = Uun + U e and Fc = Fun + F,e where Fe and −U e are the disturbed transmission and transition matrix of the controlled epidemic model under a vaccination control law with respect to those of −1 and  the uncontrolled (i.e., for the case when the vaccination control is null) one. Define R0un = ρ Fun Uun R0c = ρ Fc Uc−1 as the respective reproduction numbers in the vaccination-free and under vaccination. Assume  that the following constraints hold: ×3n (1) (−Uun ) ∈ M3n E is a stability matrix, (2) Fun  0, (3) kU e k < 1/2 kUun−1 k , 2 2     −1  −1 (4) − 1 − 1 − 1 − Fun ≺ F ≺ Fun Uun U I3n + Uun U I3n − U I3n + Uun U e e e e e − Uun1 Uun . ×3n Then, Uc ∈ M3n E is a stability matrix and the following properties hold: (i) R0c ≤ R0un . (ii) If, the conditions Equations (1)–(3) hold, Fe = − Fe ≺ 0 and the constraint equation (4) is replaced with following constraints:     −1  −1 (40 ) − Fun Uun −1 U −1 U e I3n + Uun e I3n − U e I3n + Uun−1 U e −1 Uun Uun ≺ Fe ≺ Fun . Then R0c ≤ R0un . In addition, R0c < R0un if either Fun Uun −1 or F −1 is irreducible. This property e Uun result still holds if one but not both) of the two “≺”-symbols of the above equation is replaced with “≺”. Remark 5. Note that the applicability of Theorem 7 (ii) is very feasible in practice according to the following considerations. Assume that the pairs ( Fun , Uun ) and ( Fc , Uc ) are the pairs defining the vaccination-free and vaccination cases linear dynamics around the disease-free equilibrium point which depends on the control gains such that U = Uc = Uun from (14) and (15) for the model dealt with.(Note that Theorem 7 has been worked for the more general case when Uc 6= Uun ). Now, Fe = − Fe ≺ 0 if Fc ≺ Fun , that is, if S∗ ( Fc ) < S∗ ( Fun ). This df df is directly achievable by using appropriate control gains (see Theorem 1). In the simplest case of just one patch in the model (i.e., n = 1), note that this is achievable by choosing max (V0 , K ) > 0 from Theorem 6 (iv). The choices of the values of the control gains V0 and K monitor the susceptible amounts Sd∗ f ( Fc ) at the disease- free equilibrium. Now, assume that R0un = 1. This value of the reproduction number corresponds to a certain critical disease transmission rate β cun for given remaining modeling parameters in the vaccination-free case. This fact leads to the coincidence of the disease-free equilibrium point with the attainable endemic one and the critical stability of the disease-free equilibrium point. However, under Theorem 7, and since Fe < 0, the vaccination control leads to the asymptotic stability of the modified disease-free equilibrium point and the un-attainability of the endemic one since R0c < R0un = 1. Therefore, a properly designed vaccination law increases the range of the stability boundary of the disease-free equilibrium point to reach a larger critical disease transmission rate compared to the vaccination-free case. 4. Use of Available Patch-Crossed Information in Decentralized Vaccination Control Designs The following situations can occur related to the vaccination controls monitoring actions: Symmetry 2019, 11, 430 15 of 42 (a) Centralized Vaccination Control (CVC). Each subsystem has the information available about the susceptible numbers of all the compartments and uses it for feedback vaccination control. (b) Decentralized Vaccination Control (DVC) if Kij = 0; ∀i, j(6= i ) ∈ n and Kii 6= 0; ∀i ∈ n. Each subsystem uses only self-information for control but there is no use of the susceptible number of other compartments. (c) Partially Decentralized Vaccination Control (PDVC) if Kii 6= 0; ∀i ∈ n, Kij 6= 0; ∀(i, j) ∈ n p × nq and Kij = 0; ∀(i, j) ∈ n × n\n p × nq , where n p and nq are nonempty proper subsets of n. (d) nw -Weak Decentralized Vaccination Control (nw -WDVC) if Kij = 0; ∀i, j(6= i ) ∈ n, Kii 6= 0; ∀i ∈ n p and Kii = 0; ∀i ∈ n\ n p . That is, at least one compartment of susceptible does not uses susceptible self-information for feedback in the vaccination control law which has a decentralized structure. (e) nw -Weak Partially Decentralized Vaccination Control (nw -WPDVC) if in the definition of nw -WDVC, Kij 6= 0 for some i, j(6= i ) ∈ n. Note that the various concepts of “centralized control” versus “decentralized control” refer to the complete or partial shared information between dynamic subsystems and, in particular, subsystems of the patchy model or just the use of own self- information for control rather than to the physical disposal (generic one or local for each subsystem) of the controller. This is a widely admitted principle in decentralized control of dynamic systems. See, for instance [10]. Two vaccination strategies are now discussed if the vaccination controls are assumed to be monitored via linear feedback information from the susceptible by using available information at each patch from some other patches: Strategy 1. Only the susceptible subpopulation of each patch, even if travelling population from other patches exists, is a candidate to be vaccinated while some total or partial information from the corresponding subpopulations in other patches is known and monitored for the susceptible vaccination through the crossed control gains associated with the control law (2). Such an information is used to restrict the influence of the immigration from the remaining patches into the own susceptible subpopulation of a patch in accordance with Equation (3). The control law Equation (2) is assumed to be subject to the following constraints: n n 0 ≤ Kii ≤ Mi + ∑ K ji − ∑ Kij , 0 ≤ K ji ≤ a ji , Vi0 ≤ Mi0 < Λi ; ∀i, j(6= i ) ∈ n (30) j(6=i )=1 j(6=i )=1 where Mi > 0 and Mi0 > 0are upper-bounding constant taking into account the vaccines availability at the i-th patch for i ∈ n. The first constraint of Equation (30) reflects that a fraction of the travelling susceptible populations coming from the remaining patches is vaccinated while the leaving one to other patches is not vaccinated. The second constraint takes into account that DSS in Equation (7) is an M-matrix so that its inverse exists and is positive, so that the disease-free equilibrium point is a non-negative vector of the state space and locally asymptotically stable since (− DSS ) ∈ MEn×n . Strategy 2. Only the susceptible subpopulation proper of each patch is a candidate for vaccination but there is some partial or total information from the susceptible subpopulations from other patches. The available information on the coming in and leaving travelling susceptible subpopulations from the various patches is used to control the distribution of the vaccines to be administrated between the various patches. Such an information is used to restrict the number of administered vaccines at each patch. In this case, the vaccination control law Equation (2) is modified as follows: Vi (t) = Vi0 + Ki (t)Si (t); ∀i ∈ n (31) and the vaccination control proportional gains are given by: n Ki (t) = Ki (S(t)) = Kii + ∑ Kij0 (t); ∀i ∈ n (32) j(6=i )=1 Symmetry 2019, 11, 430 16 of 42 where: Kij S j (t)    Si ( t ) if Si ( t ) > ε i Kij0 (t) Kij0  = Si ( t ) , S j ( t ) = ; ∀i, j ∈ n (33)  0 if Si ( t ) ≤ ε i  for given prefixed control gains Kij and design constants ε i ∈ R0+ ; ∀i, j ∈ n. It turns out from Equations (31)–(33) that coupled information between distinct patch pairs can be available or not in the vaccination controls. As a result, the vaccination control (31)–(33) becomes: ( Vi0 + ∑nj=1 Kij S j (t) if Si (t) > ε i Vi (t) = ; ∀i ∈ n (34) Vi0 + Kii Si (t) if Si (t) ≤ ε i The constraints Equation (30) become modified as follows for each i ∈ n allowing some negative crossed control gains: Kij ≤ 0; 0 ≤ K ji ≤ a ji ; ; ∀ j(6= i ) ∈ n (35)  S j (t)  n  S j (t)  Kii n ≥ sup max ; Kii ≤ Mi + ∑ K ji inf min (36) ∑ j(6=i)=1 ij K t∈R0+ 1≤ j(6=i )≤n S i ( t ) j(6=i )=1 t ∈ R 0+ 1≤ j (6 =i )≤ n Si ( t ) Note that Equations (35) and (36) may be jointly expressed as follows:   n  S j (t)  n  S j (t)   ∑ Kii  sup max Si ( t ) ≤ Kii ≤ Mi + ∑ K ji inf min t∈R0+ 1≤ j(6=i )≤n Si ( t ) (37) j(6=i )=1 t∈R0+ 1≤ j(6=i )≤n j(6=i )=1 provided that the following necessary condition holds:   n M  ∑ Kij  ≤  S (t)  i j  S (t)  j (38) j(6=i )=1 sup max S (t) − inf min Si ( t ) t∈R0+ 1≤ j(6=i )≤n i t∈R0+ 1≤ j(6=i )≤n Note the following facts: (1) If Si (t) = ε i and Si (t+ ) > ε i fore some i ∈ n then Vi (t) switches from a constant term to a combined constant plus a linear feedback term except if the control gains Kij = 0; ∀ j ∈ n and such a i ∈ n. In this case, the closed-loop linearized dynamic systems around any potential equilibrium points, which are defined by their corresponding Jacobian matrices at such points after absorbing the linear feedback from the susceptible subpopulations, are not time-invariant through time. (2) If either in f Si (t) > ε i or sup Si (t) ≤ ε i ; ∀ i ∈ n, then the vaccination control law does not t ∈ R0 + t ∈ R0 + switch from a combined constant plus a linear feedback term to a constant term or vice-versa at any patch and at any time instant. (3) Concerning the Centralized/Decentralized control frameworks, note that a CVC strategy is implementable if the available information allows the use of gains Kij 6= 0; ∀i, j(6= i ) ∈ n since all the susceptible subpopulation and its distribution between the various patches is known at each patch. A PDVC, or a DVC strategy is adopted when some or, respectively, all the gains Kij are zeroed; ∀i, j(6= i ) ∈ n because the global information on susceptible is not known, or not used, at each patch. The (nw -WDVC) and (nw -WPDVC) vaccination strategies are implemented if some of the self-proportional gains are not used at some patches (i.e., there is no vaccination action at some health centre on its own susceptible subpopulation) or, if, in addition some of the crossed susceptible information between the various patches is not available or simply not used. It can be convenient to adopt vaccination strategies which allow to guarantee a worst-case minimization, in some sense, of the disease-free equilibrium subpopulations in order to achieve Symmetry 2019, 11, 430 17 of 42 a corresponding maximization of the recovered subpopulation when the infection is removed. This idea is addressed in the sequel. Note that   n k DSS k1 = max Kii + dSI + ∑  2aij − Kij  (39) 1 ≤i ≤ n j(6=i )=1   n k DSS k∞ = max Kii + dSI + ∑  aij + a ji − K ji  (40) 1 ≤i ≤ n j(6=i )=1 Then, one has from (11) via Equations (6) and (7) and using the constraints (30) for Strategy 1, by taking into account the bounded relations between the matrix and vector spectral (`2 ) and `1 and `∞ norms, that the following lower-bounds stand for the disease-free equilibrium susceptible vector: max (Λi −Vi0 ) kΛS k∞ kSd∗ f k = max Sid ∗ ≥ f k DSS k∞ = 1 ≤i ≤ n h i ∞ 1≤ i ≤ n max Kii +dSI +∑nj(6=i)=1 aij + a ji −K ji( ) 1 ≤i ≤ n (41) max (Λi − Mi0 ) 1 ≤i ≤ n ≥ Mi −∑nj(6=i)=1 Kij +dSI +∑nj(6=i)=1 ( aij + a ji ) k Λ S k1 ∑in=1 (Λi −Vi0 ) kSd∗ f k = ∑in=1 Sid ∗ ≥ f k DSS k1 = h i 1 max Kii +dSI +∑nj(6=i)=1 (2aij −Kij ) 1 ≤i ≤ n ∑in=1 (Λi −V ∑in=1 (Λi −Vi0 )  (42) i0 ) ≥ h i ≥ h i max Mi +∑nj(6=i)=1 K ji +dSI +2 ∑nj(6=i)=1 ( aij −Kij ) max Mi +∑ j(6=i)=1 K ji +dSI +2 ∑nj(6=i)=1 aij n 1 ≤i ≤ n 1 ≤i ≤ n q 2 k Λ S k2 k Λ S k2 ∑in=1 (Λi −Vi0 ) q   kSd∗ f k = ∑in=1 Sid ∗2 ≥ k DSS k2 = ≥ n max 1 k DSS k1 , 1 k DSS k ∞ (43) 2 f λ1/2 T max ( DSS DSS ) Remark 6. In view of Equations (41)–(43), one concludes that available lower-bounds susceptible subpopulations at the disease-free equilibrium points can be reduced in a suboptimal worst-case design which keeps the maximum available vaccines and jointly minimizes the `1 , `∞ and `2 norms by choosing: Vi0 = Mi0 ; Kij = 0; K ji = a ji ; ∀i, j ∈ n n n n Kii = Mi + ∑ K ji − ∑ Kij = Mi + ∑ a ji ; ∀i ∈ n j(6=i )=1 j(6=i )=1 j(6=i )=1 In the case that some outsider travelers from other patches to a certain patch i ∈ n have to be vaccinated for needs of global fulfillment of objectives, one can use normalizing factors `ij ∈ [0 , 1] so that Kij = `ij aij replaces the standard strategy Kij = 0; ∀ j ∈ n. In the case that some travelers from a certain patch i ∈ n to other patches should be vaccinated, one can use normalizing factors ` ji ∈ [0 , 1] so that K ji = ` ji a ji replaces the standard strategy K ji = a ji ; ∀ j ∈ n. Note from (31) to (34) that, in the case of Strategy 2, the vaccination control parameterization is time-varying (see, for instance [20]), since there can exist switches if the susceptible subpopulation at any patch is close to zero. The following two technical results are of usefulness for Strategy 2. Lemma 1. Let A ∈ Rn×n be a stability matrix of stability abscissa −ρ a < 0 and let be A e : R0 + → R n × n a piecewise continuous uniformly bounded matrix function. Then, the matrix function B : R0+ → Rn×n , being e(t); ∀t ∈ R0+ is stable if (ρ a /Ka ) t > t A R B(t) = A + A e(τ ) dτ; ∀t ∈ R+ , guaranteed if sup k A e(t)k < 0 t ∈ R0 + ρa Ka , for some norm-dependent real constant Ka ≥ 1. Symmetry 2019, 11, 430 18 of 42 . Proof. Consider the n-th differential system z(t) = B(t)z(t); z(0) = z0 with kz0 k < ∞. It turns out that there exists Ka ≥ 1 such that  Z t  k z(t)k ≤ Ka e−ρa t k z0 k + eρa τ k A e(τ )k kz(τ )k dτ ; ∀t ∈ R0+ (44) 0 Rt so that kz(t)k ≤ Ka kz0 k e− 0 (ρa −Ka k A(τ )k)dτ which follows from (44), the constraint (ρ a /Ka ) t > e Rt 0 A ( τ ) dτ; ∀ t ∈ R+ and Gronwall’s Lemma [33] so that k z ( t )k ≤ K a k z0 k; ∀ t ∈ R0+ and z ( t ) → 0 e as t → ∞ .  Rt The condition (ρ a /Ka ) t > e(τ ) dτ of Lemma 1 may be weakened to (ρ a /Ka ) (t − t0 ) > A 0 Rt e(τ ) dτ for any t(> t0 ) ∈ R+ and some t0 ∈ R0+ . Lemma 1 yields to the following result: A t0 Theorem 8. Consider (14) and (15) with −U ∈ MEn×n a stability matrix and F ( 0) ∈ Rn×n such that ρ FU −1 < 1 and let Fe : R0+ → Rn×n be uniformly bounded piecewise continuous and asymptotically  convergent to Fee ∈ Rn×n . Then, there exists some norm-dependent real constant Ka ≥ 1 such that ρ ( F −U ) F + Fe − U : R0+ → Rn×n is stable provided that sup k Fe(t)k < Ka . t ∈ R0 + .   If, furthermore, F (t) − F; ∀t ∈ R0+ then the differential system y(t) = F + Fe(t) − U y(t) is positive e in the sense that it has a solution trajectory within the first open orthant of the state space for any initial condition y(0) = y0  0. Proof. Since −U ∈ MEn×n , F ( 0) ∈ Rn×n and ρ FU −1 < 1 then ( F − U ) ∈ MEn×n so that it has a  maximal real eigenvalue which is stable since ( F − U ) is stable since −U is stable and ρ FU −1 < 1.  Thus, the minus stability abscissa of ( F − U ) is also its spectral radius, that is, ρ a ( F − U ) = ρ( F − U ) ρ ( F −U ) and ke( F−U )t k ≤ Ka e−ρt for any t ∈ R and some Ka ≥ 1. If sup k Fe(t)k < Ka for such an existing t ∈ R0 +   norm-dependent real constant Ka , then one has that the time-varying matrix F + Fe(t) − U is stable   from Lemma 1 and it converges asymptotically to the stability matrix F + Fee − U . On the other hand, .   the differential system y(t) = F + Fe(t) − U y(t) has a unique solution for any given y(0) = y0 ∈ Rn given by: Z t   −Ut y(t) = e y0 + e−U (t−τ ) F + Fe(τ ) y(τ )dτ (45) 0 Since −U ∈ MEn×n  thenh e−Ut  0 for any t i ∈ R0+ [12]. Now, note by direct inspection of Equation (45) that [y0 0] ∧ F (t) − F; ∀t ∈ R0+ ⇒ (y(t) 0 ; ∀t ∈ R0+ ) .  e Remark 7. A practical implementation of the vaccination control law Equations (31)–(33) is to choose the design constants ε i for i ∈ n being very close to zero and to make null all the proportional vaccination gains Kij0 (t) at patch i for the crossed susceptible information from other patches j 6= i and any t ≥ ti in the event that S(ti ) < ε i at some time instant ti . In this way, the maximum number of switches is n, the last eventual one occurring in a finite time T f . Then, the stability conditions of Theorem 8 are simplified to simpler conditions for h  ρ ( F −U ) a time-invariant system on T f , +∞ by deleting the conditions sup k Fe(t)k < Ka and Fe(t) → Fee as t ∈ R0 + h  t → ∞ , since Fe(t) = Fee ; ∀t ≥ T f and the finite time interval 0 , T f is irrelevant for stability analysis, and modifying the condition ρ FU −1 < 1 to ρ ( F + Fe )U −1 < 1.   Symmetry 2019, 11, 430 19 of 42 5. Simulation Symmetry 2019, 11,Examples x FOR PEER REVIEW 20 of 43 This section contains some numerical simulation examples related to the results presented in the 5. Simulation Examples previous sections. The examples are concerned with the existence of equilibrium points along with the effectThis section of the contains vaccination somestrategies control numericalproposed simulation examples in Section related 4 on to the results the epidemic presented spreading. in In this the previous case, it will besections. shown howThethe examples are concerned vaccination controllerswith the to are able existence of equilibrium reduce the incidence ofpoints along an infection with the within effect of the vaccination control strategies proposed in Section 4 on the epidemic spreading. a population. In this case, it will be shown how the vaccination controllers are able to reduce the incidence of an Example infection 1. Consider within the SIR patchy system defined by three patches or populations, n = 3, with parameters a population. given by: Example 1. Consider  X  hthe SIR patchy system defined i by three patchesh or populations, n = i3 , with parameters given by:d = di = 1/3 1/3.1 1/3.2 years 1 ,β = [ β i ] = 3.24 3.08 3.16 × 10−2 − [ ] h i Λ = 30d, = [ ] = 1.78 1.82 1.75 d = diX = [1/ 3 1/ 3.1 1/ 3.2]years−1 , β = [β i ] = [3 . 24 3 . 08 3 . 16 ] × 10 − 2 γ γ i in units of week−1 except otherwise indicated. Λ = 30d , The symbol d X stands for any parameter dS , d I , d R . Notice that γ = [γ i ] = [1.78 1.82 1.75] it is very typical that different outbreaks of the same epidemic have different reproduction numbers [34,35] since in units of week −1of the the spreading epidemic, except and therefore otherwise indicated.itsThe severity, symbol d X on depends manyfor stands factors any such as the geographical parameter d S ,d I ,d R . distribution of the individuals, the probability of an infected individual contact a healthy Notice that it is very typical that different outbreaks of the same epidemic have different reproduction one, etc. The initial numbers conditions are the [34,35] since given by: spreading of the epidemic, and therefore its severity, depends on many factors such as the S1 (0) = 25; geographical distribution of the individuals, I1 (0) = 10; the probability 1 (0) = individual of anRinfected 0 contact a healthy one, etc. The initial conditions are given by: 2 S ( 0 ) = 30; I2 ( 0 ) = 10; R 2 ( 0 ) = 0 S3 (0) =S20; 3 (0 ( 0 ) =I25 ;I) (= = 10;R 0 )5; R3((00)) ==0 0 1 1 1 while the travel matrices are given by: S2 ( 0 ) = 30;I 2 ( 0 ) = 10;R2 ( 0 ) = 0   ( 0 ) = 20;I 3( 0 ) = 5;R3( 0 S3 )=0   0 0.2 while the travel matrices 0.3 by: are given 0 0.22 0.4 0 0.17 0.25 A =  0.16 0 0.3 , B =  0.15 0 0.05 , C =  0.3 0 0.12         0 0.2 0.3   0 0.22 0.4   0 0.17 0.25  0.35 0.14 0  0.15 0.15 0  0.3 0.2  0 A =  0.16 0 0.3  , B =  0.15 0 0.05  , C =  0.3 0 0.12   0.35 0.14 0   0.15 0.15 0   0.3 0.2 0  The dynamics of the system without vaccination is depicted in Figures 1–3: The dynamics of the system without vaccination is depicted in Figures 1–3: Figure 1. Evolution of the susceptible within each patch without vaccination. Figure 1. Evolution of the susceptible within each patch without vaccination. Symmetry 2019, 11, x FOR PEER REVIEW 21 of 43 Symmetry 2019,11, Symmetry 2019, 11,430 x FOR PEER REVIEW 20 of 21 of 42 43 Figure 2. Evolution of the infectious within each patch without vaccination. Figure 2. Evolution of the infectious within each patch without vaccination. Figure 2. Evolution of the infectious within each patch without vaccination. Figure 3. Evolution of the immune within each patch without vaccination. Figure 3. Evolution of the immune within each patch without vaccination. From Figures 1–3 it Figure can be observed 3. Evolution that the of the immune withinaboveeachparameters patch without correspond vaccination.to the case when From Figuresnumber the reproduction 1–3 it can be observed is less than unity,thatRthe01> thanunity, 1. In this case, the disease-free equilibrium 0 . In this case, the disease-free equilibrium point point is unstable and an asymptotically stable endemic equilibrium point appears. The following Figures 7–9 is unstable and an asymptotically stable endemic equilibrium point appears. The following Figures 7–9 display display the evolution of the system in this case when no vaccination is applied. the evolution of the system in this case when no vaccination is applied. Figure 7. Evolution of the susceptible in all patches when R0 > 1. Figure 7. Evolution of the susceptible in all patches when R0 >1. Symmetry 2019, 11, x FOR PEER REVIEW 24 of 43 Symmetry 2019,11, Symmetry2019, 11,430 x FOR PEER REVIEW 23 of 24 of 42 43 Figure 8. Evolution of the infectious in all patches when R0 >1. Figure 8. Evolution of the infectious in all patches when R0 > 1. Figure 8. Evolution of the infectious in all patches when R0 >1. Figure 9. Evolution of the immune in all patches when R0 > 1. Figure 9. Evolution of the immune in all patches when R0 >1. It Figure that9.theEvolution of thedo immune in all patches when R0 >1. equilibrium point It can can be be observed observed that the infectious infectious do not not vanish vanish now. The now. Theendemic endemic equilibrium point is is given by (Send1 , Send2 , Send3 ) = (7.61, 9.26, 8.36), ((IIend1 ,, IIend2,,IIend3))=(4 = ( 4.03, 3.19, .03, 3.19, 2.70 ) 2.70 ), and given by (Send1,Send2,Send3)=(7.61,9.26, 8.36) , end1 end2The end3endemic , and ( Rend1It, Rcan end2 be , R observed end3 ) = ( that 16.40, the infectious 18.89, 19.67 ) . Ado not series vanish of now. numerical experiments equilibrium are conducted point now is (Rend1,Rend2,Rend3 )=(16.40,18.89,19.67). A series of numerical experiments are conducted now to analyze the to analyze by given the effect of1,S (Send end2,Send3)=(7.61 parameters and,9.26 , 8.36)conditions initial , ( in Iendthe 1, Iend 2, Iend3 ) = location of(4the .03, 3endemic .19, 2.70) point. , Thus, and effect of parameters and initial conditions in the location of the endemic point. Thus, the initial (Rend1,Rend2,Rend3 )=(16.40,18.89,19.67). A series of numerical experiments are conducted now to analyze the the initial values of the populations are now changed to: values of the populations are now changed to: € effect of parameters and initial conditions in the location of the endemic point. Thus, the initial S1 (0) = 55; S1(0)I=155 (0;)I1=(0)15;=15;R1R(10()0=)2= 2 values of the populations are now changed to: € S2 (0) = 40; I2 (0) = 8; R2 (0) = 1 S3 (0) = 22;SS12((00))=I=355 40;;II2((00))==15 (0)1= 5; ;R 8;R2((00))= R13 (0)=2 1 =2 SS3((00))==22 0) = 40; I32(0)=8;R23(0)=12 ; I ( 5; R ( 0 ) = 2 The The evolution of the system with different initialconditions evolution of the system with different initial is shown in Figures 10–12. S3(0) =22; I3(0)conditions =5; R3(0) =2is shown in Figures 10–12. The evolution of the system with different initial conditions is shown in Figures 10–12. Symmetry 2019, 11, x FOR PEER REVIEW 25 of 43 Symmetry 2019, 11, x FOR PEER REVIEW 25 of 43 Symmetry 2019, 11, Symmetry 2019, 11, 430 x FOR PEER REVIEW 24 of 25 of 42 43 Figure 10. Evolution of the susceptible in all patches when R0 >1 and different initial conditions. Figure 10. Evolution of the susceptible in all patches when R0 >1 and different initial conditions. Figure 10. Evolution of the susceptible in all patches when R0 > 1 and different initial conditions. Figure 10. Evolution of the susceptible in all patches when R0 >1 and different initial conditions. Figure 11. Evolution of the infectious in all patches when R0 >1 and different initial conditions. Figure 11. Evolution of the infectious in all patches when R0 > 1 and different initial conditions. Figure 11. Evolution of the infectious in all patches when R0 >1 and different initial conditions. Figure 11. Evolution of the infectious in all patches when R0 >1 and different initial conditions. Figure 12. Evolution of the immune in all patches when R0 > 1 and different initial conditions. Figure 12. Evolution of the immune in all patches when R0 >1 and different initial conditions. TheFigure endemic equilibrium 12. Evolution point of the is given immune bypatches in all the same when R0 >1indicated values before. and different Thus, initial the location of conditions. the endemic equilibrium point is not altered by a change in the initial values. Afterwards, the value Figure 12. Evolution of β 3 is perturbed of the (while the immune others in allβ 2patches β 1 and remain when R0 >1 andand unchanged) different initial conditions. the location of the endemic equilibrium point for each case is provided in Table 3. Symmetry 2019, 11, 430 25 of 42 Table 3. Location of the endemic equilibrium point for different values of β 3 . β3 (Send1 ,Send2 ,Send3 ) (Iend1 ,Iend2 ,Iend3 ) (Rend1 ,Rend2 ,Rend3 ) 28.44 × 10−2 (7.51, 9.22, 7.57) (4.00, 3.10, 2.96) (16.53, 18.85, 20.40) 37.92 × 10−2 (7.29, 9.13, 5.86) (3.93, 2.91, 3.53) (16.79, 18.74, 21.95) 63.20 × 10−2 (7.00, 9.00, 3.61) (3.84, 2.67, 4.28) (17.15, 18.61, 23.99) 94.80 × 10−2 (6.84, 8.92, 2.44) (3.80, 2.55, 4.67) (17.34, 18.55, 25.06) As it can be deduced from Table 3, the location of the endemic equilibrium point changes according to the change in β 3 . To conclude this example, consider now the values of ( β 1 , β 2 , β 3 ) included in Table 4 and the corresponding endemic points. Table 4. Location of the endemic equilibrium point for β = 29.92 × 10−2 . ( β 1 ,β 2 ,β 3 ) (Send1 ,Send2 ,Send3 ) (Iend1 ,Iend2 ,Iend3 ) (Rend1 ,Rend2 ,Rend3 ) ( β, β, β) (7.56, 8.83, 8.16) (4.00, 3.27, 2.74) (16.46, 19.19, 19.91) (10β, β, β) (0.83, 8.37, 7.31) (6.09, 2.98, 2.20) (20.95, 23.52, 20.87) ( β, 10β, β) (7.01, 0.92, 7.65) (3.61, 5.23, 2.51) (17.11, 24.86, 21.24) ( β, β, 10β) (6.61, 8.42, 0.90) (3.73, 2.48, 5.16) (17.62, 18.76, 26.50) It can be observed in Table 4 how the location of the endemic point changes as the value 10β moves from one position to another one within the vector [ β 1 , β 2 , β 3 ]. Overall, it is concluded that the endemic point does not change with variations of initial conditions, but it generally does with parameter changes. Example 3. Finally, consider the Hong Kong influenza epidemic in New York City in 1968–1969. This influenza outbreak is modeled by an SIR epidemic model with the following parameters [36]: β = 3.24 × 10−7 , γ = 1.78 in units of week−1 . The patchy environment is inspired on this real case and it is composed of three cities (or patches), n = 3, with spreading parameters similar to the above ones and given by: h i h i Λ = [Λi ] = 5 4.5 5.5 × 103 , β = [ β i ] = 3.24 3.18 3.08 × 10−7   h i h i d X = diX = 1/70 1/71 1/72 years−1 , γ = [γi ] = 1.78 1.82 1.75 in units of week−1 except otherwise indicated and the symbol d X stands for dS , d I , d R . The initial conditions for the populations are given by the 1970 New York City census as: S1 (0) = 7, 960, 000; I1 (0) = 15, 000; R1 (0) = 0 while the initial conditions for the remaining patches are given, similarly, by: S2 (0) = 8, 600, 000; I2 (0) = 20, 000; R2 (0) = 0 S3 (0) = 7, 200, 000; I3 (0) = 19, 000; R3 (0) = 0 The travel matrices are defined by:       0 1.2 0.3 0 1.12 0.4 0 0.78 0.56 A = 10−2 ×  1.1 0 1 , B = 10−2 ×  1.22 0 0.85 , C = 10−2 ×  1 0 0.95        1.2 1.4 0 1 1.14 0 1.2 0.94 0 Symmetry 2019, 11, x FOR PEER REVIEW 27 of 43 Symmetry 2019, 11, x FOR PEER REVIEW 27 of 43 0 1.2 0.3  0 1.12 0.4  0 0.78 0.56  0 1.2 0.3  0 1.12 0.4   0 0.78 0.56 Symmetry 2019, 11, 430 A = 10−−22 × 1.1 0 1 , B = 10−−22 × 1.22 0 0.85, C = 10−−22 ×  1 0 0.95 26 of 42 A = 10 × 1.1 0 1 , B = 10 × 1.22 0 0.85, C = 10 ×  1 0 0.95 1.2 1.4 0   1 1.14 0  1.2 0.94 0  1.2 1.4 0   1 1.14 0  1.2 0.94 0     TheThe aimaim of this example of this is toisshow example the effect to show of the the effect ofvaccination strategies the vaccination introduced strategies in Section introduced 4. in Section The aim of this example is to show the effect of the vaccination strategies introduced in Section The4.evolution of the system without vaccination is displayed in Figures 13–15. The evolution of the system without vaccination is displayed in Figures 13–15. 4. The evolution of the system without vaccination is displayed in Figures 13–15. Figure 13. Evolution 13. Evolution of the of the susceptible subpopulation within eacheach patch. Figure Figure 13. Evolution of susceptible subpopulation the susceptible within subpopulation within patch. each patch. Figure 14. Evolution 14. Evolution of the of the infectious subpopulation within eacheach patch. Figure Figure 14. Evolution of infectious subpopulation the infectious within subpopulation within patch. each patch. Symmetry 2019, 11, 430 27 of 42 Symmetry 2019, 11, x FOR PEER REVIEW 28 of 43 Figure 15. Evolution of the immune subpopulation within each patch. Figure 15. Evolution of the immune subpopulation within each patch. As it can be observed in Figure 14, the influenza outbreak reaches a peak during the spreading of As it can be observed in Figure 14, the influenza outbreak reaches a peak during the spreading of the infection. In order to reduce the severity of the outbreak, the two vaccination strategies proposed the infection. In order to reduce the severity of the outbreak, the two vaccination strategies proposed in in Section 4 are now applied and compared. To this end, consider the control matrices given by: Section 4 are now applied and compared. To this end, consider the control matrices given by: h i K = A + Diag 10−2 , 0.6 ×  10 −2 , 0.9 × 10 −2 ; Mi = 5 × 105 ; M0 = 0.9Λ; V0 = M0 K = A + Diag 10−2 , 0.6×10−2 ,0.9×10−2  ; Mi = 5×105; M0 = 0.9Λ;V0 = M0   It can be readily seen that the above selection satisfies the constraints imposed by (30). Moreover, the thresholds It can to bebereadily used inseen Strategy that 2the are above given by selection × 106 ; the ε 1 = 4.3satisfies 5.1 × 106 ; imposed ε 2 =constraints 6. ε 3 = 4.7 ×by10(30). The Figures 16–21 display the evolution of various infectious subpopulations6in agreement 6 with the 6 Moreover, the thresholds to be used in Strategy 2 are given by ε1 = 4.3×10 ; ε2 = 5.1×10 ; ε3 = 4.7×10 . implemented vaccination controls. The Figure 16, Figure 18, and Figure 20 show the evolution of The Figures 16–21 display the evolution of various infectious subpopulations in agreement with the the infectious subpopulation at each patch without vaccination and when both vaccination strategies implemented vaccination controls. The Figures 16, 18, and 20 show the evolution of the infectious introduced in Section 4 are employed. Furthermore, the Figure 16, Figure 18, and Figure 20 show the subpopulation at each patch without vaccination and when both vaccination strategies introduced in vaccination commands generated by both strategies at each patch. It can be seen that the solution Section 4 are employed. Furthermore, the Figures 16, 18, and 20 show the vaccination commands trajectory of the infectious is non-negative and globally bounded as it is proved in Theorem 4. From generated by both strategies at each patch. It can be seen that the solution trajectory of the infectious Figure 16, Figure 18, and Figure 20 it can also be concluded that the application of a judicious is non-negative and globally bounded as it is proved in Theorem 4. From Figures 16, 18, and 20 it can vaccination campaign significantly reduces the peak caused by the outbreak. In addition, Figure 17, also be concluded that the application of a judicious vaccination campaign significantly reduces the Figure 19, and Figure 21 show that Strategies 1 and 2 generate very similar infectious subpopulation peak caused by the outbreak. In addition, Figures 17, 19, and 21 show that Strategies 1 and 2 profiles, where the plots for both cases are almost superimposed. However, the vaccination law profile generate very similar infectious subpopulation profiles, where the plots for both cases are almost through time is different for Strategies 1 and 2, fact that can be observed in Figure 17, Figure 19, and superimposed. However, the vaccination law profile through time is different for Strategies 1 and 2, Figure 21. During the first weeks, both control laws are the same but when the susceptible reach the fact that can be observed in Figures 17, 19, and 21. During the first weeks, both control laws are the corresponding prescribed threshold, the susceptible feedback term of Strategy 2’s vaccination law is same but when the susceptible reach the corresponding prescribed threshold, the susceptible switched off and only a constant vaccination is applied. The shutting down of the feedback term causes feedback term of Strategy 2’s vaccination law is switched off and only a constant vaccination is a noticeable decrease of the control command while the evolution of the infectious subpopulations is applied. The shutting down of the feedback term causes a noticeable decrease of the control similar. Consequently, the vaccination Strategy 2 is able to reduce the outbreak peak, saving vaccination command while the evolution of the infectious subpopulations is similar. Consequently, the effort. Notice that, in this experiment, each patch disposes of full information of the remaining ones vaccination Strategy 2 is able to reduce the outbreak peak, saving vaccination effort. Notice that, in since the values of the susceptible subpopulation at the others patches are used to calculate the amount this experiment, each patch disposes of full information of the remaining ones since the values of the of vaccination according to Equations (31)–(33). susceptible subpopulation at the others patches are used to calculate the amount of vaccination according to Equations (31)–(33). Symmetry Symmetry2019, 11,11, 2019, 430x FOR PEER REVIEW 28 29of of42 43 43 Symmetry 2019, 11, x FOR PEER REVIEW 29 of Figure Figure 16. Evolution 16. Evolution of of of the infectious thethe infectious subpopulation subpopulation within patch 1within under patch 1vaccination different under different strategies. Figure 16. Evolution infectious subpopulation within patch 1 under different vaccination strategies. vaccination strategies. 1400 1400 Vaccination strategy Vaccination 1 1 strategy Vaccination strategy Vaccination 2 2 strategy 1200 1200 1000 1000 Vaccination in Patch 1 Vaccination in Patch 1 800800 600600 400400 200200 0 0 0 0 2 2 4 4 6 6 8 8 10 10 12 12 14 14 Time (weeks) Time (weeks) Figure 17. Vaccination law in patch 1 for Strategies 1 and 2. Figure 17.17. Figure Vaccination law Vaccination in in law patch 1 for patch Strategies 1 for 1 and Strategies 2. 2. 1 and Symmetry Symmetry2019, 2019,11, 11,430 x FOR PEER REVIEW 29 30 of of42 43 43 Symmetry 2019, 11, x FOR PEER REVIEW 30 of Figure Figure 18.18. 18. Figure Evolution Evolution of of Evolution thetheofinfectious the infectious infectious subpopulation subpopulation within subpopulation patch within 2within patchunder patch different 2 under 2 vaccination under different vaccination different strategies. strategies. vaccination strategies. 1400 1400 Vaccination strategy Vaccination 1 1 strategy Vaccination strategy Vaccination 2 2 strategy 1200 1200 1000 1000 Vaccination in Patch 2 Vaccination in Patch 2 800800 600600 400400 200200 0 0 0 0 2 2 4 4 6 6 8 8 10 10 12 12 14 14 Time (weeks) Time (weeks) Figure Figure19. Vaccination 19.19. Vaccination law lawininpatch 2 2for patch forStrategies Strategies1 1and and 2.2. Figure Vaccination law in patch 2 for Strategies 1 and 2. Symmetry 2019, 11, 430 30 of 42 Symmetry Symmetry 2019, 2019, 11, 11, x FOR x FOR PEER PEER REVIEW REVIEW 31 31 of 43 of 43 Figure 20. Evolution of the infectious subpopulation within patch 3 under different Figure Figure 20.20. Evolution of the infectious subpopulation within patch 3 under different vaccination strategies. Evolution vaccination strategies.of the infectious subpopulation within patch 3 under different vaccination strategies. 2000 2000 Vaccination Vaccination strategy strategy 1 1 Vaccination Vaccination strategy strategy 2 2 1500 1500 Vaccination in Patch 3 Vaccination in Patch 3 1000 1000 500 500 0 0 0 0 2 2 4 4 6 6 8 8 10 10 12 12 14 14 Time Time (weeks) (weeks) Figure Figure Figure21. 21. Vaccination 21.Vaccination Vaccination lawlaw law in patch ininpatch patch3 3for3Strategies forfor Strategies 1 and 1 1and Strategies and2.2. 2. Now, Now,we Now,wewill we willchange will change the thematrix change the KKsoso matrix matrix Kthat so it ittakes that that the it takes takes thefollowing the following following upper-triangular form: upper-triangular upper-triangular form: form: 10−10 1A   10− 2 2 −20.1A      0.1012 1A12 0.1A A.12 0.10A.13 13 13  − 2 −2 0.1A K =  K0=  00.6 ×010 2 K = 0 0.6 ×.610 ×−10 1A   0.10A.23 23−23 0  0 0.9 × 10−22 −2   0 0 0 0 ×10  0.90×.910    Symmetry 2019, 11, 430 31 of 42 Symmetry 2019, 11, x FOR PEER REVIEW 32 of 43 In this case, the first patch has available information of the second and third patches, the second patch hasInonly this information case, the firstofpatch has available the third patch whichinformation has only of the second and This self-information. thirdstructure patches, implies the second for patch haspatch, the first only information for instance,of thethe that third patch which vaccination lawhas only self-information. considers an amount of 10% Thisofstructure implies individuals for the first patch, for instance, that the vaccination law considers an coming into the patch from the second and third ones in order to calculate the total administeredamount of 10% of individuals coming into vaccination. It is the patch from important the the to notice second and third difference withones in order respect to calculate to the previous the total example, whereadministered all the vaccination. It is important to notice the difference with respect to the amount of travelling individuals (coming in and going out of the patch) is considered to calculate previous example, where the all the amount of travelling individuals (coming in and going out of the patch) is vaccination. The illness evolution is displayed in the various Figures 22–27. In particular, the evolution considered to calculate theinfectious of the vaccination. underThethese illness evolution is is circumstances displayed depictedinforthe various each patchFigures in Figure 22–27. In particular, 22, Figure 24, and the evolution Figure 26. On ofthethe infectious other hand, the under these circumstances vaccination generated by is depicted each one of forthe each patchisin strategies Figures 22, displayed for 24, each patch in Figure 23, Figure 25, and Figure 27. The main conclusions drawn before regarding the for and 26. On the other hand, the vaccination generated by each one of the strategies is displayed each effect patch in an of applying Figures 23, 25,vaccination appropriate and 27. The to main conclusions individuals as welldrawn as thosebefore related regarding the effect of to the comparison applying 1anand of Strategies appropriate 2 hold here vaccination too. However, to individuals in this caseasthe wellpeakas in those the related infectious to in thereduced comparison less of by applying vaccination than in the previous example. The main reason for this issue is that with less Strategies 1 and 2 hold here too. However, in this case the peak in the infectious in reduced the by newapplying vaccination control matrix, K, thethan in the number of previous administeredexample. Theismain vaccines much reason lowerfor now this issue than is that in the with the previous new control matrix, K, the number of administered vaccines is much lower case. This fact can be observed by comparing the Figures 17 and 23, Figures 19 and 25, and Figures 21 now than in the previous case. This fact can be observed by comparing the Figures 17 and 23, Figures 19 and 25, and Figures and 27. This result shows the importance of vaccination campaigns in order to control an epidemic 21 and 27. This result shows the importance of vaccination campaigns in order to control an outbreak in a patchy environment. epidemic outbreak in a patchy environment. Figure 22. Evolution of the infectious subpopulation within patch 1 under different vaccination Figure 22. Evolution of the infectious subpopulation within patch 1 under different vaccination strategies and upper-triangular matrix K. strategies and upper-triangular matrix K. Symmetry 11, 430 2019,2019, Symmetry 11, x FOR PEER REVIEW 32 of 3342 of 43 Symmetry 2019, 11, x FOR PEER REVIEW 33 of 43 600 600 Vaccination strategy 1 Vaccination strategy Vaccination 1 2 strategy Vaccination strategy 2 500 500 400 400 Vaccination in Patch 1 Vaccination in Patch 1 300 300 200 200 100 100 0 0 0 2 4 6 8 10 12 14 0 2 4 6 8 Time (weeks) 10 12 14 Time (weeks) Figure 23. Vaccination law in patch 1 for Strategies 1 and 2 with upper-triangular matrix K. Figure23. Figure 23.Vaccination Vaccinationlaw lawininpatch patch11for forStrategies Strategies11and and22with withupper-triangular upper-triangularmatrix matrixK.K. Figure 24. Evolution of the infectious subpopulation within patch 2 under different vaccination Figure 24. Evolution of the infectious subpopulation within patch 2 under different vaccination strategies Figure 24.and upper-triangular Evolution matrix K. subpopulation within patch 2 under different vaccination of the infectious strategies and upper-triangular matrix K. strategies and upper-triangular matrix K. Symmetry 11, 430 2019,2019, Symmetry 11, x FOR PEER REVIEW 33 of 3442of 43 Symmetry 2019, 11, x FOR PEER REVIEW 34 of 43 400 400 Vaccination strategy 1 Vaccination Vaccination strategy strategy 1 2 Vaccination strategy 2 300 300 Vaccination in Patch 2 Vaccination in Patch 2 200 200 100 100 0 0 2 4 6 8 10 12 14 0 0 2 4 6 Time (weeks) 8 10 12 14 Time (weeks) Figure 25. Vaccination law in patch 2 for Strategies 1 and 2 with upper-triangular matrix K. Figure 25.Vaccination Figure25. Vaccinationlaw lawininpatch patch22for forStrategies Strategies11and and22with withupper-triangular matrixK.K. upper-triangularmatrix Figure Evolution 26. 26. Figure of the Evolution infectious of the subpopulation infectious within subpopulation patch within 3 under patch different 3 under vaccination different vaccination strategies Figure and 26. upper-triangular Evolution of the matrix K. infectious subpopulation strategies and upper-triangular matrix K. within patch 3 under different vaccination strategies and upper-triangular matrix K. Symmetry 2019, 11, 430 34 of 42 Symmetry 2019, 11, x FOR PEER REVIEW 35 of 43 400 Vaccination strategy 1 Vaccination strategy 2 300 Vaccination in Patch 3 200 100 0 0 2 4 6 8 10 12 14 Time (weeks) Figure 27.Vaccination Figure27. Vaccinationlaw lawininpatch patch33for forStrategies Strategies1 1and and2 2with withupper-triangular matrixK.K. upper-triangularmatrix 6. Conclusions 6. Conclusions This paper has considered a SIR epidemic model in a patchy environment, each patch being This paper has considered a SIR epidemic model in a patchy environment, each patch being assumed to have its own health or medical centre. It has been assumed that there are potential travellers assumed to have its own health or medical centre. It has been assumed that there are potential coming into and leaving each patch which are interchanged with the remaining patches. It has been travellers coming into and leaving each patch which are interchanged with the remaining patches. It assumed that the vaccination controls are exerted at each community health centre while either the has been assumed that the vaccination controls are exerted at each community health centre while total information or a partial information of the total subpopulations, including the interchanging ones, either the total information or a partial information of the total subpopulations, including the is shared by all the set of health centres of the whole environment under study. In this way, vaccination interchanging ones, is shared by all the set of health centres of the whole environment under study. control laws involving constant terms and feedback information on the susceptible subpopulations In this way, vaccination control laws involving constant terms and feedback information on the have been proposed and discussed to be administrated at each health centre. In the cases that not all susceptible subpopulations have been proposed and discussed to be administrated at each health the information of the subpopulations distributions at other patches is known by the health centre centre. In the cases that not all the information of the subpopulations distributions at other patches is of each particular patch, the feedback vaccination rule is considered to have a decentralized nature. known by the health centre of each particular patch, the feedback vaccination rule is considered to Since there the control laws involved crossed gains to take into account or not (if such gains are have a decentralized nature. Since there the control laws involved crossed gains to take into account zeroed) the couplings between patches, the vaccination action can be of either a centralized or of a or not (if such gains are zeroed) the couplings between patches, the vaccination action can be of (totally or partially) decentralized nature. The paper has also investigated the existence, allocation either a centralized or of a (totally or partially) decentralized nature. The paper has also investigated (depending on the vaccination control gains) and uniqueness of the disease-free equilibrium point as the existence, allocation (depending on the vaccination control gains) and uniqueness of the well as the existence of at least an attainable and stable endemic equilibrium point. A formal analytic disease-free equilibrium point as well as the existence of at least an attainable and stable endemic characterization of the potential whole set of endemic equilibrium points has also being given based equilibrium point. A formal analytic characterization of the potential whole set of endemic on algebraic mathematical tools for the solvability of algebraic systems of equations. equilibrium points has also being given based on algebraic mathematical tools for the solvability of algebraic Author systems ofM.D.S., Contributions: equations. S.A.-Q. and R.N. conceived the whole model; M.D.S. performed the theoretical analysis snd the mathematical proofs as well as the main paper elaboration and writing ; A.I. conceived the experiments and performed Author Contributions: the simulations M.D.S., S.A.-Q. andand R.N.paper corrections conceived ; S.A.-Q. the whole corrected model; M.D.S.the whole paper performed the equations theoretical and collaborated in the conceptual framework; R.N. also contributed to the practical discussions analysis snd the mathematical proofs as well as the main paper elaboration and writing ; A.I. conceived related to the the given examples in accordance with the theoretical framework. experiments and performed the simulations and paper corrections ; S.A.-Q. corrected the whole paper equations This Funding: rsearch was in and collaborated funded by the Spanish the conceptual Government framework; R.N. also through Grants contributed DPI2015-64766-R to the and practical discussions DPI2016-77271-R (MINECO/FEDER, UE), and by UPV/EHU through Grant related to the given examples in accordance with the theoretical framework. PGC 17/33. Acknowledgments: The authors are grateful to the Spanish Government for Grants DPI2015-64766-R and Funding: This (MINECO/FEDER, DPI2016-77271-R rsearch was funded byand UE), thetoSpanish UPV/EHU Government for Grant PGC through 17/33.Grants DPI2015-64766-R and DPI2016-77271-R (MINECO/FEDER, UE), and by UPV/EHU through Grant PGC 17/33. Conflicts of Interest: The authors declare that they do not have any competing interests. Acknowledgments: The authors are grateful to the Spanish Government for Grants DPI2015-64766-R and DPI2016-77271-R (MINECO/FEDER, UE), and to UPV/EHU for Grant PGC 17/33. Symmetry 2019, 11, 430 35 of 42 Appendix A Proof of Theorem 3. Note that U is nonsingular since it has non-positive off-diagonal entries with the sum of all the entries per column being positive. Thus, (−U ) ∈ MEn×n is non-singular matrix with U −1  0 so that (−U ) is a stability matrix. Property (i) has been proved. On the other hand, note that the Jacobian matrix of the linearized system solution trajectory of the infectious subpopulations around the disease-free equilibrium point is ( F − U ) where (−U ) is a Metzler stability matrix from Property (i). Therefore, such a linearized system is globally exponentially stable if F = 0, that is if β i = 0 (fully absence of illness) ; ∀i ∈ n. Since the constraints β i = 0; ∀i ∈ n remove the quadratic terms from the model dynamics, it follows also that the stability is asymptotically global for the whole model. Property (ii) has been proved. Now, note that F − U = In − FU −1 (−U ) since (−U ) is a  Metzler stability matrix, then non-singular, from Property (i). If F = 0, F − U = −U is a stability matrix and it continuous to be a stability matrix from the continuity of its eigenvalues as functions of its entries for any F 0 such that R0 = ρ FU −1 < 1. Therefore, the disease-free equilibrium point is  locally asymptotically stable if R0 < 1. It has a critically stable eigenvalue for R0 = 1 and it is unstable if R0 > 1. Property (iii) has been proved. On the other hand, decompose U = Ud + Uod , where Ud is the diagonal part of U and Uod is its off-diagonal part. Since Ud and U are non-singular, one gets:   U = Ud + Uod = Ud In + Ud−1 Uod (A1) Note also that the matrix DSS is non-singular from Theorem 1 and it can be decomposed as the sum of its diagonal DSSd , which is also non-singular, and non-diagonal DSSod , parts to yield:   −1 DSS = DSSd + DSSod = DSSd In + DSSd DSSod (A2) so that   −1   −1 U −1 = In + Ud−1 Uod Ud−1 ; DSS −1 −1 = In + DSSd DSSod −1 DSSd (A3) Assume that kUod k2 < 1/kUd−1 k2 and k DSSod k2 < 1/k DSSd −1 k2 . Then, one gets from Banach´s Perturbation Lemma [13]: −1 kUd−1 k2 k DSSd k2 kU −1 k 2 ≤ −1 ; k DSS k2 ≤ (A4) 1 − kUd−1 k2 kUod k2 −1 1 − k DSSd k2 k DSSod k2 Then, by using Equations (3), (4)–(8) and Equation (11), since F is diagonal and U −T U −1 is symmetric, the reproduction number satisfies that: R0 = ρ FU −1 ≤ k FU −1 k2 ≤ k F k2 kU −1 k2 = ρ( F ) λmax (U −T U −1 ) = ρ( F )ρ1/2 U −T U −1  p  (A5) which leads to Equation (16). One gets also from Equation (16) and Equation (A4) that: −1   k DSSd k2  2 R0 = ρ FU −1 ≤ β max ( β ir ) −1 kΛ − V0 k2 k U −1 U −T k1/4 2 (A6) 1≤ i ≤ n 1 − k DSSd k2 k DSSod k2 which leads to Equation (17). Property (iv) has been proved. Finally, note from Equation (17) that R02 is minimized if Kij = aij ; ∀i, j ∈ n\{1}, implying that DSSod = 0 for any given model parameters and constant vaccination vector V0 . On the other hand, it follows from Equation (17) that R02 = 0, then R0 = 0, if Λi = Vi0 ; ∀i ∈ n is the influx of population into the i-th patch. Property (v) is proved.  Symmetry 2019, 11, 430 36 of 42 Appendix B Proof of Theorem 5. One firstly sums up the two first equations of Equation (1), so as to primarily delete the influence of the disease transmission rates towards a linearization study. Secondly, one expands the obtained result jointly with the third equation in a single compacted algebraic system while taking into account Equation (2). Then, one gets that Equations (22)–(23), subject to Equations (20), (21) and (24), hold. From Theorem 1 (i), the limit total population N ∗ is unique for the disease-free equilibrium point and any endemic attainable existing equilibrium point and this amount is allocated as first element in the linear system Equation (27). Then, one has to solve the auxiliary linear system Ax = b in x = x (y) with A and b defined in (26) which gives the endemic equilibrium points. It is known that there is (at least) one attainable endemic equilibrium point from Theorem 4 (ii) since R0 ≥ 1. Therefore, the above algebraic system has, at least, an attainable endemic solution and, from the Rouché–Froebenius theorem from Linear Algebra, Equation (25) holds. The whole set of endemic equilibrium solutions, including the attainable and unattainable ones, has to satisfy Equation (27). But note that the above algebraic system has only a partial information on the epidemic model Equations (1)–(2) since it does not include the information on the influence of the disease coefficient rates because of summing up action on the two first equations of Equation (1) leading to cancel the nonlinear common term. Therefore, the constraints Equation (28) are got by incorporating to Equation (27) the second equation of Equation (1) including the nonlinear term excluded from Equations (22) and (23). So, the particular vector y of the general solution Equation (27) is constrained to fulfill Equation (28). Property (i) has been proved. On the other hand, if B is irreducible, one deduces from Theorem 2 (i) that at any attainable endemic equilibrium point, the limit endemic infectious subpopulations at any patches are nonzero since if they are zero then there is no endemic infection. So, Property (ii) follows from the proof of Property (i) with y being restricted to belong to the set Ya . In the same way, Property (iii) follows with y restricted to belong to Yb since B is irreducible, A − K is irreducible and positive and Vi0 = Λi ; ∀i ∈ n so that the endemic equilibrium infectious population is positive at any patch and the susceptible ones at all patches are either all of them zero and or all of them nonzero from the first part of Theorem 2 (ii). Finally, Property (iv) follows under similar arguments from the second part of Theorem 2 (ii) involving the joint irreducibility of the positive matrices B, A − K and C.  Proof of Corollary 1. If R0 ≥ 1 then y ∈ R3n always exists such that x = A† b + I3n − A† A y  0 is an  endemic attainable equilibrium point from Theorem 4 (ii). Then, x = A† b + I3n − A† A y, where y =  y + y0 for any y0 ∈ Ker I3n − A† A . Since the whole endemic infectious subpopulation being the sum of  all the infection subpopulations in all the patches is non-zero, it holds that E y + A† (b − Ay)  0, that  is, the endemic infectious subpopulation in at least one patch has to be positive. If Bis irreducible, then the infectious subpopulations at the endemic steady-state are nonzero in all patches since, otherwise, the infections total endemic equilibrium subpopulation would be identically zero (Theorem 2.1 (i)). For uniqueness, of such an equilibrium point, the constraints Equation (28) should also hold (Theorem 5 (i)). The necessity of the constraints 1 to 3 for the uniqueness of any existing stable attainable endemic equilibrium point have been proved and it is also known that such a point always exists since R0 ≥ 1 (Theorem 4 (ii) and (iii)). Now, group the constraints Equation (28), as components of a vector β = ( β 1 , β 2 , · · · , β n ) T , resulting the following vector equation: β = γ1 γ2 (y)γ3 (y) = γ1 γ2 (y)γ3 (y) (A7) for any y = y + y0 with y0 ∈ Ker I3n − A† A , where  Symmetry 2019, 11, 430 37 of 42     T γ11 1/γ31 (y)  T γ12   1/γ32 (y)      γ1 =  ..  ; γ3 ( y ) =  ..  (A8) . .         T γ1n 1/γ3n (y) γ2 (y) = Diag[γ21 (y) , γ22 (y), ... , γ2n (y)] (A9)   n n T γ1i = diI + γi + ∑ b ji enT+i − ∑ bij enT+ j ; i ∈ n (A10) j(6=i )=1 j(6=i )=1   γ2i (y) = eiT A† b + 1 − eiT A† A y; i ∈ n (A11) h   ih   iT γ3i (y) = eiT A† b + I3n − A† A y A† b + I3n − A† A y en+i ; i ∈ n (A12) If the endemic equilibrium solution x is unique for y = y + y0 then y0 ∈ Ker I3n − A† A and  the given constant vector β of coefficient transmission rates satisfies Equation (28) for any y0 ∈ Ker I3n − A† A . If the constraint 3 is fulfilled for some y0 ∈ / Ker I3n − A† A then x is not unique and   Equation (28) is violated for y = y + y0 . Therefore,   the endemic equilibrium solution is unique under the constraints 1 to 3 if and only if ∆β = ∇yT β ∆y 6= 0 for the gradient matrix:  ∂β 1 ∂β 1 ∂β 1  ∂y1 ∂y2 ··· ∂y3n  ∂β 2 ∂β 2 ∂β 2   ...  ∇yT β =  ∂y1 ∂y2 ∂y3n  (A13)  ··· ... ... ...    ∂β n ∂β n ∂β n ∂y1 ∂y2 ··· ∂y3n for any ∆y ∈ / Ker I3n − A† A . In otherwords,and from the equivalence of a logic proposition with its  contra-positive one, if and only if, Ker ∇yT β ⊆ Ker I3n − A† A . Note that    γ21 (y)/γ31 (y) γ22 (y)/γ32 (y)     γ2 (y)γ3 (y) =  ..  (A14) .     γ2n (y)/γ3n (y) ^ Thus, in order to operate with the needed gradients in a closed form, define also the vector γ2 (y) ^ associated with the matrix γ2 (y) and the matrix γ3 (y) associated with the vector γ3 (y) as follows:   γ21 (y) γ(y)22   ^   γ2 ( y ) =  .. ; γ̂3 (y) = Diag[1/γ31 (y) , 1/γ32 (y), ... , 1/γ3n (y)] (A15) .     γ2n (y) Since the transposition and Moore–Penrose inversion can be permuted for any matrix, Equation (A12) can be expressed equivalently as follows: h T  T i h   i γ3i (y) = b T A† + y T I3n − A T A† ei enT+i A† b + I3n − A† A y ; i ∈ n (A16) Symmetry 2019, 11, 430 38 of 42 Note from Equations (A8)–(A12) via Equation (A14) subject to Equation (A15) and Equation (A16) that ^   ∇yT γ1 = 0; ∇yT γ2 (y) = Diag 1 − A† A 11 , 1 − A† A 22 , · · · , 1 − A† A 3n×3n ,    h T  T  i ∇yT γ3i (y) = 2 b T A† ei enT+i I3n − A† A − y T I3n − A T A† ei enT+i A† A ,   T  T   b T A† e1 enT+1 I3n − A† A − y T I3n − A T A† e1 enT+1 A† A    b T A †T e e T  T   (A17) T † I3n − A T A† e2 enT+2 A† A    2 n+2 I3n − A A − y  ∇ y T γ3 ( y ) = 2   ..     .     T † T T †  T T † T T † b A en e2n I3n − A A − y I3n − A A en e2n A A and direct gradient calculations yield:     ∇yT β ∆y = ∇yT [γ1 γ2 (y)γ3 (y)] ∆y = γ1 .∇yT [γ2 (y)γ3 (y)]∆y (A18)     ^ ^ = γ1 γ2 (y).∇yT γ3 (y) + ∇yT γ2 (y) γ3 (y) ∆y Then, the endemic equilibrium point is unique if and only if      ^ ^   Ker γ1 γ2 (y).∇yT γ3 (y) + ∇yT γ2 (y) γ3 (y) ⊆ Ker I3n+1 − A† A (A19) provided that the constraints 1–3 hold.  Proof of Theorem 6. For the endemic equilibrium point to exist and be attainable, there exists a ∗ = νI ∗ . If n = 1 the travel matrices in Equation (1) are non-negative real number ν such that Send end zeroed and one has at the endemic equilibrium point that: ∗ 2 ∗ Λ − βν Iend − dS νIend −V = 0 (A20)   ∗2 ∗ βν Iend − d I + γ Iend =0 (A21) ∗ γ Iend − d R R∗end + V = 0 (A22) ∗ ∗ 6 = 0, since ν = ∗ Send ∗ = d I +γ (d I +γ) Iend One gets from Equation (A21) for Iend ∗ Iend that Iend βν = ∗ βSend leading ∗ = d I +γ ∗ = β(Λ−V )−dS (d I +γ) to Send . Replacing this value in Equation (A20) leads to Iend . Note that β β(d I +γ) ∗ Send > 0 and also that if ∗ Iend ≥ 0, then v > 0 and ∗ ∗ ≥ 0 (respectively, Iend Iend > 0) if β ≥ β c (respectively, β > β c ). It is direct to see that the disease-free equilibrium point is Sd∗ f = Λd−SV , Id∗ f = 0 S∗ and R∗d f = dVR , and that β ≥ β c is fully equivalent to R0 = S∗d f ≥ 1 implying the attainability end ∗ V +γIend of the endemic equilibrium point. Note from Equation (A22) that R∗end = dR which leads to I + γ V + γ β ( Λ −V )− dS d I + γ β ( d ) [ ( ) ] R∗end = . βd R (d I +γ) After replacing the calculated endemic infectious amount. Note also that: (1) If R0 = 1 then the endemic equilibrium point is confluent with the disease-free one which is locally asymptotically stable. (2) If R0 < 1 then the endemic equilibrium point is not attainable since it has negative component. (3) If R0 < 1 then the disease- free equilibrium point is locally asymptotically stable since the state-solution trajectory of the Jacobian matrix at such a point is a stability matrix. It is also globally asymptotically stable since: (a) it is the unique attainable equilibrium point which is, Symmetry 2019, 11, 430 39 of 42 furthermore, locally asymptotically stable; (b) the total population is bounded; and (c) all the subpopulations are non-negative for all time implying that all of them are bounded for all time as result; (d) if it would be potentially surrounded by some limit cycle, such a cycle should be unstable since the critical point is asymptotically stable. Λ −V ∗ Λ−V0 −KSd∗ f Λ−V0 On the other hand, if V (t) = V0 + KS(t) then Sd∗ f = dS d f = dS leading to Sd∗ f = dS +K , V∗ V0 +KSd∗ f V +K (Λ−V0 )/(dS +K ) KΛ+dS V0 and R∗d f = ddRf = dR = 0 dR leading to R∗d f = . If V0 = K = 0 then d R (dS +K ) Sd∗ f = dΛS and R∗d f = 0 The result has been fully proved after calculating the total equilibrium population by summing up the susceptible and immune equilibrium subpopulations.  Appendix C −1  0, what is obvious since (−U ) is a Metzler Proof of Theorem 7. Note that there exists Uun un     −1    − 1 stability matrix. If ρ Uun U < 1then there exists Uun + U e e −1 U = Uun I3n + Uun e −1 =   −1 −1 U I3n + Uun e −1 . Thus, Uun        −1  h  i R0c = ρ Fc Uc−1 =ρ −1 e Fun + F Uun I3n + Uun U e −1 = ρ Fun + Fe M Uun (A23)   −1     −1 U where M = I3n + Uun e = Uc−1 Uun . Since M I3n + Uun −1 U −1 U e = I3n + Uun e M = I3n then   −1 −1 e −1 e −1 e −1 e M = I3n − Uun U M = I3n − Uun U I3n + Uun U = I3n − Uun UUun (A24) Thus, the following matrix equalities hold from:     Fc Uc−1 = Fun MUun −1 + FMU e −1 −1 e −1 un = Fun + F Uun I3n − UUc e (A25) Now, one has   −1   −1   −1 −1 e −1 −1 e I3n + Uun U = Uun Uun + Uun U = Uun + U e Uun = Uc−1 Uun (A26)     −1   −1 e −1 e −1 e −1 − Fun + Fe Uun U I3n + Uun U = − Fun + Fe Uun UUc Uun (A27) and one has from Equation (A27) that: R0c ≤ R0un if     −1 e c−1 + FU −1 e c−1 ≺ Fun Uun −1 0 ≺ Fun Uun I3n − UU e un I3n − UU (A28) Note that Equation (A28) is equivalent to:   −1 −1 e −1 −1 e c−1 ≺ 0 − Fun Uun ≺ − Fun Uun UUc + FU e un I3n − UU (A29) then to     −1 e c−1 ≺ FU −1 e c−1 ≺ Fun Uun −1 e −1 − Fun Uun I3n − UU e un I3n − UU UUc (A30) The constraints Equation (A30) can be written in equality form as follows:     −1 e c−1 + M1 = FU −1 −1 e −1 − Fun Uun I3n − UU e un I3n − U e = Fun Uun UUc − | M2 | (A31) Symmetry 2019, 11, 430 40 of 42   −1 −1 U for some given real 3n matrices M1  0 and M2  0. Since I3n + Uun e = Uc−1 Uun , note that h  i−1   −1 e Uc−1 Uun = U U −1 e Uun Uun + U e =U −1 e e I3n + Uun U (A32)     −1 −1 U so that, if ρ Uun −1 U e < 1, one has that I3n + Uun e exists and   −1   −1 e c−1 Uun = I3n − U I3n − UU e Uun + U e Uun = I3n − U −1 e e I3n + Uun U (A33)      −1  − is also nonsingular if ρ UUc e 1 e − 1 = ρ U I3n + Uun Ue < 1. From Banach Perturbation Lemma [13],   −1   −1 1 −1 e kU e Uun + U e Uun k = kU e I3n + Uun U k ≤ kU ek 2 −1 <1 (A34) 2 2 1 − kU k2 kUun e k2   that is, if kU e k < 1/2 kUun −1 k which ensures both the previous condition ρ U −1 U < 1 guaranteeing un e 2 2   −1   −1    −1  −1 that I3n + Uun −1 U e e c−1 Uun exists and that I3n − UU = I3n − U −1 U e I3n + Uun e −1 Uun exist. Thus, since , Equation (A31) is equivalent to    −1  −1 −1 + M − Fun Uun I − U I + U −1 U U −1 = FU −1 1 3n 3n un un e e e un    −1  −1    −1  −1 (A35) − 1 = Fun Uun UUc e − 1 − I3n − U I3n + Uun U e 1 e Uun− 1 − | M2 | I3n − U I3n + Uun U e − 1 e − Uun1   Recovering again the matrix inequality form for Equation (A35) and U e Uc−1 = U −1 U e I3n + Uun e , since M1 and M2 are arbitrary, yields that the condition 4 is equivalent to Equation (A35), which is also equivalent to Equation (A28), since Uun −1  0 if kU e k < 1/2 kUun−1 k . Property (i) has been proved. 2 2 Now, assume that Fe = − Fe ≺ 0. Then, Equation (A35) holds, and then Equation (A28) also holds, if,   for the given pair ( Fun , Uun ), the pair Fe = − F , Ue fulfils the matrix constraints:     −1  −1 −1 e −1 e −1 e −1 − Fun Uun U I3n + Uun U I3n − U e I3n + Uun U Uun Uun ≺ Fe ≺ Fun (A36) −1  0 and | F | U −1  0 then the matrix inequalities Equation (A36) imply that Property Since Fun Uun un (ii) holds since R0c ≤ R0un and, furthermore, R0c < R0un if either Fun Uun −1 or F e Uun−1 is irreducible. In the last case, one (but not both) of the symbols 00 ≺00 might be replaced with 00 ≺00 . This result is a direct application of Corollary 1.2 in [12] since if A and B are real matrices of the same order with A  B( 6= A)  0, equivalently, A  B  0 then the maximal eigenvalue of A is larger than that of B if A is irreducible but they can be identical if A is reducible.  References 1. Li, M.Y.; Shuai, Z. Global stability of an epidemic model in a patchy environment. Can. Appl. Math. Q. 2009, 17, 175–187. 2. Wang, W.; Zhao, X.Q. An epidemic model in a patchy environment. Math. Biosci. 2004, 190, 97–112. [CrossRef] [PubMed] 3. Muroya, Y.; Enatsu, Y.; Kuniya, Y. Global stability of extended multi-group SIR epidemic models with patches through migration and cross patch infection. Acta Math. Sci. 2013, 33, 341–3612. [CrossRef] Symmetry 2019, 11, 430 41 of 42 4. Iggidr, A.; Sallet, G.; Tsanou, B. Global stability analysis of a metapopulation SIS epidemic model. Math. Popul. Stud. 2012, 19, 115–129. [CrossRef] 5. Jin, Y.; Wang, W. The effect of population dispersal on the spread of a disease. J. Math. Anal. Appl. 2005, 308, 343–364. [CrossRef] 6. Sattenspiel, L.; Dietz, K. A structured epidemic model incorporating geographic mobility among regions. Math. Biosci. 1995, 128, 71–91. 7. Takaguchi, T.; Lambiotte, R. Sufficient conditions of endemic threshold on metapopulation networks. arXiv 2015, arXiv:1410.5116v2. [CrossRef] 8. Chalub, F.A.C.C.; Costa, T.J.; Patricio, P. Migrations, vaccinations and epidemic control. arXiv 2017, arXiv:1712.07918v1. 9. Khaleghian, P. Decentralization and public services: The case of immunization. Soc. Sci. Med. 2004, 59, 163–183. [CrossRef] 10. Singh, M.G. Decentralised Control; North-Holland Systems and Control Series; North Holland Publishing Company: New York, NY, USA, 1981; Volume 1. 11. Berman, A.; Plemmons, R.J. Nonnegative Matrices in the Mathematical Sciences; Academic Press: New York, NY, USA, 1979. 12. Kaczorek, T. Positive 1D and 2D Systems; Communications and Control Engineering Series; Springer: London, UK, 2002. 13. Ortega, J.M. Numerical Analysis; Academic Press: New York, NY, USA, 1972. 14. de la Sen, M.; Agarwal, R.P.; Nistal, R.; Alonso-Quesada, S.; Ibeas, A. A switched multicontroller for an SEIADR epidemic model with monitored equilibrium points and supervised transients and vaccination costs. Adv. Differ. Equ. 2018, 2018, 390. [CrossRef] 15. Nistal, R.; de la Sen, M.; Alonso-Quesada, S.; Ibeas, A. On a new discrete SEIADR model with mixed controls: Study of its properties. Mathematics 2019, 7, 18. [CrossRef] 16. Alonso-Quesada, S.; de la Sen, M.; Nistal, R. On vaccination strategies for a SISV epidemic model guaranteeing the nonexistence of endemic solutions. Discr. Dyn. Nat. Soc. 2018, 2018, 9484121. [CrossRef] 17. Xia, W.; Kundu, S.; Maitra, S. Dynamics of a delayed SEIQ epidemic model. Adv. Differ. Equ. 2018, 2018, 36. [CrossRef] 18. Barambones, O.; Garrido, A.J.; Garrido, I. Robust speed estimation and control of an induction motor drive based on artificial neural networks. Int. J. Adapt. Control Signal Process. 2008, 22, 440–464. [CrossRef] 19. Bakule, L.; de la Sen, M. Decentralized stabilization of networked complex composite systems with nonlinear perturbations. In Proceedings of the 2009 International Conference on Control and Automation, Christchurch, New Zealand, 9–11 December 2009; Volumes 1–3, pp. 2272–2277. 20. Ibeas, A.; de la Sen, M. Robustly stable adaptive control of a tandem of master-slave robotic manipulators with force reflection by using a multiestimation scheme. IEEE Trans. Cybern. Part B-Cybern. 2006, 36, 1162–1179. [CrossRef] 21. Kiouach, D.; Sabbar, Y. Stability and threshold of a stochastic SIRS epidemic model with vertical transmission and transfer from infectious to susceptible individuals. Discr. Dyn. Nat. Soc. 2018, 2018. [CrossRef] 22. Lee, C.; Garbett, A.; Wilkinson, D.J. A network epidemic model for online commissioning data. Stat. Comput. 2018, 28, 891–904. [CrossRef] 23. Sabbar, Y.; Kiouach, D. Long-time behavior of stochastic SIQD epidemic model with intervention strategies. In Proceedings of the International Conference on Fixed Point Theory and Applications ICFPTA´18, Mohammedia, Morocco, 8 May 2018; pp. 133–136. 24. Shamsi, N.G.; Torabi, S.A.; Shakouri, H.G. An option contract for vaccine procurement using the SIR epidemic model. Eur. J. Oper. Res. 2018, 267, 1122–1140. [CrossRef] 25. Jia, N.; Ding, L.; Liu, Y.J.; Hu, P. Global stability and optimal control of epidemic spreading on multiplex networks with nonlinear mutual interaction. Phys. A Stat. Mech. Its Appl. 2018, 502, 93–105. [CrossRef] 26. Kiouach, D.; Boulaasair, L. Stationary distribution and dynamic behaviour of a stochastic SIVR epidemic model with imperfect vaccine. J. Appl. Math. 2018, 2018. [CrossRef] 27. Das, A.; Pal, M. A mathematical study of an imprecise SIR epidemic model treatment control. J. Appl. Math. Comput. 2018, 56, 477–500. [CrossRef] Symmetry 2019, 11, 430 42 of 42 28. Alonso-Quesada, S.; de la Sen, M.; Nistal, R. A state feedback vaccination strategy applied to a SISV model for avoiding endemic equilibrium points. roceedings of the 2018 15th International Conference on Control, Automation, Robotics and Vision (ICARCV), Singapore, 18–21 November 2018; pp. 466–473. 29. Brockmann, D.; Helbing, D. The hidden geometry of a complex, network-driven contagion phenomena. Science 2013, 342, 1337–1342. [CrossRef] [PubMed] 30. Pei, S.; Kandula, S.; Yang, W.; Shaman, J. Forecasting the spatial transmission of influenza in United States. Proc. Natl. Acad. Sci. USA 2018, 115, 2753–2757. [CrossRef] 31. Okongo, M.O. The local and global stability of the disease free equilibrium in a co infection model of HIV/AIDS, tubercolosis and malaria. IOSR J. Math. 2015, 11, 33–43. 32. Barnett, S. Matrices in Control Theory with Applications to Linear Programming; Van Nostrand Reinhold Company: London, UK, 1971. 33. Bellman, R. The stability of solutions of linear differential equations. Duke Math. J. 1943, 10, 643–647. [CrossRef] 34. van den Driessche, P. Reproduction numbers of infectious disease models. Infect. Dis. Model. 2017, 2, 288–303. [CrossRef] [PubMed] 35. Biggerstaff, M.; Cauchemez, S.; Reed, C.; Gambhir, M.; Finelli, L. Estimates of the reproduction number for seasonal, pandemic, and zoonotic influenza: A systematic review of the literature. BMC Infect. Dis. 2014, 14, 480. [CrossRef] [PubMed] 36. Magal, P.; Webb, G. The parameter identification problem for SIR epidemic models: Identifying unreported cases. J. Math. Biol. 2018, 77, 1629–1648. [CrossRef] © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).