Abstract

We considered a piecewise virus-immune dynamic model to investigate the effectiveness of the HIV virus loads-guided structured treatment interruptions (STIs). To better describe the biological reality, we extended the existing models by taking the carrying capacity of the virus loads into consideration to indicate the saturated growth of virus loads. We initially investigated the sliding dynamics of the proposed model and then obtained the global dynamics of the proposed model. Our main results showed that the system can exhibit very complex and diverse dynamic behaviors including a globally asymptotically stable equilibrium, bistable equilibra, and tristable equilibria, depending on the dynamics of the subsystems and the threshold level. In particular, an interesting result indicated that, with a proper threshold condition, the virus-guided therapy policy can successfully control the virus loads far below its carrying capacity and maintain the activity of the immune system for the case that the effector cells always go to zero without therapy or with continuous therapy. The finding suggested that the optimal strategy should be individual-based due to coexistence of multiple stable steady states, depending on the threshold conditions and the initial levels of viral loads and effector cells of the patients.

1. Introduction

Highly active antiretroviral therapy (HAART) plays significant roles in improving survival and reducing morbidity and mortality in HIV patients [1, 2]. Nevertheless, lifelong HAART confronts many complicated problems such as adherence difficulties and the evolution of drug resistance [36]. Structured treatment interruptions (STIs), the alternative strategies, have been designed to effectively achieve sustained specific immunity for early therapy in HIV infection. Indeed, STIs are beneficial for some chronically infected individuals who may need to take drugs throughout their lives as well as the reconstruction of patients’ immune system when they are not taking the drugs [7].

Since the first HIV patient received the structured antiretroviral therapy in 1999 [8], STIs have been widely investigated by many researchers. Basically, it includes two kinds of treatment regimes: carrying out the therapies at fixed moments or depending on the CD4 cell counts and/or viral loads. For example, some studies designed a scheduled treatment that therapy was initiated only when the CD4 cell counts decreased below 350 cells/μ [912]. Aiming to maintain CD4 cell counts higher than 350 cells/μ and HIV virus loads less than 100,000 copies/μ, Ruiz et al. designed an experiment to evaluate the safety of CD4 cell and HIV virus-guided STIs [13]. Many works have been done to compare the STIs with the continuous antiretroviral therapy (ART), but conflicting results have been reported [7, 11, 1317]. In particular, Maggiolo et al. made a conclusion that the structured treatment can result in the similar benefits for the patients with the continuous therapies [7], while the SMART study group found that STIs may also endanger the patients’ life if their immune responses are not sufficiently stimulated [18]. STI is still an important topic since lifelong ART brings great challenges to all HIV positive individuals. Moreover, as immune therapy arises, combining the ART with immune therapy (such as interleukin-2 (IL-2) treatment) may be a good choice for treating HIV patients [19]. Modelling and quantifying the efficacy of the combination of STIs and immune therapy remains unclear and falls within the scope of this study.

Mathematical models are useful tools for investigating the dynamics of HIV-1 infection [20, 21]. Although many studies modelled the continuous therapy for HIV-infected patients [2224], few attempts studied the STIs through mathematical modelling. In 2000, using a standard population model, Bonhoeffer et al. argued the benefits and risk generated from STIs [25]. In the study [26], Adams et al. analyzed the efficacy of STIs with fixed moments through the optimal control theory. Further, Rosenberg et al. pointed out that it is needed to include drug resistance and CD4 cell counts when making the structured treatment [27]. In 2006, Park et al. developed a mathematical model to investigate the impact of finite times of structured treatment guided by CD4 cell counts on HIV patients [28]. Then, Tang et al. [29] formulated a piecewise system to describe the CD4 cell-guided STIs, quantitatively explored STI strategies, and explained some controversial conclusions from different clinical studies. To further model immune therapy, Tang et al. [30] proposed a piecewise virus-immune dynamic model considering the combined antiretroviral therapy with interleukin-2 (IL-2) treatment. The model is given as follows:withwhere and represent the HIV virus loads and the effector cells, respectively. Note that a critical value of virus loads, denoted by , is set to trigger the combined therapy. That is, only when the virus loads exceed the threshold, we carry out the antiretroviral therapy and immune therapy (such as the interleukin-2 (IL-2) treatment) simultaneously for the patients, with being the elimination rate of HIV virus due to antiretroviral therapy and representing the growth rate of the effector cells due to immune therapy. Here, is the growth rate of HIV virus which incorporates both their multiplication and death, is the death rate of the effector cells, denotes the rate of binding of the effector cells to the virus loads, represents the rate of inactivation of the effector cells, and denotes that effector cell multiplication due to immune response has a maximum value as the HIV virus loads become large sufficiently.

In [30], the authors mainly discussed the global dynamics of system (1), and they further considered the efficacy of the effector cell-guided or both the effector cell and the HIV virus-guided piecewise therapies [31, 32]. Their results indicated that HIV virus can be controlled under some conditions with virus or effector cell-guided therapy, but still there are some parameter regions that the threshold policy can not prevent HIV from growing infinitely. As we can see from model (1), they simply assumed the linear growth of HIV virus. However, this is not what the thing looks like. Some clinical facts indicate that there should be “saturation effects” for the growth of HIV virus [33, 34]. The linear growth of HIV virus is unreasonable. Note that including the saturated growth of the HIV virus makes the model more natural but significantly changes the dynamics of the basic immune-virus system. Little is known about whether HIV virus-guided therapy based on this more natural description of its growth can completely control the growth of HIV and falls within the scope of this study.

The main purpose of this paper is to propose a piecewise model with the saturated growth for virus loads to describe HIV virus-guided therapy and then evaluate the effectiveness of the threshold policy through analyzing the global dynamics of the proposed model. The remaining parts of this paper are organized as follows. In Section 2, we initially propose the model, give some basic definitions of the Filippov system and also briefly conclude the dynamic behaviors of two subsystems. In Section 3, we mainly discuss the sliding mode and sliding dynamics, including the existence of the sliding domain and the pseudoequilibrium. In Section 4, we investigate the global dynamics of the proposed model. Finally, we provide the discussion and some biological conclusions.

2. Filippov Model and Preliminaries

In this paper, we extend model (1) by including the saturated growth of HIV virus, and hence we propose the following model:withHere, represents the carrying capacity of virus loads.

Letting , Then system (3) with (4) can be rewritten as the following general planer Filippov system [3537]:with and . Additionally, we describe the discontinuity boundary separating the two regions and as , where is a smooth scale function nonvanishing on . For convenience, we denote Filippov system (6) defined in region as system and that defined in region as system . Then we letwhere represents the standard scalar product and denotes the nonvanishing gradient of on . In the following, we replace with the notation for .

In order to consider the trajectory of the Filippov vector field , through a point , we distinguish the following regions on according to whether or not the vector field points towards it:

(a) Crossing region:

(b) Sliding region:

(c) Escaping region: Further, we can define a scalar differential equation on the sliding domain by employing the Filippov method. Letwith . Then, we have thatdescribes the dynamics of system (6) on the sliding segment . In what follows, we give the definitions of two types of equilibria of Filippov system (6).

Definition 1. A point is referred to as a real equilibrium of system (6) if , or . A point is called a virtual equilibrium of system (6) if , or . Both the real equilibrium and virtual equilibrium are called regular equilibria.

Definition 2. A point is referred to as a pseudoequilibrium if it is an equilibrium of the sliding mode of system (6); i.e., the equilibrium of the equation satisfies .

Before exploring the global dynamics of the Filippov system (3), it is essential to have a clear view of the dynamics of two basic subsystems. Thus, we first briefly make some conclusions on the dynamics of subsystem and subsystem .

For subsystem , there always exists a trivial equilibrium and a boundary equilibrium . The positive equilibrium of subsystem satisfies the following equations:Letand . It is easy to verify that if holds true, then has two positive roots which are given asParticularly, we know for and for or . Substituting or into the first equation of (13), we getThus, ifholds true, then ,  , are positive equilibria of system .

Next, we briefly investigate the stability of the equilibria of system . is locally stable if , which is equivalent to or , or , and . Define a function , and calculatewhere and are defined in system . Since all the parameters are positive and solutions to system initiating from are nonnegative, does not change its sign in the first quadrant. Hence according to Dulac-Bendixson criterion, we get that there are no limit cycles or homoclinic connections for system . Therefore, if there is no positive equilibrium, then the boundary equilibrium is globally asymptotically stable because system is bounded. Further, if and hold true, it is easy to verify that , and hence the boundary equilibrium is unstable. Therefore, under this scenario, positive equilibrium is globally asymptotically stable. When there are two positive equilibria of system , the dynamics of system have been shown in [38] that the equilibrium is always an unstable saddle and and are bistable. Then the dynamics of system can be concluded as follows.

Proposition 3. System always has a trivial equilibrium and a boundary equilibrium . If or and hold true, then there is no positive equilibrium, and the boundary equilibrium is globally asymptotically stable; If and , there is one and only one positive equilibrium which is globally asymptotically stable; If and , there are two positive equilibria and , and and are bistable.

If and , system has similar properties to system . Denote andSimilarly, ifholds true, then , , are positive equilibria of system . Then we can conclude the dynamics of system as follows.

Proposition 4. System has a trivial equilibrium and a boundary equilibrium . If or and , system has no positive equilibrium, and the boundary equilibrium is globally asymptotically stable; if and hold true, there is one and only one positive equilibrium of system which is globally asymptotically stable; if and hold true, there are two positive equilibria and , and both and are bistable.

3. Sliding Dynamics

According to the definition of in Section 2, we have Let , and one yieldsIt is reasonable to assume that the therapy should be done before virus growing and approaching its carrying capacity; therefore, in this study, we always assume that holds true. Then, if   , the sliding segment is ; if , then the sliding segment is .

Next, we employ Utkin’s equivalent control method introduced in [39] to examine the sliding dynamics in the region . It follows from and the first equation of system (3) thatSolving (23) with respect to yields Then substituting into the second equation of system (3) gives

Then, the sliding dynamics of system (3) can be described as (25). Let , and we obtain two rootsThis indicates that there are two possible pseudoequilibria and . In order to verify whether they are located at the sliding segment, we consider two situations and .

It is easy to see that can not be a pseudoequilibrium when , while it is always a pseudoequilibrium for . Thus, we only need to analyze the existence of . First, we prove the following lemmas.

Lemma 5. If and hold true, then we have .

Proof. LetTo ensure the function to be well defined, we need (i.e., ). Further, we know Definitely, we have , and hold true. Therefore, the domain for is , which is a continuous interval. Through simple calculation, we obtainThen, we haveBecause all the parameters are positive, always holds true in its domain. Then, according to , , and , we have . Similarly, letwhich is also defined in the interval . Note that is equivalent to . Further, it is obvious that . This means that always holds true in its domain. Thus, we obtain that because , and . In conclusion, we have if and . The proof is completed.

Lemma 6. If (or ) holds true, then we have   . If and , then we have that for or ; for or ; and for .

Proof. Let , and we know that if . Therefore, it follows from the formula of that when , then . If , we get for , and for or . This means that, when , we get for , and for or . Rearranging the formula of yieldsTherefore, applying similar analysis we have that when , then ; when , then for , and for or . This completes the proof.

Lemma 7. If and , then we have for .

Proof. Taking the derivative of with respect to givesLet , we have .
Thus, if , is increasing for and decreasing for ; if , then is decreasing for . Therefore, when , there may be three cases: (a) is always decreasing as increases, which means ; (b) is always increasing as increases, then ; (c) is initially increasing and then decreasing as increases, and hence . This completes the proof.

It is easy to verify that if , then we have , which means that there can not be when . According to the existence of possible equilibria of subsystems or , we consider the following cases to discuss the existence of the pseudoequilibrium .

Case  A1. There is no positive equilibrium for subsystem  . Then, we have three subcases:

(i) There is no positive equilibrium for subsystem   (denoted by Case  A11). Hence, there could beIt follows from Lemma 6 that for Case there is always for . This means that there is no pseudoequilibrium for . Then, we mainly consider . Now, there is . Therefore, we have for the first subcase of Case A11 (i.e., , ), which means is not feasible. For the other two subcases, if , is not feasible as well. However, if , then the sign of is not determined for the whole interval . We denoteTherefore, is feasible for , while it is not feasible for .

(ii) There is a unique equilibrium for subsystem   (denoted by Case A12); that is, For such case, there is no pseudoequilibrium for . There are two pseudoequilibria and for and , and there is only one pseudoequilibrium for or or and .

(iii) There are two equilibria for subsystem   (denoted by Case  A13); that is, Similar to Case A12, there is no pseudoequilibrium for . is feasible for , , or , .

Case  A2. There is one positive equilibrium for subsystem  . Then, consider the following:

(i) There is no positive equilibrium for subsystem   (denoted by Case  A21). Then, there should beIt follows from Lemma 6 that if   , then we have for . Further, because , according to Lemma 6, there is no pseudoequilibrium for while is feasible for . Similarly, for , the sign of can not be determined. Thus, when , there might be either two pseudoequilibria and for , or one pseudoequilibrium for . Especially, if , then is not feasible for .

(ii) There is a unique positive equilibrium for subsystem   (denoted by Case  A22); that is,For such case, it follows from Lemmas 6 and 7 that is feasible for . Especially, if , then there are two pseudoequilibria and for , and there is no pseudoequilibrium for .

(iii) There are two positive equilibria for subsystem   (denoted by Case  A23); that is, Here, we consider three subcases: , , or in terms of the relationship of to and . If , then is feasible for or or and . If , then is feasible for or and . If , then it follows from Lemma 7 that is feasible for or and .

Case A3. There are two positive equilibria for subsystem . Under this case, there can not be a unique positive equilibrium for subsystem according to Lemma 5. Then, consider the following:

(i) There is no positive equilibrium for subsystem   (denoted by Case  A31); that is, It follows from Lemma 6 that is feasible for .

(ii) There are two positive equilibria for subsystem   (denoted by Case  A32); that is,Similarly, according to Lemma 6, is a pseudoequilibrium for or .

Having discussed the existence of the pseudoequilibria, we take the derivative of with respect to to analyze the stability of the pseudoequilibria. Note that Then, we have and . Therefore, if holds true, then the pseudoequilibrium is locally stable on the sliding segment when it is feasible, and is unstable; if holds true, then becomes infeasible and the pseudoequilibrium is locally stable on the sliding segment when it exists.

4. Global Dynamics of System (3)

In this section, we will investigate the global dynamics of Filippov system (3). Initially, we shall rule out the existence of limit cycle for system (3). Here, we take Case A11 as an example to illustrate that there is not any limit cycles. The detailed proof is given in the appendix. Note that, for the other cases, we can prove the nonexistence of limit cycle for system (3) using the similar methods, and here we omit the detailed proof.

Lemma 8. For Case A11, there is no limit cycle for Filippov system (3).

Corresponding to the discussion of the sliding dynamics, we then consider the following cases to examine the global dynamics of system (3).

Case  A1. There is no positive equilibrium for subsystem  .

(i) There is no positive equilibrium for subsystem (i.e., Case A11). Then, we discuss the global dynamics of the Filippov system by considering the following subsituations:

(a) If , then we know that is a real equilibrium which is stable while is a virtual equilibrium which is unstable. As we mentioned above, there is no limit cycle for system (3), and moreover there is not any other stable equilibrium, thus the local stability of indicates that is globally asymptotically stable (shown in Figure 2(a)).

(b) If and hold true, then we have that becomes a virtual equilibrium. It follows from the discussion in Section 3 that the pseudoequilibrium is stable on the sliding segment . Similarly, system (3) does not have any limit cycle. Therefore, the local stability of indicates that is globally asymptotically stable (shown in Figure 2(b)).

(c) If and hold true, then there are two pseudoequilibria and on the sliding segment . As we proved in Section 3, is locally stable while is unstable. Similar to subcase (b), we have that is globally asymptotically stable (shown in Figure 2(c)).

Based on the above discussions, the global dynamics of system (3) can be concluded as follows.

Theorem 9. For Case A11, if  , then the boundary equilibrium is globally asymptotically stable; if , then the pseudoequilibrium is globally asymptotically stable for () being satisfied, while the pseudoequilibrium is globally asymptotically stable for () being satisfied.

(ii) There is a unique positive equilibrium for subsystem (i.e., Case A12). If , and are virtual equilibria, and is a real equilibrium which is stable. As mentioned above, under this scenario, there is no limit cycle or no pseudoequilibrium for Filippov system (3), and thus is globally stable. If , becomes a virtual equilibrium, and there emerges the stable pseudoequilibrium on the sliding segment . If , the same results hold compared with the case if holds true, while there are two pseudoequilibria (unstable) and (stable) on if holds true. Furthermore, If , becomes a real equilibrium which is locally stable. There is a pseudoequilibrium on which is unstable. Then, the global dynamics of Filippov system (3) are concluded as follows.

Theorem 10. For Case A12, if , then is globally asymptotically stable (shown in Figure 3(a)); if , then is globally asymptotically stable (shown in Figure 3(b)); if , then is globally asymptotically stable for () being satisfied and is globally asymptotically stable for () being satisfied (shown in Figure 3(c)); if , then is globally asymptotically stable (shown in Figure 3(d)).

(iii) There are two positive equilibria for subsystem (i.e., Case A13). If there is only one steady state (boundary equilibrium or positive equilibrium or pseudoequilibrium), then it is globally stable. We can analyze the dynamics of system (3) under the same scenario using the similar methods. However, there is a different situation for such case. That is, if holds true, then is a real equilibrium; hence it is locally stable. Furthermore, is a locally stable pseudoequilibrium if holds true, but if holds true, lose its stability and becomes a locally stable pseudoequilibrium. In conclusion, if , is bistable with for being satisfied. Based on the discussion above, the global dynamics of system (3) for Case can be concluded as follows.

Theorem 11. For Case A13, if  , then is globally asymptotically stable (shown in Figure 4(a)); if , then is globally asymptotically stable (shown in Figure 4(b)); if , then (or ) is globally asymptotically stable for ()  (or ()) being satisfied (shown in Figure 4(c) (or Figure 4(d))); if , then is globally asymptotically stable (shown in Figure 4(e)); if , then and are bistable for () being satisfied and and are bistable for () being satisfied (shown in Figure 4(f)). In particular, and are always bistable for .

Cases A2  (A21, A22, and A23) and A31. All the dynamics of system (3) for four cases are similar to the former three cases, and thus we just omit the proof here. The global dynamics of these cases are summarized in Table 1 with the dynamical behaviors shown in Figures 58, respectively.

Case  A32. There are two positive equilibria for both subsystems and , respectively. And the horizontal compartments of these equilibria satisfy . In terms of the existence and stability of the positive equilibria, the dynamics of the subsystems are similar to the focused situation in study [30]. Thus, based on the relationships of with respect to , , , and , we can analyze the dynamics of system (3) through the following subcases.

If the threshold condition is relatively small, i.e., , only , , and are real equilibria. It follows from the dynamics of system that both and are stable, and thus they are bistable for the Filippov system (3). If increases and satisfies , becomes a virtual equilibrium while the locally stable pseudoequilibrium emerges, and thus and are bistable. Further, we can conclude that is bistable with when or , while it is bistable with when . If , and are stable while there also exists pseudoequilibrium which is locally stable as well. Thus, based on the above discussions, we have the following conclusions.

Theorem 12. For Case A32, if , then equilibrium and the boundary equilibrium are bistable (shown in Figure 9(a)); if , then and are bistable (shown in Figure 9(b)); if , then and are bistable (shown in Figure 9(c)); if , then , , and are locally tristable (shown in Figure 9(d)); if , and are bistable (shown in Figure 9(e)); if , then and are bistable (shown in Figure 9(f)).

5. Conclusion and Discussion

Based on the basic immune-virus model, we considered the threshold policy depending on virus loads for the HIV-infected patients and obtained a nonsmooth system. Compared with the model in [30], we extend it by taking the saturated growth of virus into consideration, and hence our model is better in line with the biological reality. Meanwhile, our model also includes a piecewise elimination rate of HIV virus and growth rate of the effector cells to indicate that the antiretroviral therapy and interleukin-2 (IL-2) treatment are carried out simultaneously once the HIV virus loads increase above a threshold level.

We first briefly concluded the dynamics of two subsystems. Based on the properties of the subsystems, we discussed the sliding mode and sliding dynamics with different threshold conditions. In more detail, the sliding segments are , for , and for . We note that there is always a pseudoequilibrium on the sliding segment . And we further discussed the existence of the other pseudoequilibrium by considering eight different cases. Correspondingly, we then investigated the global dynamics of system (3) for various cases. Our main results show that the system can exhibit very rich dynamic behaviors: (a) one of the equilibria is globally asymptotically stable, which can be the pseudoequilibrium or the pseudoequilibrium or the boundary equilibrium of subsystem or the positive equilibrium or of subsystem or ; (b) two equilibria are bistable, particularly, which can be that the positive equilibrium is bistable with (or or ), or the boundary equilibrium is bistable with the pseudoequilibrium (or the positive equilibrium ); (c) three equilibria are locally tristable; that is, the pseudoequilibrium , the boundary equilibrium , and the positive equilibrium are stable for .

Global dynamics of our proposed model means that there is the case that the boundary equilibrium is globally stable for both the subsystems. Note that the immune-virus system without any therapy is actually system (the free system) while the immune-virus system with continuous therapy can be described by system (the control system). Comparing system with system , we can find that the continuous therapy can shrink the final density of the virus (i.e., ). It follows that the continuous therapy failed to maintain the activity of the immune system since the effector cell goes to zero finally. However, threshold policy indicated that, with a proper threshold condition, the STIs can not only successfully control virus loads to a previously given level which can be far below its natural carrying capacity (i.e., ) but also maintain the activity of the immune system (corresponding to the case that the pseudoequilibrium is globally stable), as Theorem 9 suggested.

Moreover, according to the dynamics of subsystem for Case we can see that virus loads will finally approach its carrying capacity if there is no therapy for the patients. With the continuous therapy, the virus loads of the patients can be successfully controlled far below its carrying capacity and their immune system can always maintain the activity; that is, equilibrium is globally asymptotically stable for system . For the virus-guided threshold policy, if the threshold level is set to satisfy , then the pseudoequilibrium is globally stable, and will tend to as . This means that the threshold policy can obtain the similar effects for the patients compared with the continuous therapy in terms of virus loads being controlled, which is in agreement with the clinical study [7].

It is worth mentioning that if we also assume that the two subsystems have two positive equilibria (i.e., Case A32), our model results, due to the saturated growth in model equations, show that the virus can not grow to infinity. Under this situation, we also found that the dynamics of system (3) are similar to those of the control system (system ) when the threshold condition is relatively small (i.e., ); that is, the positive equilibrium and the boundary equilibrium are bistable. And if the threshold is relatively high (i.e., ), the dynamics of system (3) are similar to the free system . Moreover, for middle values of the threshold (i.e., ), then equilibria , , and are locally tristable. Therefore, the bistability and the tristability suggest that the optimal strategy should also be individual-based by taking the initial conditions of the patients and the threshold conditions into consideration.

Appendix

Proof of Lemma 8

As mentioned in Section 2, we have precluded the existence of limit cycle for system or system . Then, we will first prove that no limit cycle exists containing part of the sliding segment . Thus, we consider the following situations:

(a) If , there is no pseudoequilibrium and the sliding dynamic equation (25) satisfies , so the trajectory on the sliding segment moves from the top to the bottom. Therefore, in the region in Figure 1, we need to show that any orbit starting from cannot hit again; if it does not hold, it contradicts the stability of of the subsystem . Therefore, no closed orbit containing part of the sliding segment exists.

(b) If and , we have that is locally asymptotically stable on the sliding segment . The local stability of implies that there is no limit cycle containing part of the sliding segment .

(c) If and hold true, we have that is locally asymptotically stable on the sliding segment . The local stability of on the sliding segment implies that there is no limit cycle containing part of the sliding segment .

Next, we preclude the limit cycle for Filippov system (3) which surrounds the sliding segment . Denote the right-hand side of (3) in the region (or ) by , where . We consider that C is a limit cycle surrounding the sliding segment . Let () be its part on the left (right) side of the line . Denoting that the limit cycle C and the line intersect at and , the intersection points of C and the two auxiliary lines and by , respectively, are shown in Figure 1, where    is any sufficiently small number. Denote the region delimited by () and the segment () by (), and the boundaries of and are denoted by and . Define the Dulac function by as before. From the discussions given above and Green’s theorem, we haveLet andsince , and, furthermore,and we getSuppose that the ordinates of the points are , where is continuous and satisfies , . Thus, we getSimilarly, we haveTherefore,which contradicts (A.4). Thus, we preclude the existence of the limit cycle C surrounding the sliding segment . The proof is completed.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This project was partially supported by the National Natural Science Foundation of China (NSFC, Grant nos. 11631012 (Yanni Xiao) and 11571273 (Yanni Xiao)) and by the National Mega-project of Science Research (Grant no. 2017ZX10201101- 002-002 (Yanni Xiao)). The authors thank Professor Robert Cheke for his valuable comments and for his generous help in polishing the manuscript.