key: cord-0688548-2226moqe authors: nan title: The Basic Reproduction Number as a Loop Gain Matrix date: 2021-02-02 journal: IEEE Control Syst Lett DOI: 10.1109/lcsys.2021.3056616 sha: c9838aea5947d26cd314a803a7e466e09a46726e doc_id: 688548 cord_uid: 2226moqe The COVID-19 pandemic and the disordered reactions of most governments made the importance of mathematical modelling and model-based predictions evident, even outside the scientific community. The basic reproduction number [Formula: see text] quickly entered the common jargon, as a concise but effective tool to communicate the spreading power of a disease and estimate, at least roughly, the possible outcomes of the epidemic. However, while [Formula: see text] is easily defined for simple models, its proper definition is more subtle for larger, state-of-the-art models. Here we show that it is nothing else than the spectral radius of the gain matrix of a linear system, and that this matrix generalizes [Formula: see text] in the computation of the vector-valued final epidemic size and epidemic threshold, in a large class of finite-dimensional SIR-like models. I N DECEMBER 2019 a wave of pneumonia cases spread through the city of Wuhan (China) [1] , [2] . These were soon attributed to a new strain of coronavirus, SARS CoV-2, a virus related to the strains that caused the MERS and SARS epidemics in recent years [3] , but with different pathology and epidemiology that make it a relatively unpredictable pathogen. In a year the infection spread to affecting over 100 million people worldwide [4] . Chinese authorities reacted by setting up a massive campaign of testing and social restrictions. As the pathogen spread through other countries, governments were for the most part initially reticent to implement restrictive measures along the line of what was done in China, only to follow suit under the pressure of exponentially increasing infections and consequent fatalities. An unexpected by-product of the disease's disruptive dynamics was a surge in public consciousness about the importance of reliable modelling of the disease's dynamics, to provide predictions and guide the response to the pandemic. Common mathematical models, such as the SIR model by Kermack and McKendrick [5] , became common parlance, and the basic reproduction number R 0 entered the pages of daily newspapers. A finite-dimensional SIR model is a system of three differential equations describing the total number of susceptible (S), infected (I) and removed (R) individuals in a population. In virtually all studies on COVID-19 a simple variant of the model is used, which assumes no endemic equilibrium (all infected eventually heal or die), and no vital dynamics (birth and death for reasons other than the disease is neglected). Despite the simplicity of the setting, these SIR-like models can have hundreds of differential equations and thousands of parameters to model the geographic distribution or social stratification of the population [6] - [9] , and have been successfully used to model the diffusion of past epidemics [10] , [11] . In these cases, the time evolution of the S, I, and R compartments can be written as:Ṡ = −D(S)TI, where S ∈ R m + , I ∈ R n + , R ∈ R p + . Here and in the following we use the notation D(x) for a diagonal matrix with elements of vector x in the diagonal. Symbols A,Ã, B, and T instead denote matrices of suitable dimensions. We assume that (i) A is stable (all eigenvalues have strictly negative real part), (ii) A is Metzler (all elements except diagonal ones are nonnegative), and all elements ofÃ, B, C, and T are nonnegative, (iii) columns of B and C sum to 1, and columns of A +Ã sum to 0, (iv) all rows of T have at least one strictly positive element. Assumption (i) guarantees the absence of endemic equilibria. Indeed, if such an equilibrium existed, by definition it would satisfyṠ = 0 andİ = 0 for I = 0. The first equation requires D(S)TI = 0, which would imply that, in the second equation, AI = 0. This contradicts (i). Assumption (ii) follows from the nature of the right-hand side of (1), which describes flows between and within compartments, while assumption (iii) ensures mass conservation, i.e., S(t) 1 + I(t) 1 + R(t) 1 = 1, for all t. Assumptions (ii) and the structure of (1) further imply that the system is positive. Finally assumption (iv) ensures that all subcompartments of S can be infected. c IEEE 2021. This article is free to access and download, along with rights for full text and data mining, re-use and analysis. The matrix product BD(S)T in (1) is known in the epidemiological literature (e.g., [12] , [13] ) as the transmission matrix: it models all epidemiological events that lead to a new infection, through the flow D(S)TI from compartment S to I. Matrix A is instead known as the transition matrix: its elements encode the flows of individuals within compartment I, as well as from compartment I to R. The flow from I to R, in particular, is given byÃI, and is distributed among the subcompartments of R through matrix C. This family of compartmental models covers most of the finite dimensional, deterministic epidemic models that have been proposed and validated for the modelling of COVID-19, including socially stratified, spatially explicit (spatially structured) and networked models, though it excludes vital dynamics. As an illustrative example, consider the case of a toy metapopulation consisting of two SEIR subpopulations, where the exposed individuals (E) are not infective, while the infective individuals (I) can come into contact with susceptibles of either subpopulation. The model equations for one subpopulation arė with i = 1, j = 2, the other subpopulation model is obtained by swapping indices. The model has 8 state variables. To cast it into form (1) we collect them into three vectors: . The model matrices then become T := 0 a 1,1 0 a 1,2 0 a 2,1 0 a 2,2 , C := 1 1 0 0 0 0 1 1 , The sign of the elements of vectors and matrices plays an important role in the analysis of this system, so we need to introduce some notation to avoid confusion. Given a vector x or a matrix M, we write x ≥ 0 or M ≥ 0 if all elements of x or M are nonnegative. Such a vector or matrix is called nonnegative. In the case that all entries of x or M are strictly positive, we say that x or M is positive, writing x > 0 or M > 0. The contributions of this letter are as follows. In Section II, we rewrite system (1) to highlight the feedback structure of its I compartment's dynamics, and we use this structure to compute a loop gain matrix, called GD(S), which we show to share the same spectral properties as the Next Generation Matrix (NGM), the tool typically used to define the basic reproduction number in ODE-based epidemic models (see, e.g., [13] ). Our main contributions are in Theorems 2 and 3, which establish a relation between GD(S) and the asymptotic distribution of non-infected individuals among the compartments of S, and prove this distribution to be the globally asymptotically stable equilibrium of a discrete-time map over a compact set, and in Theorem 5, which proves that, for the class of systems we consider, a decrease of any term of matrix T (the contact rates) implies a decrease of the final epidemic size. These results suggest how GD(S) (more naturally than the NGM) plays the role of a matrix-valued basic reproduction number in the general class of systems described by (1) . We also show through Theorem 4 and its corollary how the epidemic threshold (the threshold of susceptible individuals, distributed among compartments of S, above which a new infection can trigger and epidemic) is related to the spectral radius of GD(S). In the light of the link between GD(S) and the Next Generation Matrix, which we established, this is a known result. The relevance stands in the link itself: if the epidemic threshold depends only on the dynamics of the linear feedback system (4), the vast body of results about the robust stability and control of linear feedback systems can be used to assess, for instance, the robustness of a sub-threshold population in uncertain models with nominal equations of the form (1). AS A LOOP GAIN MATRIX The basic reproduction number, R 0 , is defined as the expected number of new infections that a typical infected individual in (1) will cause, in a fully susceptible population. It was defined algebraically in [13] , [14] as the spectral radius of the NGM (or, more precisely, the NGM with large domain), which in the notation of (1) is written as where S 0 is the value of S at time 0. Note that the NGM is a nonnegative matrix, a fact that follows from B, D(S 0 ) and T being nonnegative, from our assumption of A being stable, and from the following result. Lemma 1 (See [15] ): Let A be a Metzler matrix whose eigenvalues have negative real part. Then A −1 ≤ 0. The Perron Frobenius theorem (see, e.g., [16] ) thus ensures that its spectral radius equals its largest positive eigenvalue. Starting from (3), we can give a further, more intuitive representation of R 0 as the spectral radius of the gain matrix of the linear system that models the instantaneous diffusion of the disease. We begin from the following algebraic result. Lemma 2: Let M and N be matrices such that MN and NM are square. Then the set of nonzero eigenvalues of MN is equal to the set of nonzero eigenvalues of NM. Now, we consider once again system (1), and specifically the equation ofİ. We can rewrite the dynamics of the infected as the following feedback looṗ The upper block is a linear time-invariant system. The input vector u, which has the same size as S, represents the (normalized) instantaneous flow of newly infected individuals, which are distributed through the compartment I as dictated by matrix B. The output vector y, of the same size as S and u, is the fraction of susceptible individuals in each subcompartment of S that come in contact with an infected, per unit time. The gain matrix of this linear system is and describes the ratio of y over u at steady state. The lower block is a time-dependent gain, simply multiplying the vector y by the square, time-variant matrix D(S), that is, translating the number y of infectious contacts into the flow u of newly infected individuals according to the current distribution of susceptible individuals. It follows that the loop gain matrix of the feedback loop for a fixed value of S is equal to If we evaluate it at S = S 0 , where S 0 encodes the structure of the metapopulation in the absence of the disease, then GD(S 0 ) indeed describes how a given distribution y of new infections in a fully susceptible population initially propagates under the dynamics of (1). It easily follows from Lemma 2 and the definition of the NGM that Theorem 1: R 0 is the spectral radius of the loop gain matrix GD(S 0 ). Proof: From Lemma 2, setting M := BD(S 0 ) and N := TA −1 . Note that, in the common case of a vector I of dimension larger than S (i.e., in metapopulations where each infected compartment is divided in more than one subcompartment), matrix GD(S 0 ) has smaller dimension than −BD(S 0 )TA −1 . It may resemble, in this sense, the NGM with small domain from [13] . The two matrices are however not identical, in general, even though they share the same set of nonzero eigenvalues. We see in the next section how the gain matrix GD(S 0 ) generalizes usage of R 0 in the computation of typical epidemic quantities for all models of the form (1). A typical use of R 0 is in estimating, based on initial conditions and the structure of the population, the final epidemic size, i.e., the total number of individuals that will be infected during the course of the epidemic. As long as I vanishes asymptotically, by mass conservation the asymptotic values of S or R, which we call S ∞ and R ∞ , can be used interchangeably for this purpose. This is a simple exercise in the case where S, I and R are scalar. First notice that, using (3), assuming S 0 = 1, and taking B = 1 according to our hypotheses, the basic reproduction number for the scalar system is taking the first two equations in (1) we can write Assuming that S 0 is asymptotic to 1 (therefore I 0 is asymptotic to 0) we obtain C = 1. Finally, assuming lim t→∞ I(t) = 0 we obtain We can generalize the above calculation as follows, using vector with i := S i,∞ /S i,0 to represent the final epidemic size weighted on the initial population distribution, that is, the fraction of individuals in each subcompartment of S that remain uninfected throughout the epidemic. Theorem 2: The asymptotic value satisfies where 1 is a vector with elements equal to 1 and ln( ) is the elementwise natural logarithm of . We know thatṠ = −D(s)TI and, by the assumption of no endemic equilibrium, that I ∞ = 0, therefore Multiplying both sides by TA −1 we obtain Finally, from the law ofṠ we have Substituting this in (7) and simplifying we obtain the desired relation. We can further prove existence and uniqueness of the solution of (6), and a means to compute it, as follows. Take x ∈ R m and define the map where e −x is the element-wise exponential of −x. Call X the box {x : 0 ≤ x ≤x} for somex > −TA −1 I 0 + GD(S 0 )1. The following theorem generalizes the result of [17, Th. 5.4 ] to the class of systems in (1), using a similar proof structure. Theorem 3: Defining := e −x and provided that I 0 > 0, the value of solution of (6) exists and is the unique fixed point of F(x) in X, and is the limit of the sequence x (k) := F(x (k−1) ), starting from any x (0) ∈ X. Proof: X is positively-invariant under F. Consider the sequences The map F is order-preserving, F(0) ≥ 0, and F(x) 0 (since which is a contradiction. Therefore F has a unique fixed point in X. In the statement of Theorem 2, I 0 appears as an additive term in the implicit definition of . This means that the distribution of individuals in S ∞ is minimally affected by the distribution of the initial seed I 0 of infected individuals: as I 0 1 → 0, S ∞ tends to a distribution independent of I 0 and implicitly given by the relation This is indeed the generalization of (5) to arbitrary models of the form (1), and GD(S 0 ) now plays the role of R 0 in the implicit expression of the final epidemic size. The above formula (or its more general version (6)) can be used, given an estimate of the initial population distribution S 0 , to evaluate the most effective among a set of policies affecting matrix G (e.g., distancing or isolation of given population classes), by solving an optimization problem (8) where Cost( ) is a suitable cost function. As a toy example, consider again the two-population model (2), with a 1,1 = 2.4, a 1,2 = 1.2, a 2,1 = 1, a 2,2 = 2, b 1 = b 2 = 2, c 1 = 1, c 2 = 1.2, S 1 (0) = 0.5, S 2 (0) = 0.5 (two populations of identical size, population 1 has slightly higher probability of contracting the disease from contacts with members of either population, and slightly lower recovery rate). Assume that a fixed amount of protective devices (e.g., face masks) are to be distributed among the two populations, and their effect is to reduce the transmission rates, proportionally to the amount of devices that are distributed. Assume that we aim to minimize the final epidemic size, so that Cost( ) := S 0 = S ∞ 1 . Calling α 1 and α 2 the effect of the distributed devices in the two populations, respectively, elements of T change as a i,j → a i,j (1 − α i ). Now, assuming that the total amount of available devices binds α 1 + α 2 = 0.4, we can construct the set G in (8) as G := {−T α 1 ,α 2 A −1 B : α 1 + α 2 = 0.4}, where −T α 1 ,α 2 denotes matrix T after rescaling of its elements. A nonlinear solver (in this case, MATLAB fmincon with default options) finds α 1 = 0.4, α 2 = 0 as the optimal resource allocation, corresponding to S ∞ 1 0.73. Note that an even resource allocation (α 1 = α 2 = 0.2) would have resulted in S ∞ 1 0.66: the optimal resource allocation, as opposed to an even allocation, can in this case spare about 7% of the population from contracting the disease. An epidemic threshold is a parameter value, in an epidemic model, below which the introduction of an infinitesimal number of infected individuals in an otherwise disease-free population cannot trigger an epidemic. In the scalar SIR model without vital dynamics (e.g., [18] ), this is given by When S 0 is below threshold, herd immunity prevents an epidemic outbreak. Mathematically speaking, this means that I 0 is an asymptotically stable equilibrium of (4) with S = S 0 . An alternative interpretation of the epidemic threshold in the same model is as the upper bound to the value of S ∞ . Once an epidemic has started, the set of individuals not affected by the disease at the end of the epidemic wave tends to a value below the epidemic threshold: In vector models such as (1) the fact that a given S be above or below the epidemic threshold depends not only on S 1 , but also on the relative value of its elements. Characterizing the set is therefore a bit more involved. We can generalize the relations (9) and (10) as follows. Call 1 the set of S vectors such that (4) with input S is asymptotically stable, and call 2 the set of S vectors such that there exists R ≥ 0 for which (S, I, R), with I = 0 is a limit point of at least one orbit of (1) with initial condition (S(0), I(0), R(0)) : I(0) > 0. The two sets coincide and are characterized as follows. Theorem 4: The proof of this theorem uses the following result Lemma 3 (See, e.g., Corollary 1.5, Chapter 2, in [19] ): Given a Metzler matrix M, its dominant eigenvalue λ pf is a weakly nondecreasing function of any positive perturbation of the elements of M Proof of Theorem 4: We consider the first two equations in (1) . The value of S is nonincreasing, and since S is bounded above 0 it must tend to some limit value S ∞ ≤ S(t), for all t ≥ 0. We can thus rewrite the equation of I aṡ which we can see as the linear systeṁ with a time-dependent perturbation d ≥ 0. The value of I instead must go to 0 as t → ∞ by the assumption of no endemic equilibrium. We now show by contradiction that (12) does not admit any solution I(t) → 0 with I(0) > 0, unless matrix (A + BD(S ∞ )T) is stable. First, note that (A + BD(S ∞ )T) is Metzler. Assume that its dominant eigenvalue is nonnegative, call v(0) a vector lying in its stable eigenspace, and assume that v(0) > 0. Then we can pick two finite perturbations ±p(0) such that v(0) ± p(0) ≥ 0 and such that p(0) is a linear combination of the eigenvectors of (A + BD(S ∞ )T) relative to null or positive eigenvalues. Consider how the flow of an initial condition v(0) ± p(0) evolves under the vector field (12) with d = 0. We have lim t→∞ v(t) ± p(t) = lim t→∞ ±p(t), with lim t→∞ |p(t)| > 0. Since ±p(t) cannot belong to the positive orthant at the same time, but the positive orthant of (12) is positively invariant, we conclude that the stable eigenspace of (A + BD(S ∞ )T) cannot contain any vector v > 0, unless the dominant eigenvalue of (A + BD(S ∞ )T) is negative. Now we take, as by assumption, an arbitrary I(0) > 0. By the above reasoning this vector must have nonzero component in the center or unstable eigenspaces of (A + BD(S ∞ )T), if they exist. This implies that if (A + BD(S ∞ )T) is not stable then the solution of the linear equation (12) cannot go to 0 for any d ≥ 0, and the same must hold for the solution of (11). We thus conclude that (A + BD(S ∞ )T) is stable at any limit point S ∞ . This means that 1 = 2 . We can now characterize the set . By the matrix determinant lemma, we have To prove that the converse also holds, note that by the above mentioned matrix determinant lemma det(I − GD(cS)) > 0, for all c ∈ [0, 1], implies that det(A + BD(cS)T) = 0, for all c ∈ [0, 1]. Since A+BD(cS)T is Metzler, the Perron Frobenius theorem ensures its dominant eigenvalue is real. We know that A + BD(cS)T is stable for c = 0, so it can only lose stability through a zero crossing of its real dominant eigenvalue, which would imply det(A + BD(cS)T) = 0. Therefore we have We have seen before that S ∞ must be such that A + BD(S)T is stable, so the above statement is equivalent to The result in Theorem 4, which generalizes the epidemic threshold, can also be restated as follows: This means that the dominant eigenvalue of GD(S) must be smaller than 1. As this is also the spectral radius by the Perron Frobenius theorem, we have ρ(GD(S)) < 1. The results above prove the almost global asymptotic stability of the states with S ∈ and I = 0, and characterize the set by the inequality ρ(GD(S)) = R 0 < 1. Equivalent characterizations have been given before albeit with different proofs see, e.g., [20] (Theorem 2). What makes the above formulation novel, besides the different proof, is that it elucidates the relation between these results and the classical theory of linear feedback systems: calling F the loop transfer function of (4) for a fixed S, ρ(GD(S)) ≤ F ∞ , where the latter denotes the H ∞ norm of the transfer function. Therefore, the stability of (4) by the Small Gain Theorem implies S ∈ . This means that questions related to S ∈ , such as the ability of a new disease to trigger an epidemic or the minimum number of infected at the end of an epidemic wave, when the nonlinear model (1) is subject to parametric or modelling uncertainties, can be attacked using linear robust stability techniques [21] . To conclude, we remark that in the first few months of the COVID-19 pandemic, when data was insufficient to reliably identify any predictive model, numerous publications stated that a reduction of social contacts (in (1), a reduction of the value of the elements of T) was nonetheless a fundamental tool to reduce the impact of the pandemic, as numerical simulations showed that it contributed to reducing the final epidemic size. While this is certainly a true and important message to convey to decision makers and to the public, its truth did not need the use of numerical simulations to be supported, as it holds for any, arbitrary SIR-like model of the form (1): Theorem 5: Let I 0 > 0, and let T jk be element (j, k) of matrix T. Then, dS ∞ dT jk ≤ 0 and dS ∞ dT jk = 0. Proof: We have seen that S ∞ and T are related through the vector-valued equation We have that dH(S ∞ ,T) . We know that S ∞ > 0, otherwise it would not satisfy H(S ∞ , T) = 0 due to the log term, so D −1 (S ∞ ) exists. Also, from Theorem 4 we know that (I − GD(S ∞ )) is nonsingular. This implies that ∂H(S ∞ ,T) Putting the expression for H(s ∞ , T) in the above formula we obtain dS ∞ dT jk = D(S ∞ )(I − GD(S ∞ )) −1 δ(j, k)A −1 I 0 where δ(j, k) is a matrix with all elements equal to zero, except element (j, k), which is equal to 1. Now, by assumption we have I 0 > 0, while by Lemma 1 A −1 ≤ 0. Since all other terms in the expression (δ(j, k)A −1 I 0 + δ(j, k)A −1 B(S 0 − S ∞ )) are nonnegative and A −1 is nonsingular, we have Furthermore, having A −1 ≤ 0 and nonsingular, and having I 0 > 0, imply that A −1 I 0 < 0, therefore We then know, from Corollary 1, that I − GD(S ∞ ) is stable and invertible, and we can easily verify that −(I −GD(S ∞ )) is Metzler, so by Lemma 1 we know that (I − GD(S ∞ )) −1 ≥ 0. Therefore, Nonsingularity of (I − GD(S ∞ )) −1 implies Finally, we have already seen that D(S ∞ ) ≥ 0 and is nonsingular, therefore dS ∞ dT jk ≤ 0, dS ∞ dT jk = 0. Taking a large class of SIR-like models, which covers many modern state-of-the-art networked, spatially explicit, and socially stratified models, we have shown how the basic reproduction number, commonly termed R 0 , is linked to the gain matrix of a linear system. This highlights a stronger connection than commonly acknowledged between a large class of epidemic models and the linear-feedback systems. We have shown how this matrix (and not just its spectral radius), plays the role of R 0 in the computation of the final epidemic size and the epidemic threshold, and how the epidemic threshold is essentially a simple by-product of the small gain theorem. Finally, we have formally proved that for all systems in the above-mentioned class, a reduction in the rate of contacts (at any time during the epidemic) reduces the final epidemic size. The COVID-19 epidemic Phase-adjusted estimation of the number of coronavirus disease 2019 cases in Wuhan, China COVID-19, SARS and MERS: Are they closely related? WHO Coronavirus Disease (COVID-19) Dashboard A contribution to the mathematical theory of epidemics Social contagion models on hypergraphs Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2) The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures Model for disease dynamics of a waterborne pathogen on a random network A systematic review of early modelling studies of Ebola virus disease in West Africa Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation (Wiley Series in Mathematical and Computational Biology) The construction of next-generation matrices for compartmental epidemic models On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations Sign properties of Metzler matrices with applications Positive Linear Systems: Theory and Applications On the dynamics of deterministic epidemic propagation over networks Unraveling R 0 : Considerations for public health applications Nonnegative Matrices in the Mathematical Sciences Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission