A discrete two time scales model of a size-structured population of parasitized trees

1.   Introduction

In this work, we aim to build a basic modeling framework of the dynamics of parasitized trees. Although focused on a particular isolated, parasite/tree relationship, the study may have further importance due to the basic position that plants in general occupy in communities with various trophic levels ^[1].

By parasitized trees we understand those that are infected by one or more parasites, which leads to damage at different levels of their population dynamics. Our approach will initially be very general to end up analyzing certain characteristics of particular cases. The type of model we propose is purely based on population dynamics. We assume the parasite to be a macro-parasite that affects a population of trees structured by size.

We consider discrete-time models. In the literature on parasitized trees, there are both continuous-time models (see ^[2,3] or more recent examples like ^[4] or ^[5]) and discrete-time models ^[6]. Population dynamics models with nonoverlapping generations are usually represented as discrete-time models ^[7]. In this type of population, the processes act concentrated in small time intervals (egg laying, synchronized reproduction, etc.). Insects in our work that can play the role of parasites are typically included in this category. On the other hand, although in the dynamics of a tree population both time and size can be considered continuous variables, it is also common to consider discrete time and the population structured in a finite number of diameter classes ^[8]. This is done because the data held in forest inventories fits within that framework.

Some of the cited references ^[2,3,5] share a characteristic that we also include in our assumptions. These models are a better representation of real systems if the different timescales at which parasites and trees act are taken into account. Obviously, in our model the generations of parasites follow one another much faster than those of trees ^[9].

Systems of difference equations with two timescales have been proposed in ^[10,11]. The approach consists of choosing the time unit as the one associated with the slow process and assuming that the fast process acts $l$ times while the slow process acts only once. The parameter $l$ would represent the relationship between the timescales. Calling $S$ and $F$ the maps that represent the slow and fast processes, respectively, $X$ the vector with the state variables, and $t$ the time variable on the slow scale, the system would have the following expression: $X(t + 1) = S\big(F^l(X(t)\big)$. The superscript $l$ represents the $l$-th iterate of the corresponding map. In references ^[10,11], there are results that allow us to study some aspects of the asymptotic behavior of the solutions of these systems with the help of reduced systems. The basic idea is to let the fast process act first and, if it stabilizes, make the slow process act. The rigorous justification of this way of acting is developed in the above cited references.

The type of reduced system that we arrive at falls within the category of nonlinear matrix models ^[12,13]. These are of great use in the study of population dynamics, especially of the so-called structured populations. The parasite/tree models that we will initially propose cannot be considered matrix models. However, this will be the case of the associated reduced models through which their analysis will be carried out. The advantage of matrix models is that they have a well-developed theory that we exploit in our work. In particular, we use the characterization of population extinction or viability represented by the inherent net reproductive number. In ^[12], there is a general theorem that demonstrates the loss of stability of the extinction equilibrium when $n$ crosses a value of 1. This is completed with a general bifurcation result from the extinction equilibrium toward positive equilibria, which can be stable or unstable. In both ^[12] and ^[13], there are also results, complementary to the previous ones, related to the persistence and permanence of the model.

In the construction of our model we do not specify the type of parasite. We focus mainly on the negative effects they have on trees, which increases model generality. The dynamics of trees is defined, like the one of any size-structured population, by the specific survival, growth and fertility rates of each class. The vital rates that define the dynamics of the parasite are assumed to be dependent on the state vector of the tree population. Likewise, the negative influence that parasite density exerts on the trees, is reflected in their reduced vital rates.

In the model analysis we are able to find conditions for the parasite/tree community to be viable and, in that case, for it to tend to steady states or, alternatively, to persistently fluctuate. From these conditions, some conclusions are obtained about the influence of the parasite on the changes between the viability and extinction, and between the equilibrium and oscillation regimes. Evidently, the parasite can cause the tree population to go from viable to endangered. In this situation, distinguishing whether the parasite acts primarily on young or adult trees, differences are observed. By acting on young trees it is always possible that extinction is attained while this is not the case by acting on adults. In this same line, we see that acting on young trees favors the stabilization of the community. In circumstances where the non-parasitized population of trees fluctuates, the negative effect of the parasite on young trees can bring the community toward an equilibrium. On the contrary, the parasite acting mainly on adult trees can favor the destabilization of the tree population, in the sense of going from equilibrium when non-parasitized to oscillations when parasitized.

The organization of the rest of the work is as follows. Section 2 presents a general discrete-time model of a size-structured tree population living in a community with a macro-parasite. It is assumed that the dynamics of the parasite acts on a faster timescale than the one of the trees. Parasite dynamics are represented by the Beverton-Holt model with growth rate and carrying capacity broadly dependent on the tree population. Also in Section 2, we proceed to the construction of the reduced model that allows us to apply the results of ^[11]. Section 3 is devoted to specifying the properties of the different vital rates, both of parasites and trees, still in a very general way. Even so, results are obtained on the extinction/viability of the parasite/tree community. To obtain precise results of asymptotic behavior in situations of community viability in Section 4, the model is specialized to two size-classes for the trees and the forms of the different vital rates are specified. The results prove the existence of a transcritical bifurcation from the ecosystem vanishing, ${\bf 0}$) equilibrium, to positive and stable equilibrium points and, furthermore, for constant fertility rate, of a supercritical period-doubling bifurcation. In the Discussion, the results of the previous sections are interpreted, finding conditions for the long-term behavior of the tree population which changes whether in the presence or in the absence of the parasite. A brief conclusion presents some perspectives of the work, which ends with some appendices where the results stated in the previous sections are rigorously proved.

2.   General discrete-time model of a parasitized size-structured tree population

We consider a tree population divided into $q$ different size classes. We denote $n_j(t)$, $j\in\{1, \ldots, q\}$, the size of class $j$ at time $t = 0, 1, 2, \ldots$. The tree population vector at time $t$ is

The tree population is affected by a parasite that feeds on it. Let us denote $p(t)$ as the size of the parasite population at time $t$.

First, we give the basic form of the structured model. We use a general nonlinear Usher matrix model ^[12]. We call $\hat{s}_j\in (0, 1)$ the fraction of individuals of class $j$, $j\in\{1, \ldots, q\}$, alive at time $t$ that survive to time $t+1$, $\hat{m}_j\in (0, 1)$ the fraction of surviving individuals from class $j$, $j\in\{1, \ldots, q-1\}$, at time $t$ that grow up and belong to class $j+1$ at time $t+1$ and, finally, $\hat{f}_j\in [0, \infty)$, $j\in\{1, \ldots, q-1\}$, $\hat{f}_q > 0$, to the number of new recruits in class 1 at time $t+1$ for each individual in class $j\in\{1, \ldots, q\}$ at time $t$. We assume that the last class is fertile. In its most general formulation, all the parameters $\hat{s}_j$, $\hat{m}_j$ and $\hat{f}_j$ can be dependent on ${\bf n}$ and $p$.

The corresponding discrete-time system is the following one:

Calling ${\bf U}$ the $q\times q$ matrix in (2.1), and indicating its possible dependence on ${\bf n}(t)$ and $p(t)$, the system can be written for short in the following form

We assume that ${\bf U}\in C^1\big({\bf R}_+^{q+1}\big)$.

We could build the dynamics of the parasite from any of the discrete single-population models that exist in the literature ^[14,15,16]. This would give greater generality to the model and increase the difficulty of its analytical treatment. To cover the objectives of our work, it is enough to use the Beverton-Holt discrete-time model, whose dynamics are similar to the one of the continuous logistic model. In order for the parameters to maintain the usual interpretation of growth rate, $r$, and carrying capacity, $K$, we prefer to write the Beverton-Holt model ^[14] as follows:

To represent parasite-tree interactions in the general form, both parameters $K > 0$ and $r$ can be considered dependent on the population vector of the trees ${\bf n}$. We assume that $K, r\in C^1\big({\bf R}_+^{q}\big)$. The dependence on $p$ is already contemplated in the Beverton-Holt equation itself. For each value of ${\bf n}$, the map $g:{\bf R}_+^{q+1}\to {\bf R}_+$ that defines the passage of a unit of time in the parasite population dynamics is defined as follows

To include in the discrete-time model the dynamics of the tree and parasite population at different timescales we follow the procedure established in ^[10]. Roughly speaking, we consider the timescale of the joint model that is associated with the slow dynamics of the tree population and we assume that the fast dynamics of the parasites acts a number $l$ of times while that of the trees acts only once.

We represent by means of map $F:{\bf R}_+^{q+1}\to {\bf R}_+^{q+1}$ what happens in the fast timescale with the population of parasites assuming that the tree population is constant, and by means of a map $S:{\bf R}_+^{q+1}\to {\bf R}_+^{q+1}$ what happens in the slow timescale with the population of trees assuming constant parasite population.

Following ^[10], the discrete-time model with two timescales that we propose is

where $F^l$ is the $l$-th iterate of map $F$, that reflects the $l$ times the parasite dynamics acts during a slow time unit.

This type of system is generally difficult to study analytically. In ^[10], a general system reduction is proposed, which is further extended in ^[11]. In Appendix A, Theorem 2 in ^[11] is applied to justify the reduction and to detail the information on the complete model (2.5) that can be obtained from the reduced model (2.6).

The key point in the reduction is in the existence, for any constant population of trees ${\bf n}$, of a globally asymptotically stable (GAS) equilibrium point in the parasite equation

Assuming, henceforth, $r({\bf n}) > 0$, this equilibrium is $K({\bf n})$ ^[14], that is, nothing but the carrying capacity for the corresponding value of ${\bf n}$. Thus, for any initial population of parasites $p(0) > 0$, we have

The reduced system (2.6) associated with system (2.5) has ${\bf n}$ as the state variables vector, and it is obtained by substituting in the tree population dynamics system (2.1) the variable $p$ by the GAS equilibrium of the parasite equation $K({\bf n})$. As a consequence, the reduced system can be expressed as follows

The results in ^[11] allow us to extend, with better approximation the greater the ratio between timescales, some properties of the asymptotic behavior of the solutions of the reduced system (2.6) to those of the initial system (2.5). Fundamentally, they are those that have to do with belonging to the basin of attraction of a hyperbolic asymptotically stable equilibrium or periodic solution.

3.   Some general results on model (2.5)

In this section, we first detail the dependence on the population densities of the different parameters that appear in (2.5).

In the fast part of the system, the one that defines the dynamics of the parasite population, two parameters, $r$ and $K$, appear. We will consider the growth rate $r$ as constant. Its effect is that the population of parasites approaches the carrying capacity $K$ more or less quickly but, since we consider this dynamic faster than the one of the trees, whether it is a specific value or if it varies or not, it plays a superfluous role. On the other hand, it seems evident that the carrying capacity K has to be dependent on the size of the different classes of trees. The higher the density of trees, the greater the carrying capacity of the parasite.

The dynamics of the tree population are governed by three types of rates: survival, $\hat{s}_j$, maturation, $\hat{m}_j$, and fertility, $\hat{f}_j$. In the variation of these rates with respect to the densities of the populations, we want to include, on one hand, the parasites effect of feeding and, on the other one, the two types of competition that can occur between the trees. Trees compete for soil resources and also for light ^[17]. We assume that the level of competition for soil resources is a function of a weighted total population of trees $N_w = \sum_{j = 1}^{q}w_jn_j$, where $w_j > 0$ for $j = 1, \ldots, q$. Competition for light is asymmetric. Trees of a certain size are assumed to be affected only by those that are larger. Thus, for class $j$, the variable indicating the level of competition for light is defined as a weighted sum of trees in classes $i$ with $i > j$:

The weighting coefficients, $\bar{w}_{ji} > 0$, are considered dependent on both the class that is harmed by competition and the class that causes it.

Therefore, we assume that each of the three rates associated with each class of trees $j$ is a function of the three variables $p$, $N_w$, and $N^j$. In general, we only require that an increase in any of the three variables causes any of the three rates to decrease. Increased quantities of parasites or tree competition reduce survival, maturation/growth, and fertility.

The comments above lead to the following conditions on model parameters (note that the inequalities between vectors should be understood component-wise):

Parasite dynamics

Trees dynamics

Even though model (2.5) with hypotheses (3.1) and (3.2) that we have just established is still very general, we can already find some results about its asymptotic behavior. As we mentioned in Section 2, we can study the reduced system (2.6) to obtain information about the asymptotic behavior of the solutions of system (2.5).

System (2.6) is an example of a nonlinear matrix system ^[12,13]. A first question that we can study is the one related to the stability of the extinction equilibrium ${\bf 0}$, that is, finding conditions for it to be locally asymptotically stable (LAS) or even for it to be GAS. Obviously, ${\bf 0}$ is an equilibrium of system (2.6) and its local stability can be analyzed by linearization. The Jacobian matrix of the map $\bar{S}$ at ${\bf 0}$ is the matrix ${\bf U}\big({\bf 0}, K({\bf 0})\big){ = {\bf U}\big({\bf 0}, 0\big)}$, which is called the inherent projection matrix in the context of nonlinear matrix systems.

where $s_j = \hat{s}_j(0, 0, 0)$ and $f_j = \hat{f}_j(0, 0, 0)$ for $j = 1, \ldots, q$, and $m_j = \hat{m}_j(0, 0, 0)$ for $j = 1, \ldots, q-1$.

Equilibrium ${\bf 0}$ is LAS if the spectral radius of matrix $\bar{{\bf U}}_0$, $\rho(\bar{{\bf U}}_0)$, is less than 1. As the matrix $\bar{{\bf U}}_0$ is primitive, the Perron-Frobenius theorem ensures that $\rho(\bar{{\bf U}}_0)$ is an eigenvalue and it is strictly dominant.

In general, it is not possible to obtain an expression for $\rho(\bar{{\bf U}}_0)$ from the coefficients of $\bar{{\bf U}}_0$, which would allow us to have an explicit expression of the local stability condition of ${\bf 0}$. However, this can be done through the so-called net reproductive number $r_0$ associated with a certain decomposition of the matrix $\bar{{\bf U}}_0$. If we call ${\bf T}_0$ the matrix that contains only the terms of $\bar{{\bf U}}_0$ corresponding to survival and maturation, and $\textbf{F}_0$ the one that contains only those related to fertility, we have the following decomposition:

in which the rates of transition between size classes have been separated from those of recruitment of new individuals. The net reproductive number $r_0$ is defined as the spectral radius of matrix ${\bf F}_0({\bf I}-{\bf T}_0)^{-1}$, where ${\bf I}$ is the $q$-dimensional identity matrix. It is still an eigenvalue of a matrix, but easy to calculate since ${\bf F}_0$ has all its rows equal to zero except the first one. Furthemore, it is useful for our purposes in view of the fact that it satisfies the following equation ^[12],

Adapting the calculation of $r_0$ performed in ^[12] for the Usher model, we obtain

Now we can state a result of local stability of trivial equilibrium ${\bf 0}$ of system (2.6) and, with the help of an additional condition, also of global stability.

Proposition 1. Consider the nonlinear matrix system (2.6). Assume that conditions (3.1) and (3.2) hold.

(ⅰ) If $r_0 < 1$, then ${\bf 0}$ is LAS.

(ⅱ) If $r_0 < 1$ and, in addition, the maturation rates are constant, $\hat{m}_j = m_j$ ($j = 1, \ldots, q-1$), then ${\bf 0}$ is GAS, i.e., $\lim\limits_{t\to\infty}{\bf n}(t) = {\bf 0}$ for all ${\bf n}(0)\in {\bf R}_+^q$.

(ⅲ) If $r_0 > 1$, then ${\bf 0}$ is unstable.

Proof. It follows the proof of Th. 1.2.1 in ^[12], keeping in mind that the assumption on the $\hat{m}_j$ together with conditions (3.1) and (3.2) imply that $\bar{{\bf U}}_0\geq \bar{{\bf {\bf U}}}\big({\bf n}\big)$ for all ${\bf n}\in {\bf R}_+^q$.

As explained in detail in the Appendix A, equilibrium $({\bf 0}, 0)$ for system (2.5) has analogous properties to those established in Proposition 1 regarding equilibrium ${\bf 0}$ for system (2.6). In the following proposition, we specify it with the help of the Theorem 8 in the Appendix A.1.

Proposition 2. Consider the system (2.5). Assume that conditions (3.1) and (3.2) hold.

(ⅰ) If $r_0 < 1$, there exists $\delta_{0} > 0$ such that for every $\delta > 0$ satisfying $0 < \delta < \delta_{0},$ there exists $l_{\delta}\in {\bf N}$ such that for all $l\geq l_{\delta}$, $\bar{B}\big(({\bf 0}, 0)), \delta\big)$ is a trapping region for mapping $S(F^l)$.

(ⅱ) Assume that $r_0 < 1$ and that the maturation rates $\hat{m}_j$ are constant ($j = 1, \ldots, q-1$). Then, for any $({\bf n}(0), p(0))\in {\bf R}_+^q\times {\bf R}_+$ and every $\delta > 0$, there exist $t_{\delta}, l_{\delta}\in{\bf N}$ such that $\big(S(F^l)\big)^t ({\bf n}(0), p(0))\in B\big(({\bf 0}, 0), \delta\big)$ for all $l\geq l_{\delta}$ and $t\geq t_{\delta}$.

(ⅲ) Assume that $r_0 > 1$. Then, there exists $\delta_{0} > 0$ such that for every $\delta > 0$ satisfying $0 < \delta < \delta_{0}$ and any point $({\bf n}, p)\in {\bf R}_+^q\times {\bf R}_+$ such that $\left\Vert ({\bf n}, p)-({\bf 0}, 0)\right\Vert = \delta,$ there exists $l_{\delta}\in{\bf N}$ such that $\big(S(F^l)\big)({\bf n}, p)\notin \bar{B}\big(({\bf 0}, 0), \delta\big)$ for all $l\geq l_{\delta}$.

Roughly speaking, if $r_0 < 1$ and $\hat{m}_j$ are constant, any solution of (2.5) tends to be as close to $({\bf 0}, 0)$ as we like as long as the ratio between timescales $l$ is large enough. From expression (3.4), we obtain an extinction criterion, $r_0 < 1$, that is easily calculable in terms of the model parameters.

4.   Model (2.5) structured in two size classes

In this section we analyze a two size classes case of the general model (2.5) presented in Section 2. The assumptions we are making about the different rates involved in the model will allow us to apply the results developed in Section 3.

For the parasite dynamics, we assume that the growth rate $r$ is constant and that the carrying capacity $K$ is a linear combination of the densities of each of the two size classes of trees, $n_1$ and $n_2$, that we call from now on young and adult trees.

where $r$, $k_1$, and $k_2$ are positive constants. The latter represent the carrying capacity that each tree of the corresponding class contributes to the parasite population.

For tree dynamics, we must specify survival, maturation, and fertility rates.

The maturation rate is assumed to be constant $\hat{m}\big(n_1, n_2, p\big) = m\in (0, 1)$.

The survival rate of young trees, $\hat{s}_1\big(n_1, n_2, p\big)$, depends on their intrinsic survival rate, on parasites, and on asymmetric competition caused by adult trees. We denote by $s_1\in (0, 1)$ the intrinsic young trees survival rate. To take into account the influence of parasites and adult trees, we consider this survival rate as a decreasing function of both the number of parasites supported by each tree and the number of adult trees. Since each young tree represents a part $k_1$ of the parasite carrying capacity and each adult tree a part $k_2$, we assume that the fraction of parasites supported by each young tree is $k_1/(k_1 n_1+k_2 n_2)$ and, thus, the number of supported parasites per young tree is $pk_1/(k_1 n_1+k_2 n_2)$.

Further, we define

where $a_1$ and $b$ are positive parameters.

The survival rate of adult trees is defined analogously to that of young trees, eliminating asymmetric competition. It is the following decreasing function of the number of parasites supported by each one of them:

where $s_2\in (0, 1)$ is the adult tree intrinsic survival rate and $a_2$ is a positive parameter.

Finally, we assume that the young trees are not yet fertile and that the fertility of the adult trees is not affected by the parasites. On the other hand, adult fertility may depend on tree density, reflecting competition for soil resources. We consider it as a simple decreasing function of the number of young trees and adult trees:

where $c_1$ and $c_2$ are nonnegative parameters weighting the differences in resources consumption between young and adult trees.

The complete parasite-tree model with two timescales has the form of system (2.5)

where the map $F$, (2.4), is defined from the map $g(n_1, n_2, p)$:

and the map $S$ (2.4) from the $2\times 2$-matrix ${\bf U}(n_1, n_2, p)$:

For the three-dimensional system (4.3), the associated two-dimensional reduced system (2.6) is

with

where we have used as equilibrium for $p$ the value of $K({\bf n})$ defined in (4.1).

So, the reduced system reads as follows:

From now on, we analyze the system (4.5) in order to obtain information about the asymptotic behavior of the solutions of system (4.3). Basically, if the asymptotic behavior in (4.5) is defined by an equilibrium ${\bf n}^* = \big(n_1^*, n_2^*\big)$ (or an asymptotically stable periodic solution), then $\big({\bf n}^*, K({\bf n}^*)\big) = \big(n_1^*, n_2^*, k_1n_1^*+k_2n_2^*\big)$ (or the corresponding periodic solution) approximates the asymptotic behavior in (4.3). The larger the ratio of the timescales $l$, the better the approximation.

To begin, we calculate $r_0$ for system (4.5) and then apply Proposition 2 to system (4.3). In this case, the matrix $\bar{{\bf U}}_0$ (3.3) is

and its net reproductive number $r_0$ (3.4) is as follows:

Since system (4.5) is a particular case of system (2.6), Proposition 1 allows us to characterize the stability of equilibrium ${\bf 0}$ with the help of $r_0$. Furthermore, when $r_0 > 1$, we can establish the permanence of system (4.5), which roughly means that for positive solutions the total population neither explodes nor tends to zero (see Def. 1 in ^[18]).

Proposition 3. Consider system (4.5).

(ⅰ) If $r_0 < 1$, then ${\bf 0}$ is GAS.

(ⅱ) If $r_0 > 1$, then ${\bf 0}$ is unstable and (4.5) is permanent.

Proof. It is based on Theorems 1.2.1 of ^[12] and Theorem 3 of ^[18] (see Appendix A.1).

The global extinction condition of the population of both trees and parasites, $r_0 < 1$, can be expressed as:

When $r_0 > 1$, the system is permanent and we can turn to look for nontrivial asymptotic behaviors. To do this, we study the existence and stability of positive equilibria of the reduced system (4.5).

A first result that we can use is that of the bifurcation of positive equilibrium solutions of system (4.5), a nonlinear matrix system, when $r_0$ crosses the value 1. We follow ^[12] and use the inherent net reproductive number $r_0$ as the bifurcation parameter. We can make $r_0$ appear explicitly in the projection matrix $\bar{{\bf U}}\big({\bf n}\big)$ by decomposing it as the sum of the matrix of transitions $\bar{{\bf T}}\big({\bf n}\big)$ and the normalized fertility matrix $\bar{{\bf \Phi}}\big({\bf n}\big)$ times $r_0$:

where

and

Nonnegative or positive equilibria of system (4.5) satisfy the equilibrium equation

We are interested in the values of $r_0$ for which the Eq (4.9) has nonnegative solutions other than ${\bf 0}$. The following proposition ensures the existence of a supercritical bifurcation of positive and stable equilibrium points from the bifurcation critical value $r_0 = 1$, which is where ${\bf 0}$ loses its stability.

Proposition 4. Consider system (4.5).

(ⅰ) For each value of $r_0 > 1$, a positive equilibrium point ${\bf n}_{r_0}$ can be chosen such that the set of pairs $\mathcal{C}^+ = \{\big(r_0, {\bf n}_{r_0}\big):r_0\in [1, \infty) \}$, considering ${\bf n}_{1} = {\bf 0}$, is closed and connected, and the set of equilibrium points $\{{\bf n}_{r_0}:r_0\in [1, \infty) \}$ is unbounded.

(ⅱ) For every value of $r_0 < 1$, there exists no positive equilibrium point.

(ⅲ) There exists $\delta > 0$ such that for every $r_0\in (1, 1+\delta)$, the equilibrium point ${\bf n}_{r_0}$ is LAS.

Proof. It is based on Theorems 1.2.7 and 1.2.8 of ^[12] (see Appendix A.2).

The last proposition assures that when the inherent net reproductive number $r_0$ grows past 1, the stability of the equilibrium point ${\bf 0}$ disappears while a branch of stable positive equilibrium points appears (transcritical bifurcation), as shown in Figure 1. When we move away from $(1, {\bf 0})$ along $\mathcal{C}^+$, the positive equilibria may lose their stability. For this to happen, it is enough that the modulus of some eigenvalue of the Jacobian matrix of map $\bar{S}$ (4.5) at ${\bf n}_{r_0}$ become greater than 1 when we advance through $\mathcal{C}^+$. As it will be proven in the next subsection, in the example in Figure 1, the loss of stability occurs for values of $r_0 > r_{fb}\approx 1.367787$, at which a flip bifurcation occurs.

The results of Theorem 8 in Appendix A.1 allow us to extend the results of Proposition (4) to system (4.3), in the sense that the points $\big({\bf n}_{r_0}, K({\bf n}_{r_0})\big)$, for every $r_0\in (1, 1+\delta)$, act approximately as equilibria LAS in the asymptotic behavior of its solutions. Figure 2 compares the asymptotic behavior of solutions of system (4.3) with their approximations obtained from the reduced system (4.5). The comparison is established from the variable $p$ of (4.3), which is approximated by $k_1n_1+k_2n_2$, with $n_1$ and $n_2$ variables of system (4.5). In Figure 2a, the ratio between timescales has been taken to be $l = 5$, and $l = 10$ in Figure 2b. As can be seen, the approximation is good in both cases for the values of $r_0$ in which the asymptotic behavior is defined by an equilibrium. On the other hand, when 2-cycles appear, for $r_0 > r_{fb}$, we see that if $l = 10$, the approximation is still good, while it fails if $l = 5$.

In what follows, we try to specify the asymptotic behavior of the solutions of system (4.5) when $r_0$ grows beyond 1.

4.1.   Stability of positive equilibrium points of system (4.5) for constant fertility rate

The analysis in this section is performed with the added hypothesis that the fertility rate, (4.2), is constant, $\hat{f} (n_1, n_2) = f$. Thus, the reduced model (4.5) from which we start is the following:

The detailed analysis of this system is contained in Appendix A.3, where to perform the calculations the number of parameters of system (4.10) is reduced to five; see (A.11). Here, we only define a new parameter

which can be interpreted as the survival rate of adult trees taking into account the influence of parasites.

The net reproductive number (4.7) is expressed using $\sigma_2$ as

and, thus, condition (4.8) of population extinction can be expressed by the parameter $L$ (A.15):

When $f > L$, system (4.10) has a unique positive equilibrium $E^* = (n_1^*, n_2^*)$ (A.14) that can be written, using $\sigma_2$ and $L$, as

Just as $L$ gives us the existence condition of $E^* > {\bf 0}$, to express conditions of stability/instability of $E^*$ we need the parameter $M$ (A.21):

The results of the following proposition are included in Proposition 9 in Appendix A.3.

Proposition 5. Let $f > L$ hold in system (4.10), and $E^*$ be its positive equilibrium (4.12).

(ⅰ) If either $\sigma_2\geq 1/3$, or $\sigma_2 < 1/3$ and $f < M$, then $E^*$ is hyperbolic and LAS.

(ⅱ) If $\sigma_2 < 1/3$ and $f > M$, then $E^*$ is hyperbolic and unstable.

As we have noted in previous sections, using the results of Appendix A.1, the asymptotic results in Proposition 5 on the two-dimensional reduced system (4.10) admit an approximate extension to the corresponding three-dimensional starting system with two timescales (4.3).

As stated in the Figure 3, $E^*$ loses its stability when crossing the $f = M$ curve, and stable 2-cycles appear as asymptotic elements. In Appendix A.3, it is proven that system (4.10) undergoes a period-doubling bifurcation on $f = M$. We reproduce below Proposition 8 in the appendix for ease of reading.

Proposition 6. For $\sigma_2 < 1/3$, system (4.10) undergoes a period-doubling bifurcation at $f = M$.

Again, the results of Appendix A.1 allow us to make, although approximately, the same types of statements about system (4.3).

5.   Discussion

We have proposed a general model (2.5) of a tree size-structured population affected by a parasite. The model is a discrete-time one and it includes two timescales, a fast one for parasite dynamics and another slow one for the structured population dynamics. The model (2.5) is analyzed from the reduced model (2.6). The latter is a nonlinear population matrix model. For this type of model there is a well-developed theory ^[12] that allows us to study important features of the asymptotic behavior of its solutions.

Considering the model (2.5) in all its generality, Propositions 1 and 2 establish the conditions for the tree-parasite community to be viable or not. This is peformed with the help of the net reproductive number $r_0$, which can be calculated explicitly based on the model parameters (3.4). Thus, it is possible to analyze the influence of the different terms of the model on the viability of the population.

In order to achieve deeper results, in Section 4 the model is particularized by reducing the number of tree classes to two, young and adult, and proposing concrete expressions for the different growth, survival, maturation and fertility rates involved. The initial model becomes (4.3), which is studied with the help of the reduced model (4.5), in which the variable corresponding to the parasite does not appear and its asymptotic behavior is reflected in the new parameters.

As we mentioned before, the expression (4.7) of $r_0$ of system (4.5) allows us to decide about extinction or viability of the tree-parasite community. If $r_0 < 1$, we have extinction, and if $r_0 > 1$, no extinction. From this, the effect of the parasite on the viability of the tree population can be quantified. To do this, we write (4.7) using the $\sigma_2 = s_2/(1+a_2k_2)$ notation and, for the sake of symmetry between tree classes, also $\sigma_1 = s_1/(1+a_1k_1)$, which we had not used previously. $\sigma_1$ has an analogous interpretation as that of $\sigma_2$. It is a sort of survival rate of young trees taking into account the influence of parasites. The expression for $r_0$ becomes

This expression is strictly increasing with respect to any of the four variables on which it depends. In particular, it is in $\sigma_1$ and $\sigma_2$. As $\sigma_1\leq s_1$ and $\sigma_2\leq s_2$, we can ensure that

where $\bar{r}_0$ represents the net reproductive number of system (4.5) in the case of no parasites, $a_1 = a_2 = 0$. The condition that the tree population is viable without parasites and nonviable with parasites can then be expressed as $r_0 < 1 < \bar{r}_0$, or, isolating the fertility rate $f$, as

From this inequality, we see that, for any fertility rate $f$ large enough to allow viability of the tree population in the absence of parasites, sufficiently small values of $\sigma_1$ can be found that lead the parasitized population to extinction. A sufficiently large negative effect of parasites on the survival of young trees can be lethal to both tree and parasite populations. The same does not necessarily occur if the negative effect of the parasite on the survival of adult trees increases. As a function of $\sigma_2$, the last term of inequality (5.2), which is the parameter L in (4.11), is less than $\big(1-\sigma_1(1-m)\big)/(\sigma_1m)$, so if $f$ is greater, the parasitized population remains viable regardless of the greater or lesser damage that the parasite causes to adult trees.

We discussed the viability of the tree population in relation to the existence of the parasite. This discussion has revolved, in one way or another, around the relationship between $f$ and $L$ and the results in Propositions 3 and 4. We now proceed to carry out an analogous discussion on the influence of the parasite on whether the tree population balances toward a steady state or, on the contrary, oscillates toward a 2-cycle. This discussion is based on Propositions 5 and 6, and has as its axis the relationship between $f$ and the parameter $M$ (4.13), which, replacing $s_1/(1+a_1k_1)$ with $\sigma_1$, is expressed as follows:

From now on, we assume that the parasite effect preserves the viability of the tree population, i.e., $f > L$. We focus on finding situations in which the parasite changes the balance/oscillation tendency of the non-parasitized tree population.

Starting from (5.3), we define the function

so that the value of $M$ in the case of absence of parasite is $\mathcal{M}(s_1, s_2)$, and in the case of parasitized trees it is $\mathcal{M}(\sigma_1, \sigma_2)$. It is simple to verify that $\mathcal{M}$ is a strictly decreasing function in $x$ and strictly increasing in $y$. This allows us to easily address two particular cases, those in which the parasite affects only one of the two classes of trees. We will see that each one of them allows us to find conditions that ensure one of the two regimes shifts, stable equilibrium/oscillation, between parasitized and non-parasitized tree population that we are looking for.

Suppose that the parasite is harmful only to young trees and that the non-parasitized tree population tends to oscillate. This means that $a_2 = 0$ and, as a consequence of Proposition 5, that $s_2 < 1/3$ and $f > \mathcal{M}(s_1, s_2)$. From $a_2 = 0$ we get $\sigma_2 = s_2$ and, since $s_1 > \sigma_1$ and $\mathcal{M}$ is strictly decreasing in $x$, we have that $\mathcal{M}(s_1, s_2) < \mathcal{M}(\sigma_1, s_2) = \mathcal{M}(\sigma_1, \sigma_2)$, so the conditions for the parasite to make an oscillating population of trees dampen toward a stable equilibrium are, apart from $\sigma_2 = s_2 < 1/3$,

In situations where a high fertility rate coupled with a low survival rate of adult trees causes a non-parasitized population of trees to fluctuate, a parasite sufficiently damaging only young trees can balance the population. Arguments similar to those just presented allow us to conclude that a parasite affecting only young trees cannot change the trend of a population of non-parasitized trees from stable to oscillating.

If we now assume that the parasite only harms adult trees and maintain that the non-parasitized tree population tends to oscillate, we have $a_1 = 0$, $s_1 = \sigma_1$, $\sigma_2 < s_2 < 1/3$, and $f > \mathcal{M}(s_1, s_2)$. This implies $\sigma_2 < 1/3$ and, since $\mathcal{M}$ is strictly increasing in $y$, $\mathcal{M}(\sigma_1, \sigma_2) < \mathcal{M}(s_1, s_2) < f$, so Proposition 5 yields that the parasitized population also oscillates. On the other hand, if the non-parasitized tree population tends to a stable equilibrium, there are conditions that ensure that the parasite causes the population to oscillate. According to Proposition 5, without parasite the population settles to a stable equilibrium if either $s_2\geq 1/3$ or $s_2 < 1/3$ and $f < \mathcal{M}(s_1, s_2)$, and with a parasite there is an oscillation trend if $\sigma_2 < 1/3$ and $f > \mathcal{M}(\sigma_1, \sigma_2)$. If $a_2$ is large enough, we can obtain $\sigma_2 < 1/3$ and, provided that $f > (s_1(1-m)+1)/(s_1 m) = \lim_{\sigma_2\to 0^+}\mathcal{M}(\sigma_1, \sigma_2)$, also that $f > \mathcal{M}(\sigma_1, \sigma_2)$. If the parasite is sufficiently harmful ($\sigma_2$ small enough) and only to adult trees, and the fertility rate satisfies $\mathcal{M}(s_1, 0) < f < \mathcal{M}(s_1, s_2)$, it destabilizes a population that instead would settle to a stable equilibrium without parasite, inducing persistent oscillations.

6.   Conclusions

In this work we have presented a general framework for modeling a structured population that serves as support for another population, the latter exerting a negative influence on the former. We have talked about trees and parasites, but obviously the range of application is greater. We had in mind basic questions concerning plant diseases, or pests, including not only forest but also crop diseases. Plant diseases that, even though they share some characteristics of diseases in humans or animals ^[19], need to be treated specifically ^[1,20]. In particular, we have avoided the use of compartmental models of disease evolution.

We have managed to establish relationships between the intensity of damage caused by parasitization and different long-term population trends. We have taken advantage of the size structure to distinguish the effects of parasitism by class. It seems clear, however, that there is a long way to go to achieve realistic applications on this topic, which we plan to address in future publications hoping that the established framework is sufficiently broad and flexible.

Finally, we mention some aspects to include in the presented framework to get closer to real applications ^[20]. We think about applications that study plant diseases from the point of view of their economic impact and their control ^[19], including their effects on both biodiversity and ecosystem services ^[21]. The model used to represent parasite dynamics can be more dynamically rich, so that it better reflects changes in the host population, including those produced by harvesting or control measures themselves ^[22]. Cases in which the host population is harvested or is a deciduous plant that hosts a parasite feeding on its leaves must be modeled by including this periodic annual variation in its dynamics. Furthermore, its impact on the parasite must be able to be realistically represented through the dynamics proposed for it. In all cases, one objective of the model will be to allow theoretical comparison of the effectiveness of possible disease control methods.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Authors are supported by Ministerio de Ciencia e Innovación (Spain) (MCIN/AEI/ 10.13039/501100011033), Project PID2020-114814GB-I00.

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix

A.   Reduction of system (2.5)

The asymptotic behavior of the solutions of system (2.5) and the reduced system (2.6) can be related by making use of results in ^[11]. We need to find the limit $\bar{F}$ of the iterates $F^l$ of the function $F$ (2.4), and prove that the convergence of $F^l$ to $\bar{F}$ is uniform on compact sets.

Since ${\bf n}$ remains constant through $F$, to obtain $\lim_{l\rightarrow \infty}F^{l}$, it is enough to calculate the limit of the iterates of the function $g_{{\bf n}}(p) = g({\bf n}, p)$, which we know (Beverton-Holt model, ^[14]) it is $K({\bf n})$. Thus, we define

Lemma 7. Let $F$ and $\bar{F}$ be respectively defined in (2.4) and (A.1); the following limit exists uniformly on compact sets of $\Omega: = \big({\bf R}_+^q-\{{\bf 0}\}\big)\times (0, \infty)$

Proof. Let $H$ be any compact set included in $\Omega$ and choose any $\varepsilon > 0$. We need to prove that we can find $l_{0}\in {\bf N}$ such that

For any $({\bf n}', p')\in H$, since $K({\bf n}')$ is a GAS equilibrium of equation $p(t+1) = g({\bf n}', p(t))$, there exists $l_{({\bf n}', p')}\in {\bf N}$ such that

As the convergence of the solutions of equation $p(t+1) = g({\bf n}, p(t))$ to $K({\bf n})$ is monotonous, we also have, for any $({\bf n}, p)\in H$ and any $0\leq l_1\leq l_2$, that

Using (A.2) and the fact that $|K({\bf x})-g_{{\bf x}}^{l_{({\bf n}', p')}}(y)|$ is a continuous function for $({\bf x}, y)\in \Omega$, we can find an open neighborhood $U_{({\bf n}', p')}$ of ${({\bf n}', p')}$ such that

and using (A.3), we obtain that

To finish the proof, we perform a standard compactness argument.

$\{U_{({\bf n}', p')}\}_{({\bf n}', p')\in H}$ is an open covering of the compact set $H$, therefore, it possesses a finite sub-covering $\{U_{({\bf n}'_1, p'_1)}, \ldots, U_{({\bf n}'_m, p'_m)}\}$ and, choosing

we have

as we wanted to prove.

The next theorem, which relates the asymptotic behavior of systems (2.5) and (2.6) for big enough values of parameter $l$, is nothing more than Theorem 2 of ^[11] adapted to our setting.

Theorem 8. Let $m\in\mathbb{N\ }$ and let ${\bf n}^{\ast}\in\mathbb{R}_+^{q}$ be a hyperbolic $m$-periodic point of system (2.6).

(ⅰ) Assume ${\bf n}^{\ast}$ is LAS and let $({\bf n}(0), p(0))\in\Omega$ be such that $\lim_{t\rightarrow \infty}\bar{S}^{mt}({\bf n}(0)) = {\bf n}^{\ast}$. Then, it follows:

(1.a.) There exists $\delta_{0} > 0$ such that for every $\delta > 0$ satisfying $0 < \delta < \delta_{0},$ there exists $l_{\delta}\in {\bf N}$ such that for all $l\geq l_{\delta}$, $\bar{B}\big(({\bf n}^{\ast}, K({\bf n}^{\ast})), \delta\big)$ is a trapping region for mapping $\big(S(F^l)\big)^{m}$ and, moreover, $\big(S(F^l)\big)^{m}$ has at least a fixed point $({\bf n}_{l}^{\ast}, p_l^{\ast})\in B\big(({\bf n}^{\ast}, K({\bf n}^{\ast})), \delta\big)$.

(1.b.) There exists $\delta_{0} > 0$ such that, for every $\delta > 0$ satisfying $0 < \delta < \delta_{0},$ there exist $t_{\delta}, l_{\delta}\in{\bf N}$ such that $\big(S(F^l)\big)^{mt+1}({\bf n}(0), p(0))\in B\big(({\bf n}^{\ast}, K({\bf n}^{\ast})), \delta\big)$ for all $l\geq l_{\delta}$ and $t\geq t_{\delta}$.

(ⅱ) Assume ${\bf n}^{\ast}$ is unstable. Then, there exists $\delta_{0} > 0$ such that for every $\delta > 0$ satisfying $0 < \delta < \delta_{0}$ and any point $({\bf n}, p)$ such that $\left\Vert ({\bf n}, p)-({\bf n}^{\ast}, K({\bf n}^{\ast}))\right\Vert = \delta,$ there exists $l_{\delta}\in{\bf N}$ such that $\big(S(F^l)\big)^{m}({\bf n}, p)\notin \bar{B}\big(({\bf n}^{\ast}, K({\bf n}^{\ast})), \delta\big)$ for all $l\geq l_{\delta}$. In particular, $\bar{B}\big(({\bf n}^{\ast}, K({\bf n}^{\ast})), \delta\big)$ is not a stable set for $\big(S(F^l)\big)^{m}({\bf n}, p)$.

The theorem is stated in a general way for periodic solutions, although in most of the results it will be applied to equilibrium points. What this theorem tells us in the case of equilibrium points is the following. Always for $l$ large enough, if we have found a hyperbolic LAS equilibrium ${\bf n}^\ast$ of the reduced system (2.6) with a certain basin of attraction $B$, we can ensure that the solutions of the complete system (2.5) with initial conditions ${\bf n}(0)\in B$ and any $p(0) > 0$ will tend toward a neighborhood of point $\big({\bf n}^\ast, K({\bf n}^\ast)\big)$, which we can make as small as we want provided $l$ (timescales ratio) is big enough.

Although this is only approximately so, we will assume that if ${\bf n}^\ast$ helps us to define the asymptotic behavior of solutions with certain initial values ${\bf n}(0)$ in system (2.6), then $\big({\bf n}^\ast, K({\bf n}^\ast)\big)$ plays the same role in the complete system for all solutions that share these ${\bf n}(0)$, whatever $p(0)$ is.

A.1.   Proof of Proposition 3

Proof. Statement (ⅰ) is a direct consequence of Proposition 1 (ⅱ).

The instability of ${\bf 0}$ in statement (ⅱ) follows from Proposition 1 (ⅲ).

To prove the permanence of system (4.5), we apply Theorem 3 in ^[18]. Some of the hypotheses of the theorem are immediate to check:

- System (4.5) is continuous.

- $\bar{{\bf U}}_0$ (4.6) is irreducible, in fact, it is a positive matrix.

- ${\bf R}_+^2\setminus {{\bf 0}}$ is positively invariant for (4.5) since $\bar{{\bf U}}({\bf n})$, (4.4), is a positive matrix for all ${\bf n}$, and so if ${\bf n}\in {\bf R}_+^2\setminus {{\bf 0}}$, their product is positive.

- As $r_0 > 1$, the dominant eigenvalue of $\bar{{\bf U}}_0$ is greater than one.

It only remains to prove that (4.5) is dissipative to finish the proof. Following the definition 2 in ^[18], we have to prove that there exists a compact set $L\subset {\bf R}_+^2$ such that for all ${\bf n}(0)\in{\bf R}_+^2$, there exists a $\bar{t} = \bar{t}({\bf n}(0))$ satisfying ${\bf n}(t)\in L$ for all $t \geq \bar{t}$.

Starting from (4.5), we write the quotient $n_1(t+1)/n_2(t+1)$ in the following form:

where we have set $\sigma_2: = s_2/(1+a_2 k_2)$, and obtain the next inequalities:

and, substituting into the second equation of (4.5),

Let us now consider the scalar difference equation

and, for any initial conditions $\big(n_1(0), n_2(0)\big)\in {\bf R}_+^2$, let $\left\{\big(n_1(t), n_2(t)\big)\right\}$ be the corresponding solution of system (4.5) and $\left\{x(t)\right\}$ be the solution of Eq (A.7) that satisfies $x(0) = n_2(0)$. Without loss of generality, we assume that $n_2(0) > 0$ (if this is not the case, $n_2(1) > 0$ will be satisfied for $\big(n_1(0), n_2(0)\big)\in {\bf R}_+^2$). Since $g(x)$ is strictly increasing for $x\geq 0$, inequalities (A.5) and (A.6) imply, for $t = 1, 2, \ldots\ $,

The equation $g(x) = x$, apart from 0, has only one solution

If $x^*\leq 0$, as $\lim_{x\to\infty}(x-g(x)) = \infty$, then $g(x) < x$ for all $x > 0$. So, all solutions of (A.7) with $x(0) > 0$ would be decreasing and tend to 0. This would imply, using (A.8), that solutions of system (4.5) with $\big(n_1(0), n_2(0)\big)\in {\bf R}_+^2$ would tend to $(0, 0)$. Since this cannot be if $r_0 > 1$, we can ensure that $x^* > 0$. Now, using that $g(x)$ is increasing and $g'(x)$ decreasing, it is straightforward to see that for $0 < x < x^*$ we have $x < g(x) < g(x^*) = x^*$, and for $x > x^*$ we have $x^* < g(x) < x$. Then, Theorem 2.8 in ^[14] applies and, thus, all positive solutions of (A.7) are monotonic and tend to equilibrium $x^*$.

Finally, the previous calculations allow us to ensure that if we choose any $l > x^*$, the compact set $L$ whose existence we had to prove can be defined as

A.2.   Proof of Proposition 4

Proof. To prove the first assertion of the proposition, we need to apply Th. 1.2.7 and part (ⅰ) of Th. 1.2.8 in ^[12]. So, we have to check that the hypotheses (1.34) (a)–(e) are met, and that $\lim_{|{\bf n}|\to\infty, {\bf n}\in{\bf R}_+^2}\nu({\bf n}) = 0$. We start with the hypotheses (1.34) (a)–(e):

$({\rm{a}})$ $\bar{{\bf T}}({\bf n}) = \big(t_{ij}({\bf n})\big)$, $t_{ij}\in C^\infty (\Omega, (0, 1))$ ($\Omega$ is an open set of ${\bf R}^2$ containing ${\bf R}_+^2$, and $t_{1j}({\bf n})+t_{2j}({\bf n})\leq 1$ for $j = 1, 2$. It readily holds taking $\Omega = \{{\bf n}\in {\bf R}^2:1+a_1k_1+bn_2 > 0\}$.

$({\rm{b}})$ ${\bf I}-\bar{{\bf T}}({\bf n})$ is invertible for all ${\bf n}\in{\bf R}_+^2$. Immediate.

$({\rm{c}})$ $\bar{{\bf \Phi}}({\bf n}) = \big[\phi_{ij}({\bf n})\big]$, $\phi_{ij}\in C^\infty (\Omega, {\bf R}_+)$ holds taking $\Omega = \{{\bf n}\in {\bf R}^2:1+c_1n_1+c_2n_2 > 0\}$.

$({\rm{d}})$ $\bar{{\bf T}}({\bf n})+r_0\bar{{\bf \Phi}}({\bf n})$ is irreducible and primitive for all ${\bf n}\in{\bf R}_+^2$ and $r_0 > 0$. It is clear since all the entries of these matrices are positive.

$({\rm{e}})$ Straightforward computations show that, for all ${\bf n}\in{\bf R}_+^2$, the matrix $\bar{{\bf \Phi}}({\bf n}) \left({\bf I}-\bar{{\bf T}}({\bf n})\right)^{-1}$ has a positive, strictly dominant, simple eigenvalue

$\nu({\bf 0}) = 1$, with a nonnegative eigenvector $u({\bf n}) = \left(\begin{array}{c} 1 \\ 0 \end{array}\right)\geq {\bf 0}$, such that

Finally, from the expression for $\nu({\bf n})$, it is easy to deduce that

The second and third assertions of the proposition are a direct consequence, respectively, of parts 3 and 2 of the aforementioned Th. 1.2.8 in ^[12]. To verify that they can be applied, it is enough to prove that $\nu({\bf n}) < 1$ for all ${\bf n} > {\bf 0}$, which is immediate since

A.3.   Analysis of system (4.10)

Nine parameters appear in system (4.10), but we can group them and reduce their number to five. We leave $m$ and $f$ as they are, and define

The expression of system (4.10) with the new parameters is

The equilibria $E^* = (n_1^*, n_2^*)$ of system (A.12) are the solutions of the following system:

From the second equation, and equating $n_1/(q_1+ b_1 n_2)$ in both equations, we get

and, discarding equilibrium ${\bf 0}$, we obtain equilibrium $E^* = (n_1^*, n_2^*)$, where

$E^*$ is a feasible equilibrium whenever $(1-m-q_1)(1-\sigma_2)+mf > 0$, that is, if the following condition holds:

This condition is equivalent to the net reproductive number $n_0$, (4.7), being greater than one. From now on, we assume that $f > L$ holds.

We study the stability of $E^*$ through linearization. Direct calculations lead us to the following Jacobian matrix:

We let $\text{tr}\big(J(E^*)\big)$ and $\det\big(J(E^*)\big)$ denote the trace and determinant of matrix $J(E^*)$, respectively, and apply the Jury test ^[14]. We look for conditions so that the following double inequality is satisfied,

In the following, we will repeatedly use the conditions set on the parameters: $m\in (0, 1)$, $f > 0$, $q_1 > 1$, $b_1 > 0$, and $\sigma_2\in (0, 1)$.

The expression of $\det\big(J(E^*)\big)$ simplifies to

and, since $1-\det\big(J(E^*)\big) = \frac{(1-m)(1-\sigma_2)^2+(2-\sigma_2)mf}{(1-m)(1-\sigma_2)+mf} > 0,$ it follows that $1+\det\big(J(E^*)\big) < 2$.

The trace $\text{tr}\big(J(E^*)\big)$ can be expressed as

To check the first inequality in (A.17), we need to know the sign of $\text{tr}\big(J(E^*)\big)$. The only term that can be negative in expression (A.19) is $2\sigma_2-1$, so if $\sigma_2\geq 1/2$, then $\text{tr}\big(J(E^*)\big) > 0$. On the other hand, if $\sigma_2 < 1/2$, $\text{tr}\big(J(E^*)\big)\geq 0$ if, and only if,

In short, $\text{tr}\big(J(E^*)\big)\geq 0$ if either $\sigma_2\geq 1/2$, or $\sigma_2 < 1/2$ and $f\leq \bar{L}$, and, thus, $\text{tr}\big(J(E^*)\big) < 0$ if $\sigma_2 < 1/2$ and $f > \bar{L}$.

Let us first assume that $\text{tr}\big(J(E^*)\big)\geq 0$, and we then have

due to condition (A.15). So, if $\text{tr}\big(J(E^*)\big)\geq 0$, the feasible equilibrium $E^*$ is LAS.

Now assuming that $\text{tr}\big(J(E^*)\big) < 0$, i.e., $\sigma_2 < 1/2$ and $f > \bar{L}$, let us study the sign of

The denominator is positive, and the numerator is also positive whenever either $\sigma_2\geq 1/3$ or $\sigma_2 < 1/3$, and the following condition is met

In all these cases, $E^*$ is also LAS.

We can assure that $E^*$ is unstable if $\text{tr}\big(J(E^*)\big) < 0$, $\sigma_2 < 1/3$, and $f > M$. Since

we have that $\sigma_2 < 1/3$ implies that $L < \bar{L} < M$ and, therefore, the existence and instability conditions of the feasible equilibrium $E^*$ are simplified to $\sigma_2 < 1/3$ and $f > M$.

We bring together in the following proposition what we have seen so far about the positive equilibrium of system (A.12).

Proposition 9. System (A.12) possesses the feasible equilibrium $E^*$ (A.13) whenever $f > L$ (A.15).

(ⅰ) If either $\sigma_2\geq 1/3$, or $\sigma_2 < 1/3$ and $f < M$ (A.21), then $E^*$ is hyperbolic and LAS.

(ⅱ) If $\sigma_2 < 1/3$ and $f > M$, then $E^*$ is hyperbolic and unstable.

(ⅲ) If $\sigma_2 < 1/3$ and $f = M$, then $E^*$ is non-hyperbolic and the eigenvalues of $J(E^*)$ are

Proof. Claims (ⅰ) and (ⅱ) have been previously proven.

To prove claim (ⅲ), in (A.16), just replace $n_2^*$ with its value (A.14) and make $f = M$:

and a simple calculation gives us the result.

The last assertion reveals the possible existence of a period-doubling bifurcation when $f$ crosses the threshold value $M$, an issue that we justify below.

We assume that $f > L$ so that system (A.12) has the only feasible equilibrium $E^* = (n_1^*, n_2^*)$ (A.13). We apply the procedure developed in ^[23] (Section 5.4.2).

We first move $E^*$ to the origin. Let us apply the following change of variables, $X = (x_1, x_2) = (n_1-e_1, n_2-e_2)$, which converts system (A.12) into

with

This system can still be written in the following form:

where $A = Jg({\bf 0})$, the Jacobian matrix of map $g$ at equilibrium ${\bf 0}$. Thus, $G(X) = g(X)-Jg({\bf 0})X = O\big(\|X\|^2\big)$, and it can be expressed in the following form:

with $B(X, Y)$ and $C(X, Y, Z)$ multi-linear functions, defined, for $i = 1, 2$, as follows:

Following Proposition 9 (ⅲ), to check the existence of a period-doubling bifurcation, we assume from now on that $\sigma_2 < 1/3$ and $f = M$. In this case, matrix $A$ in (A.23) is the Jacobian matrix (A.22).

Let ${\bf v}$ and ${\bf w}$ be right and left eigenvectors, respectively, of the matrix $A$ associated with the eigenvalue -1, and normalized so that their scalar product is 1, $\langle {\bf v}, {\bf w}\rangle = 1$. For example:

In ^[23], from these vectors, a linear transformation of system (A.23) is carried out. With the help of the Taylor expansion (A.24), the center manifold of the system at equilibrium ${\bf 0}$ is approximated. After calculating the restriction of the system on this manifold, which is a scalar equation, it is transformed to its normal form. The structure of the normal form is $x(t+1) = -x(t)+c(x(t))^3+O((x(t))^4)$ and the coefficient $c$ determines the nondegeneracy of the period-doubling bifurcation and its direction. When $c > 0$, the bifurcation is supercritical, to the right, and, therefore, when the parameter $f$ becomes greater than $M$, the ${\bf 0}$ equilibrium is no longer stable and a stable 2-cycle appears in its place.

The calculation of $c$ can be performed using expression (5.49) in ^[23]

where $B$ and $C$ are defined in (A.25), ${\bf v}$ in (A.26), ${\bf w}$ in (A.27), $A$ is matrix (A.22), and $I$ the second order identity matrix.

Straightforward calculations give us

and

and, finally,

Taking into account the range of values of each parameter (A.11) and the assumption that $\sigma_2 < 1/3$, it is obvious that $c > 0$. Thus, we have proved the following proposition.

Proposition 10. Assume $\sigma_2 < 1/3$. Then, system (A.12) undergoes a period-doubling bifurcation at $f = M$.