key: cord-0843654-t7yrjgvu authors: Huang, Yu-Jhe; Huang, Hsuan Te; Juang, Jonq; Wu, Cheng-Han title: Multistability of a Two-Dimensional Map Arising in an Influenza Model date: 2021-12-28 journal: J Nonlinear Sci DOI: 10.1007/s00332-021-09776-4 sha: ae0268c829a9fe7807576d63ffa721f8e5dda9b0 doc_id: 843654 cord_uid: t7yrjgvu In this paper, we propose and analyze a nonsmoothly two-dimensional map arising in a seasonal influenza model. Such map consists of both linear and nonlinear dynamics depending on where the map acts on its domain. The map exhibits a complicated and unpredictable dynamics such as fixed points, period points, chaotic attractors, or multistability depending on the ranges of a certain parameters. Surprisingly, bistable states include not only the coexistence of a stable fixed point and stable period three points but also that of stable period three points and a chaotic attractor. Among other things, we are able to prove rigorously the coexistence of the stable equilibrium and stable period three points for a certain range of the parameters. Our results also indicate that heterogeneity of the population drives the complication and unpredictability of the dynamics. Specifically, the most complex dynamics occur when the underlying basic reproduction number with respect to our model is an intermediate value and the large portion of the population in the same compartment changes in states the following season. The theory of dynamical systems is a major mathematical discipline closely intertwined with most of the main areas of mathematics. Furthermore, its concepts and methods can be applied to a wide variety of topics such as epidemiology, biology, medicine, physics, chemistry, finance, and more (Strogatz 2001; Jackson and Radunskaya 2015; Agarwal 1995; Swishchuk and Islam 2013) . The investigation of the dynamical systems has also revealed an interesting phenomenon, namely multistability behaviors or coexisting attractors (Natiq et al. 2018 Rahim et al. 2019; Martins and Gallas 2008; Casas and Rech 2012; Zhang and Luo 2018; Bao et al. 2021; Li et al. 2020a; Gilardi-Velázquez et al. 2017; Anzo-Hernández et al. 2018; Escalante-González and Campos 2020) . Multistability is the characteristic of a system presenting two or more mutually exclusive stable states. Bistable systems, for example, enable the implementation of logic gates (Santos-Moreno et al. 2020 ) and therefore computation. The applications of multistability can also be found in the field of image processing (Morfu et al. 2007 ) and morphing aircraft (Weisshaar 2013; Mattioni et al. 2008) . In this paper, we investigate a two-dimensional map that consists of both linear and nonlinear dynamics depending on where the map acts on its domain. It can be used to describe the dynamics of the season-to-season epidemic for an influenza. The derivation of the map is given in appendix. The two-dimensional map is also capable of generating very rich dynamics. For instance, such map displays a variety of qualitatively different dynamics: fixed points, period solutions, chaotic attractors, or bistable states. Surprisingly, bistable states include not only the coexistence of a stable fixed point and stable period three points but also that of stable period three points and a chaotic attractor. The goal of the paper is to classify the effect of the parameters on the dynamics of the model and to provide their theoretical analysis. Our theoretical main results contain the following. First, the existence, uniqueness, and stability of the fixed point are established. Moreover, we also show that if the adjusted basic reproductionr 0 is small, sayr 0 ≤ 1, then the disease dies out. On the other hand, ifr 0 is sufficiently large, then the epidemic returns every year. Second, the existence and nonexistence of a special type of period n points, epidemic returns once every n years, are obtained. Note that such special type, denoted by 3(0,0,+) type, of stable period three points may coexist with the nontrivial stable fixed point as well as with a chaotic attractor at certain range of parameters, respectively. We shall term, respectively, the above-mentioned phenomena as (1,3(0,0,+))-and (3(0,0,+),C)stability. Third, we prove the existence of (1,3(0,0,+))-stability over a wide range of parameters and also provide the possible range of parameters for producing the (3(0,0,+),C)-stability. It should be remarked that the verification of chaotic attractors is done by computing their Lyapunov exponents numerically. In summary, our theoretical and numerical results seem to suggest that great heterogeneity of the population contributes to the occurrence of the complication and unpredictability of the dynamics. Otherwise, the dynamics is much less complicated. We conclude the introductory section by mentioning the organization of the paper. Our main results are contained in Sect. 2. We also provide the bifurcation diagram of one of the state variables versusr 0 . Some interesting and complicated dynamics could be identified including but not limited to, other type of multistability, at the end of Sect. 2. The derivation of the model map is placed on the appendix. Some concluding remarks concerning the physical discussions and interpretations of the dynamics of the model as well as some related open problems are stated in Sect. 3. Of concern is a discrete iterative map defined in the following. where Here p n , x n , and y n are, respectively, the proportion of the infected population, the proportion of the quickly and slowly recovered populations at time n. For p n = 0 (resp., > 0), it is said that the epidemic dies out (resp., occurs) at time n.r 0 , a positive number, denotes the adjusted basic reproduction numbers. If 0 1, let Proposition 2.2 The map F defined in (1a) exists only one fixed point in T . In particular, ifr 0 > 1, then such fixed point exists in Tr 0 . Proof Ifr 0 ≤ 1, then the map F reduces to the linear map C, which has a unique fixed point in the origin. We next consider the case thatr 0 > 1. For such case, the only possible candidate for the fixed point of F lies in Tr 0 . Let (x, y, p) = (x * , y * , p * ), where (x * , y * ) is a fixed point of the map F. Then, x * = c 2 1−c 1 p * and y * = c 2 +c 3 1−c 1 p * . Let Upon using (1c), we have that p * satisfies the equation 1 = (1 − cp * )(r 0 −r 2 0 2 p * ). Solving the above equation, we get, for c = 0, Clearly, x * + y * (= cp * ) should be less than one for any feasible fixed point (x * , y * ) ∈ Tr 0 . However, cp + > 1 for any c ∈ [0, ∞). To see this, we compute ∂(cp + ) ∂c . The resulting calculation yields that cp + is strictly increasing in c. Since cp + | c=0 = 1, we have that cp + > 1. Thus, (x + , y + ) is not a feasible fixed point for the map defined in (1), where x + = c 2 1−c 1 p + and y + = c 2 +c 3 1−c 1 p + . Likewise, we define x − and y − similarly. To complete the proof, it remains to show that (x − , y − ) is a feasible fixed point. Some direct calculation would yield that cp − < 1 − 1 r 0 . Now, if c ≥r 0 2(r 0 − 1) or equivalently 2cr 0 ≥ 2c +r 0 , then p − < 1 as desired. If 0 ≤ c 1. For c > 1, we must have thatr 0 ≤ 2. Otherwise, c 1, then the inequality in (4) holds true as well. We just completed the proof of the proposition. Remark 2.1 (a) The fixed point mentioned in Proposition 2.2 has the following form x * = c 2 1−c 1 p − , and y * = c 2 +c 3 1−c 1 p − , where p − is defined in (3b). (b) Forr 0 ≤ 1, the map F reduces to the linear map C. Consequently, the epidemic dies out each season. We next seek to find the conditions for the existence of a special form of period n point. Specifically, a period n solution when there is an epidemic once in every n years is to be sought, i.e., the infected population sizes are all zero at all n seasons except one. For necessity, we set x k,m,n and y k,m,n to be the size of R φ (= x) and R 1 (= y), respectively, at the k-th season with the nonzero infected population size occurring at the m-th season, where 0 ≤ k, m ≤ n − 1. We further assume that its infected population size p m at the m-th season is equal top m,n . Should no ambiguity arise, we shall writep m,n =p n , x k,m,n = x k,m and y k,m,n = y k,m . Then, x 0,m y 0,m = C n x 0,m y 0,m +p n C n−m−1 c 2 c 3 , and so Now, Here Consequently, it follows from (5a) that Set m = 0. Then,p n satisfies the following equation Consequently,p Note that another root of equation (5e) is not feasible. Using (5a), we have that x i,0 = x 0,n−i and y i,0 = y 0,n−i i = 1, 2, . . . n − 1. Now, to verify the existence of the period n solution with p 0 =p n , p 1 = p 2 = · · · p n−1 = 0, it remains to show that and We have used (6b) to justify the first equality in (6d). The following proposition is to show that for fixed n, t n,n−i is decreasing in i. Hence, if (6d) is satisfied for i = n − 1, then (6d) holds true for all i = 1, 2, . . . , n − 1. For fixed n and 0 < c 1 < 1, the inequality t n,m > t n,m−1 holds true for m = 1, 2, . . . , n − 1. Proof Using (5d), we have that The above term is positive provided that c 1 (1−c n 1 )+n −m −1+(m +1)c n 1 −mc n+1 1 − (n − m)c 1 > 0. Rearranging the right hand side of the above inequality, we get that , which is, indeed, greater than zero whenever 0 < c 1 < 1. The proof of the proposition is complete. Clearly, (6c) is equivalent to Some direct calculation from the above inequality would yield that (6c) is always satisfied for any n as long asr 0 ≥ 1. For i = n − 1, we have, via Proposition 2.3, (5c) and (5d) that (6d) becomes For (7a) to be satisfied, its left-hand side has to be positive, which is equivalent to 1 r 0 > 2t n,0 − t n,1 2t n,0 (t n,1 + 1) . In particular, the above inequality is satisfied provided that 2t n,0 − t n,1 < 0 or 2t n,0 − t n,1 ≥ 0 and Suppose (7b) is satisfied. Then, (7a) can be further simplified to the following form. Consequently, the inequality in (7a) holds true provided that (7b) and is satisfied. In fact, only (7c) is needed to imply that (7a) holds true. To see this, it suffices to show thatr +,n ≤ k n whenever t n,0 < t n,1 < 2t n,0 , which can be verified by direct calculation. For (7c) to be feasible, it requires that r +,n > 1, or, equivalently, 2(t n,1 − t n,0 ) > 1. We are now in a position to state the conditions for the existence/nonexistence of such special type of period n points. (i) For n = 2, suppose Then, the corresponding (7d) is satisfied and the above-mentioned period two points exist whenever 1 1 for all 0 < c 1 < 1 and a 3 > 1 for 0 < c 1 < c e,3 , see Remark 2.2(i) for more details, we conclude that there are parameter regions for which these special types of period two or three points exist. We have just completed the proof of the assertions in (ii) and (iii). We next show that the assertion in (iii) holds true. Suppose a n < 1 for all n ≥ 4 and 0 ≤ c 1 ≤ 1. Note that b n < 1 for all n ≥ 3. Then, the equality in (7d) means that a n c 2 + b n c 3 > 1, which, in turn, implies that c 2 + c 3 > 1, for all n ≥ 4 and 0 ≤ c 1 ≤ 1, a contradiction to the assumption that c 2 + c 3 ≤ 1. Hence, to complete the third assertion of the theorem, it remains to show that (10) holds true. For n ≥ 5, it is easy to see that, by expanding n−1 k=0 c k 1 2 , there are no less than or equal to 2(n − 2) terms containing c k 1 , 0 ≤ k ≤ n − 3. Moreover, the total number of the terms containing c k 1 , n − 3 < k ≤ 2n − 3, is greater than 4. Hence, for n ≥ 5, (10) holds as claimed. For n = 4, it is obvious that We just proved the third assertion of the theorem. Now, (8a) and (8b) can be directly obtained from (9a) and (9b). For n = 2, clearly b 2 > 1 and a 2 < 1. Hence, for each 0 < c 1 < 1, there exist feasible c 2 and c 3 such that (8a) is satisfied. For n = 3, b 3 < 1, a 3 = 2(1+2c 3 1 ) 1+c 1 +c 2 1 . Some elementary calculations would yield that there exists a c e,3 such that a 3 > 1 (resp. < 1) whenever c ∈ (0, c e,3 ) (resp. (c e,3 , 1)). Hence, we have just completed the proof the theorem. Remark 2.2 (i) Numerically, c e,3 ≈ 0.3509168319. The notation c e,3 is so set up that its subscripts e and 3 represent the existence of period three points. Similar notations will be employed throughout the paper. The equations in (8a) and (8b) are to be denoted by e,2 and e,3 , respectively. For fixed c 1 , we illustrate the parameter region in the c 2 − c 3 plane for which the existence of period 2 or 3 points exists in Fig. 1 . (ii) The intersection point I 2,3 of e,2 and e,3 lies above the line (c 2 + c 3 = 1). Here . The verification of the above is cumbersome but elementary. This implies that the coexistence of such period two and three points is impossible. (iii) The lines and e,2 intersect at ( 1+c 1 2 , 1−c 1 2 ). (iv) Using (ii) and (iii), we also arrive at the conclusion that c 2 > c 3 is a necessary condition for the existence of such period three points. (v) Clearly, there are three types of period three points, depending on the number of the infected population being zero at the year end in a three-year span. We shall use the notations 3(+,+,+), 3(0,+,+), 3(0,0,+) to distinguish one from the other. In particular, the special type of period three points described in Proposition 2.4 is of 3(0,0,+) type. Similar notations are to be used for period n points if needed. (1+c 1 ) 2 c 2 + 2 (1+c 1 ) c 3 = 1 and e,3 : 2(1+2c 3 1 ) (1+c 1 +c 2 1 ) 2 c 2 + 2c 1 1+c 1 +c 2 1 c 3 = 1 with c 1 = 0.1. The green (resp., blue) region represents the parameter region for the existence of period two (resp., three) points mentioned in Proposition 2.4. These two regions do not intersect with each other (colour figure online) Our first main result deals with the stability of the fixed point. Forr 0 > 1, let (x * , y * , p * ), be defined as in Remark 2.1. The Jacobian matrix of F at (x * , y * ) has the following form. Here C is defined in (1b) and We have used the fact, via (1c), that p * = 2 to justify the last equality in (11a). Consequently, Some direct calculations would yield that the characteristic polynomial of Then, λ ± are complex (resp., real) roots whenever b ,1 < b * ,1 < 0 (resp., b * ,1 ≤ b ,1 ). We next compute the equivalent conditions for the following three inequalities and Since b * ,1 < 0 < 1−c 1 2c 2 +c 3 , the inequality λ + < 1 is satisfied for all feasible parameters. Hence, if b * ,1 ≤ b ,1 , then λ ± are real roots and so Then, (11e) is equivalent to b * ,1 > b c,1 . On the other hand, if c 1 ≥ c 2 2c 2 +c 3 , then (11e) is also satisfied and so is the inequality (λ + )(λ − ) < 1. We are now in a position to summarize the stability of the fixed point. (II) For 0 < c 1 < 1, we denote by A r ,c 1 (resp., A c,c 1 ) the interception point of the lines b r ,1 = −2 (resp., b c,1 = −2) and c 2 + c 3 = 1. Moreover, let I r ,c 1 (resp, I c,c 1 ) be the c 3 (resp., c 2 ) intercept of the line b r ,1 = −2 (resp., b c,1 = −2). Then, the (x * , y * ) is locally stable provided that the parameter point (c 2 , c 3 ) is in the interior of one of the following two regions. ( To see this, we get that b ,1 > b r ,1 if and only if g 1 (c 1 ) > 0. Indeed, The proof that condition (I-a) implies the stability of (x * , y * ) was provided on the paragraph in between (11e) and (11f). If c 1 ≥ c 2 2c 2 +c 3 , then b ,1 > 0, and so, b * ,1 < b ,1 due to the fact that b * ,1 < 0 as indicated in (11b). As a result, λ ± are real roots. It then follows from (I-a) that (x * , y * ) is stable as asserted. On the other hand, if c 1 < c 2 2c 2 +c 3 , then (x * , y * ) is stable as long as b * ,1 > b c,1 holds. Hence, if condition (Ic) holds, then (x * , y * ) is stable as claimed. We next show that if (I-d) is satisfied, then (x * , y * ) is stable. To see this, we break b ,1 into the following cases. Let b ,1 ≥ −2. It follows from (11a) that b * ,1 > max{b r ,1 , b c,1 }. Now, suppose b * ,1 > b ,1 . Then, (x * , y * ) is stable by (I-b) and (I-c). On the other hand, if b * ,1 ≤ b ,1 , then the stability of (x * , y * ) is guaranteed by (I-a). For the case that b ,1 < −2, we see that b * ,1 > max{b ,1 , b r ,1 , b c,1 }. Hence, (x * , y * ) is also stable by (I-b) and (I-c). We just completed the proof concerning (I-d). We now move to prove the second part of the theorem. In fact, if the algebraic condition (I-d) is converted into the geometric region in the (c 2 , c 3 ) plane, then the shape of corresponding region has two types, which are described in (II-a) and (II-b). To see this, we first note that b r ,1 ≤ −2 and b c,1 ≤ −2 are, respectively, equivalent to the following inequalities. It should be mentioned that the c 2 intercept of the line b c,1 = −2 is equal to one provided that c 1 = 2 − √ 3. It is then easy to see that the c 2 −intercept is greater than or equal to one (resp., smaller than one) provided that c 1 < 2 − √ 3 (resp., 2 − √ 3 < c 1 < 1 2 ). For c 1 ≥ 1 2 , b c,1 is always less than or equal to −2 since both the c 2 and c 3 intercepts are negative. We just completed the proof of the second part of the theorem. To prove the last part of the theorem, we see that lim r 0 →∞ b * ,1 = 0. It then follows from (11c) that λ ± are complex roots provided thatr 0 is sufficiently large. Upon using (11e), we get that |λ ± | < 1 wheneverr 0 is sufficiently large. where c is given in (3a). Some direct calculations would yield that ∂ ∂r 0 b * ,1 is increasing inr 0 and decreasing in c 1 . Moreover, it is easy to see that lim for all feasible c 1 andr 0 . This amounts to saying that one may chose suitable parameters so that b * ,1 can be arbitrarily close to −2 from the right. (b) The assumptions on the (II-a) and (II-a) of Theorem 2.1 are only sufficient conditions for (x * , y * ) to be stable. To see this, let c 1 = c 2 = 0,r 0 = 2 and c 3 > 1 2 , the corresponding b * ,1 , b ,1 and b r ,1 are, respectively, By choosing c 1 and c 2 sufficiently small and r 0 sufficiently close to 2, we have, via Theorem 2.1-(I-b), that the corresponding fixed point is locally stable. In the following, we state a sufficient condition for the stability of the special type of period two point described in Proposition 2.4. (8a) is satisfied and sor +,2 > 1. For 0 < c 1 < 1, if in addition, c 2 and c 3 satisfy the following inequality then the period two points of 2(0,+) type is locally stable for allr 0 ∈ (1,r +,2 ). Fig. 3 is where (c 2 , c 3 ) satisfies both (8a) and (12a). Proof For n = 2, let (x 0,0,2 , y 0,0,2 ,p 2 ) and (x 1,0,2 , y 1,0,2 , 0) be the period two points of 2(0,+) type. Then, the local stability of such period two points is determined by the size of the spectral radius of the product of two Jacobian matrices F (x 0,0,2 , y 0,0,2 ) and F (x 1,0,2 , y 1,0,2 ), which has the following form where b * ,2 = −2 r 2 0 (1−x 0,0,2 −y 0,0,2 ) 2 . Note that −2 ≤ b * ,2 ≤ −2 r 2 0 . The subscript 2 of b * ,2 indicates the quantity is associated with this special type of period two points. The corresponding characteristic polynomial of the matrix F 2 has the following form λ 2 − α 2 λ + β 2 = 0, where and Let 2 = α 2 2 − 4β 2 . If 2 < 0 and β 2 < 1 or 2 ≥ 0 and 1 + β 2 > |α 2 |, then such special type of period two points is locally stable. Since (8a) is satisfied, (c 2 , c 3 ) lies in the triangular region for which its three vertices are (0, 1), (0, 1+c 1 2 ) and ( 1+c 1 2 , 1−c 1 2 ). It is clear that the maximum of d 1 as defined in (12c) occurs at ( 1+c 1 2 , 1−c 1 2 ). Consequently, d 1 < 1 2 (1 − 2c 1 − c 2 1 ) for all c 1 ∈ (0, 1). Hence, β 2 < c 4 1 + c 2 1 (1 − 2c 1 − c 2 1 ) = c 2 1 (1 − 2c 1 ) < 1 regardless of the sign of 2 . The verification of the sign of 1+β 2 −α 2 is obvious since 1+β 2 −α 2 = (1−c 2 1 ) 2 −((1+c 2 1 −2c 3 1 )c 2 + (c 1 − c 3 1 )c 3 )b * ,2 , which is clearly positive for all 0 < c 1 < 1 regardless of the sign of 2 . To complete the proof of the theorem, it remains to show that 1 + β 2 + α 2 > 0 whenever (8a) and (12a) are satisfied. Let d 2 := (1 − c 2 1 + 2c 3 1 )c 2 + (c 1 + c 3 1 )c 3 . Then, 1 + β 2 + α 2 = (1 + c 2 1 ) 2 + d 2 b * ,2 , which is greater than zero whenever The above inequality is satisfied provided that b r ,2 ≤ −2, or equivalently, (12a) holds true. Clearly, the c 3 (resp., c 2 ) intercept of the line b r ,2 = −2 is greater (resp., less) than one. Moreover, the point I r ,c 1 = ( 1+c 1 2 , 1−c 1 2 ), which is the intersection of the line c 2 + c 3 = 1 and the line defined by the equality presented in (8a), and the origin are on the opposite side of the half planes, determined by b r ,2 = −2. Hence, the green region in Fig. 3 is a stable region. We have completed the proof of the theorem. Note that I r ,c 1 and the origin are on the opposite side of the half planes, divided by b r ,2 = −2. However, it should be pointed out that for the parameters chosen from the blue region, see Fig. 3 , the corresponding period two points are still possibly be stable depending on the range ofr 0 . We next investigate the stability of the period three points described in Proposition 2.4 and the possibility of bistable states. The numerical simulations suggest the that the existence of bistable states can be observed for certain ranges of the parameters. We then aim to provide inside as to what range of parameters would generate those bistable states. Consequently, the inequality b r ,3 ≤ −2 has the following form. 2(2c 1 − c 2 1 − c 4 1 + 2c 5 1 )c 2 + 2(c 2 1 + c 5 1 )c 3 ≤ (1 + 2c 3 1 + c 6 1 ). The stability region of the period three points of 3(0,0,+) type is addressed in the following. (I-a) There exists a c s 1 ,3 (< c e,3 ) such that if 0 < c 1 ≤ c s 1 ,3 , then the existence of such period three points implies its stability for anyr 0 ∈ (1,r +,3 ). Here c s 1 ,3 is the number so obtained that the c 2 -intercept of the line b r ,3 = −2 is equal to 1 when c 1 = c s 1 ,3 . Numerically, c s 1 ,3 ≈ 0.31949. A stable region in the c 2 -c 3 plane for such range of c 1 is depicted in the blue region in Fig. 1. (I-b) There exists a c s 2 ,3 such that, for c s 1 ,3 < c 1 < c s 2 ,3 andr 0 ∈ (1,r +,3 ), its corresponding stability region, determined by (8b) and (13b), is nonempty. Here c s 2 ,3 is the number so obtained that the three lines , see Fig. 1 , b r ,3 = −2, and e,3 have a common intersection when c 1 = c s 2 ,3 . Numerically, c s 2 ,3 ≈ 0.33747. A stability region placed on the c 2 -c 3 plane is illustrated as the green region in Fig. 6 . F exhibits (1,3(0,0,+) ) stability for anyr 0 ∈ (1,r +,3 ) if one of the following three conditions is fulfilled. (II-a) 0 < c 1 < 2 − √ 3, c 2 and c 3 satisfy (8b) and b c,1 ≤ −2, see the green region in Fig. 4 . c s 1 ,3 , c 2 and c 3 satisfy (8b), see the green region in Fig. 5 . (II-c) c s 1 ,3 < c 1 < c s 2 ,3 , c 2 and c 3 satisfy (8b) and (13b), see the green region in Fig. 6 . In particular, for parameters chosen from the region described in (II-b) , the existence of such period three points implies (1,3(0,0,+) )-stability. Proof For n = 3, let (x 0,0,3 , y 0,0,3 ,p 3 ), (x 1,0,3 , y 1,0,3 , 0) and (x 2,0,3 , y 2,0,3 , 0) be the period three points of 3(0,0,+) type. Then, the product of Jacobian matrices F at those three period three points, respectively, has the following form . Some direct calculation would yield that the characteristic polynomial of F 3 has the following form λ 2 − α 3 λ + β 3 = 0, where and β 3 = c 6 1 + c 4 1 (2c 1 − 1)c 2 + c 1 c 3 b * ,3 =: c 6 1 + c 4 1 d 4 b * ,3 We then follow the similar argument as those in proving the last theorem. To this end, we first show that |β 3 | < 1. Using Remark 2.2(i) and (iv), we have that −c 2 ≤ d 3 ≤ 1 2 c 2 . Hence, |β 3 | < c 6 1 + 2c 4 1 < 1 for c < 1 2 . Now, which is positive for any 0 < c 1 < 1. Next, we have that 1+β 3 +α 3 > 0 provided that b * ,3 > b r ,3 . The above inequality is satisfied as long as −2 ≥ b r ,3 or, equivalently, (13b) holds true. It is easy to check that the c 3 −intercept of the line b r ,3 = −2 is greater than 1 for 0 < c 1 ≤ c e,3 (< 1 2 ). The c 2 −intercept of the line b r ,3 = −2 is greater than or equal to 1 if and only if 1 − 4c 1 + 2c 2 1 + 2c 3 1 + 2c 4 1 − 4c 5 1 + c 6 1 =: g(c 1 ) > 0. It is easy to set that g is decreasing on (0, 1 3 ). Since g(0) > 0 , and g( 1 points implies its stability. To complete the assertion in (I-b), we need to show that the inequalities (8b), (13b) and 0 < c 2 + c 3 ≤ 1 have a nonempty (stability) region in the c 2 -c 3 plane. To this end, we first show that the slope m e of e,3 is larger than that m r of b r ,3 = −2 whenever c 1 < 0.4. Some direct calculation would yield that if and only if 1 + c 1 + c 2 1 − 5c 3 1 + c 4 1 + c 5 1 > 0, which clearly holds for all feasible c 1 . We next compute the intersection (c 2 ,c 3 ) of the lines 3 and b r ,3 = −2. After some tedious calculation, we have that 1 + 2c 5 1 − 10c 4 1 + 2c 3 1 + 2c 2 1 + 2c 1 and c 3 = 3c 8 1 + 4c 7 1 − 5c 6 1 + 3c 2 1 + 2c 1 − 1 2c 7 1 + 2c 6 1 − 10c 5 1 + 2c 4 1 + 2c 3 1 + 2c 2 1 and thatc 2 +c 3 ≤ 1 if and only if (c 2 1 − c 1 + 1)(2c 5 1 + 4c 4 1 − 4c 3 1 − 5c 2 1 − c 1 + 1) ≥ 0. Solvingc 2 +c 3 = 1 numerically, we have that c 1 ≈ 0.33747 =: c s 2 ,3 or −2.44783 or 1.25025. Clearly, c s 2 ,3 is the only feasible solution. Using all the above information, we arrive at the conclusion thatc 2 +c 3 < 1 (resp., ≥ 1) whenever c 1 ∈ (0, c s 2 ,3 ) (resp., [c s 2 ,3 , 1]). Upon using the fact that 0 > m e > m r , we further conclude that the stability region is nonempty (resp., empty) for c 1 ∈ (c s 1 ,3 , c s 2 ,3 ) (resp., [c s 2 ,3 , c e,3 )). For 0 < c 1 < 2 − √ 3, see Fig. 2b , we have, via Theorem 2.1-(II-b) and Theorem 2.3-(I), that the map F exhibits (1,3(0,0,+))-stability provided that assumptions in (IIa) hold true. Such region, see Fig. 4 , is nonempty provided that the c 2 −intercept of b c,1 = −2 is greater than that of 3 . We omit the elementary verification of such claims. Suppose 2 − √ 3 ≤ c 1 ≤ c s 1 ,3 . Then, the c 2 −coordinate of A r ,c 1 = 1+c 1 2 and that of the interception point A 3,c 1 of 3 and c 2 + c 3 = 1 is 1+c 2 1 +c 4 1 2(1−c 1 ) 2 (1+c 1 ) (> 1+c 1 2 ). Hence, A r ,c 1 is to the left of A 3,c 1 . It then follows from Theorem 2.1-(II-a) and Theorem 2.3-(I) that the map F exhibits the (1,3(0,0,+))-stability as long as such period three points exist. For c s 1 ,3 < c 1 < c s 2 ,3 , it is clear that the stability of period three points of 3(0,0,+) type implies (1,3(0,0,+))-stability of the map. We just completed the proof of the theorem. We conclude the section by mentioning that, for a certain range of the parameters, other type of multistability such as (3(0,0,+),C)-stability and (3(0,0,+),8)-stability can be observed numerically. Specifically, the parameters chosen from the blue region in Fig. 4 may exhibit interesting and complicated dynamics asr 0 varies. To see this, we first fix c 1 = 0.1, c 2 = 0.8 and c 3 = 0.2 and pick two sets of initial conditions. Their eventual states, see Fig. 7 , are then colored by red and blue, respectively. For 1 r +,3 . Hence, as compared to Fig. 7 , the existence of the (3(0,0,+),C) stability is no longer there. In fact, chaotic dynamics can only be observed wheneverr +,3 1 (resp., < 1), then Eq. (16) has a unique positive (resp., negative) solution. The season-to-season map is to be constructed in terms of two recovered classes R φ and R 1 with the assumption that those in R φ compartment have the stronger immunity than those in R 1 compartment. Specifically, we assume that c 0 and 1 − c 0 portions of those in R φ compartment at the beginning of the previous season remain in R φ compartment and move to R 1 compartment, respectively, at the beginning of the following season. Here 0 ≤ c 0 < 1. Furthermore, c 1 and 1−c 1 portions of those in R 1 compartment at the beginning of the previous season remain in R 1 compartment and move to S compartment, respectively, at the beginning of the following season. Here 0 ≤ c 1 < 1. Finally, c 2 , c 3 , and 1 − c 2 − c 3 portions of those in P compartment at the end of the previous season move to R φ , R 1 , and S at the beginning of the following year, respectively. Here 0 ≤ c 2 , c 3 < 1 and c 2 + c 3 ≤ 1. Moreover, denote by R φ n+1 and R 1 n+1 the size of R φ and R 1 , respectively, at time n + 1. We then arrive at the following equations. Herer 0 = r 0 1 + k . The newly defined termr 0 is to be called the adjusted basic reproduction number. Clearly, (17a), (17b) and (17c) define an algebraic-difference equation modeling the dynamics of a between-season influenza. We next seek to find an approximate solutionP n to equation (17c). To this end, we replace e −r 0 P n in (17c) by its Taylor polynomial 1 −r 0 P n +r 2 0 2 P 2 n of degree 2. Then, the approximate equation to equation (17c) reduces to (1 − R φ n − R 1 n )(r 0 −r 2 0 2 P n ) = 1. Denote the solution to the above approximate equation byP n . We then get the following approximate model. To save notations, we shall denote R φ , R 1 , andP by x, y, and p, respectively. Writing (18) into a compact form with the new notations, we then arrive at the iterative map defined in (1a)-(1c). The reasons to use such approximation are twofold. First, each flu season lasts only finite time. However, (17c) or (16) is derived from within-season model (14) as time approaches to infinity. It should also be noted that the size of the infected population is increasing in time. Therefore, in real situation, the size of the infected population at the end of each season should be smaller than the one produced from (16). It can be easily verified thatP n , as defined in (18c), is smaller than those obtained from (17c). Hence, it is more feasible to consider approximate model (1). Second, it is obvious that the map defined in (1) is easier to analyze analytically than that of defined in (17). Nevertheless, to analyze map F is still a nontrivial matter, because that the map F consists of both linear and nonlinear dynamics and which one of the types is to be acted on at time n + 1 depends upon the position of the previous iterate (x n , y n ). Specifically, if x n + y n ≤ 1 − 1 r 0 , then its corresponding the dynamics is nonlinear. Otherwise, the map F reduces to the linear map C, as given in (1b). Dynamical Systems and Applications On multistability behavior of unstable dissipative systems Initials-boosted coexisting Chaos in a 2-D sine map and its hardware implementation Multistability annihilation in the Hénon map through parameters modulation Mathematical Tools for Understanding Infectious Disease Dynamics Multistable systems with hidden and self-excited scroll attractors generated via piecewise linear systems Multistability in piecewise linear systems versus eigenspectra variation and round function Applications of Dynamical Systems in Biology and Medicine Two-dimensional memristive hyperchaotic maps and application in secure communication Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2) Superspreading and the effect of individual variation on disease emergence Multistability phase diagrams and statistical properties of the kicked rotor: a map with many coexisting attractors The application of thermally induced multistable composites to morphing aircraft structures On the use of multistability for image processing Self-excited and hidden attractors in a novel chaotic system with complicated multistability Can hyperchaotic maps with high complexity produce multistability? Dynamics of a new hyperchaotic system and multistability A simple influenza model with complicated dynamics Multistable and dynamic CRISPRi-based synthetic circuits Nonlinear Dynamics and Chaos: with Applications to Physics Random Dynamical Systems in Finance Morphing aircraft systems: historical perspectives and future challenges Multistability of a three-degree-of-freedom vibro-impact system Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations The following derivation of the map follows closely to the work done in Roberts et al. (2019) . To obtain our seasonal model, we begin with the derivation of the relationship between the concerned population groups within a season. Let the proportion of the population that is susceptible to infection at time t be S(t) and the infected proportion be I (t), to explore the effect of recovering from symptoms at different rates on an influenza epidemic as well as the effect of waning immunity. Such different rates of recovering may be due to the facts that people just recover from the flu differently or use non-drug approaches or even take wrong drug unintentionally or get two cooperative diseases, such as HIV virus and flu, in an overlapping period. We then further divide the proportion of the recovered population into two compartments, R φ (t) and R 1 (t) with the implication that the people in the R φ (t) compartment recover from the disease relatively quicker than those in R 1 (t) compartment. The population size is assumed to remain a constant, and so, S(t) + I (t) + R φ (t) + R 1 (t) = 1 for all t > 0. The motion of the dynamics then reads as follows.