On the global dynamics of a forest model with monotone positive feedback and memory
Franco Herrera Sergei Trofimchuk








a]Franco Herrera
a]Sergei Trofimchukmycorrespondingauthor
[mycorrespondingauthor] e-mails addresses:   franco.herrera@utalca.cl (F. Herrera); trofimch@inst-mat.utalca.cl (S. Trofimchuk, corresponding author)  

[a]Instituto de Matemáticas, Universidad de Talca, Casilla 747, Talca, Chile 






We continue to study (see <cit.>) a renewal equation  ϕ(t)= Fϕ_t proposed in <cit.> to model  trees growth. This time we are 
considering the case when the per capita reproduction rate β(x) is a non-monotone (unimodal) function  of tree's height x. Note that  
 the height of some species of trees can impact negatively seed viability, cf. <cit.>, in a kind of autogamy depression. 
Similarly to previous works, it is also assumed that the growth rate  g(x) of an individual of height x is a strictly decreasing function.   
 As in <cit.>, we analyse the connection between dynamics of the associated one-dimensional map F(b)=  Fb, b ∈_+, and the delayed (hence infinite-dimensional) model ϕ(t)= Fϕ_t.  
 Our key observation is that this model is of monotone positive feedback type since F is strictly increasing on _+ independently on the monotonicity properties of β. 
 
 size-structured model, semiflow, global attractor, stability, Allee effect 










§ INTRODUCTION
   

In this note, we consider the renewal equation 

    b(t) = Fb_t:= ∫_0^∞β( x_m + ∫_0^a  g( e^-μ(τ -a)∫_a^∞  e^-μ s b(t-s) ds ) dτ) e^-μ a b(t-a)da,

where  b_t(s):=b(t+s),  s ≤ 0, and functional  F is given by 
 
    Fϕ:=∫_0^∞β( x_m + ∫_0^a g( e^-μ(τ -a)∫_a^∞ e^-μ sϕ(-s) ds ) dτ) e^-μ aϕ(-a)da,
 This equation was proposed recently in <cit.> to model  trees growth and presents one of possible approaches to analyse the forest population dynamics, several other  models are mentioned in <cit.>.  In equation (<ref>), function  b(t) corresponds to the tree population birth rate at time t; continuous strictly decreasing function g(x), g: [0,+∞)→ (0, +∞), g(+∞)=0, is the growth rate  of an individual of height x;   μ >0 represents the per capita death rate;  globally Lipschitzian function β(x), β:[x_m,+∞)→ [0, +∞), is the per capita reproduction rate (depending only on height x).  This model postulates the minimal height x_m ≥ 0 for all newborn trees. By  <cit.>, b(t) gives the influx at x_m originating from reproduction by the extant population. To ensure good persistence properties of solutions, it suffices to assume that β(x) > 0 a.e. on [x_0,+∞)  for some x_0≥ x_m <cit.>. 


Clearly, (<ref>) can be considered as an autonomous functional equation with infinite delay  
once  it is  provided with non-negative initial conditions 
b(s)=ϕ(s) ≥ 0,  s ≤ 0, where ϕ belongs to the cone K of non-negative measurable functions of the Banach space 

    L^1_ρ(_-)= {ϕ: |ϕ|_ρ < ∞},      |ϕ|_q:= ∫_-∞^0|ϕ(s)|e^qsds),

with ρ < μ. It was proved in <cit.> that 
the renewal equation  (<ref>) generates a continuous semiflow 
 by  ( S^t ϕ)(s) = b(t+s), s ≤ 0. It has a zero 
steady state b(t) ≡ 0 while  the positive equilibria of S^t are defined from the equation 
 
    1= R(b):=∫_0^∞β( x_m + ∫_0^a g( e^-μτ/μb) dτ) e^-μ ada.

Note that  equation (<ref>)  has at least one  positive solution b_* if  
R(0) >1 and R(+∞)= β(x_m)/μ <1.  If additionally  β increases on _+, then R(b) is strictly decreasing so that solution b_* is unique, cf.  <cit.>.  Moreover, this unique positive steady state is locally stable and globally attracting solution of (<ref>), see  
<cit.>. The key argument in proofs here is that semiflow F^t ϕ is monotone  <cit.> once  β is increasing on _+. 


In some cases, however, the assumption of monotonicity of β is not realistic. For example,  the height of some species of trees can impact negatively seed viability, cf. <cit.>, in a kind of autogamy depression. Thus it would be interesting to study the dynamics in (<ref>) in the case when 
 the per capita reproduction rate β(x) is a non-monotone  (for instance, unimodal) function  of tree's height x. Now, a general description of this dynamics is known from  <cit.>, where it was proven that under all assumed (except monotonicity of β) restrictions on β, g and 
 0 ≤β(x_m)< μ,  R(0) >1,  the semiflow  F^t ϕ has a compact global attractor 𝒜 
attracting each non-zero solution. Moreover, each element ψ  of 𝒜 satisfies 
0 <θ_1 ≤ψ(s) ≤θ_2,    |ψ'(s)| ≤θ_3,  s ≤ 0,
for some  universal positive numbers θ_1 ≤θ_2, θ_3. 

In the sequel we present several analytical and numerical arguments shedding some light on possible dynamical structures of the global attractor 𝒜 without assumption of the monotonicity for β. As in <cit.>, we will analyse the connection between dynamics of the associated one-dimensional map F(b)=  Fb, b ∈_+, and the delayed (hence infinite-dimensional) model ϕ(t)= Fϕ_t.  






§ MONOTONICITY OF THE FUNCTION F(B)= BR(B),  B ∈_+. 
   
 
Assume all conditions  (except monotonicity of β) imposed on β and g in Section 1. Then function F(b)= bR(b) 
is strictly increasing on _+.  

 Assuming that ϕ(s) >0 a.e. on _- and setting  h(s)=ϕ(s)e^μ s and θ = H(u)=∫_-∞^uh(s)ds,  we find that 
    Fϕ=∫_-∞^0 β( x_m + ∫^0_u g( e^-μ sH(u))ds)h(u)du =
    ∫_0^H(0)β( x_m + ∫^0_H^-1(θ) g( e^-μ sθ)ds)dθ.

Clearly,  H(u) is strictly increasing and positive on ,  with H(-∞)=0, 

    d/dθH^-1(θ)= 1/h(H^-1(θ))   (0, H(0)).
  
If ϕ_2(s) > ϕ_1(s) >0, almost everywhere on _-, then 

    H_2(u)=∫_-∞^uh_2(s)ds > H_1(u)=∫_-∞^uh_1(s)ds,    u ≤ 0,

so that 

    H_2^-1(θ) < H_1^-1(θ) ≤ 0,       0 <θ≤ H_1(0).

Since g is strictly decreasing, we find that the function 

    γ(θ, ϕ):=x_m + ∫^0_H^-1(θ) g( e^-μ sθ)ds =x_m + ∫^μ^-1ln(1/θ)_H^-1(θ)+μ^-1ln(1/θ) g( e^-μ s)ds,

 is continuous and strictly decreasing in θ for each fixed ϕ >0, 

    ∂γ(θ, ϕ)/∂θ= -1/μ(g(θ)/θ - g(θ e^-μ H^-1(θ))[1/θ- μ/h(H^-1(θ))]  ) <0.

Also, for a fixed θ and ϕ_2(s) > ϕ_1(s) >0, it holds that 

    γ(θ, ϕ_2)> γ(θ, ϕ_1).

 Thus the equation u=γ(θ, ϕ) can be solved with respect to θ, we denote this solution by θ=σ(u,ϕ).  Clearly, σ is strictly decreasing in the first argument and  is strictly increasing in ϕ. Therefore g(σ(u,ϕ))/σ(u,ϕ) is strictly decreasing in ϕ
 
 In this way, 

    Fϕ = μ∫_x_m^+∞β(u) (g(θ)/θ - g(θ e^-μ H^-1(θ))[1/θ- μ/h(H^-1(θ))])^-1du,   θ=σ(u,ϕ).

Particularly, 

    F(b)=  Fϕ = μ∫_x_m^+∞β(u)θ/g(θ)du,   θ=σ(u,b).

so that F(b) is a strictly increasing function independently on the monotonicity properties of β. 

Theorem <ref> suggests to explore a possible connection between dynamical structures for the delayed equation (<ref>) and seemingly simpler equation 

    b(t)=  Gb_t:= ∫_0^∞ f(b(t-a-h)) e^-μ ada,

where F= μ^-1f:→ is a strictly increasing smooth function, f'(x)>0 for all x ∈,  having three equilibria  e_1 < e_2:=0 < e_3, F(e_j)= e_j.  Equation (<ref>) describes bounded solutions from the global attractor for the semiflow generated by the delayed equation 

    b'(t)= -μ b(t)+ f(b(t-h)),   μ >0.

Specific applications of  such a model can be found  in neural network theory <cit.>. Obviously, each equilibrium for  (<ref>) is determined from the equation 

    b=  Gb = μ ^-1f(b),

where the right hand side is given by a strictly increasing function. The global attractor 𝒜
 for (<ref>) has relatively simple structure described in <cit.> (if e_2 is "sufficiently unstable", then 𝒜 is 
a smooth solid 3-dimensional spindle separated by the invariant disk into the basins of attraction toward the tips e_1 and e_3). See also <cit.>
for further results  and references on this topic. 




§ CHARACTERISTIC EQUATIONS AND PRINCIPLE OF LINEARIZED STABILITY
  
Formal linearization of the operator F at an equilibrium b(t) ≡ b yields the following action 

    F'(b)ψ=∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ) e^-μ aψ(-a)da+

 
    ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)[g( b/μ)- g( e^-μ ab/μ)]∫_a^∞ e^-μ sψ(-s) ds da=

  
    β(x_m)∫_0^∞ e^-μ sψ(-s) ds+   g( b/μ) ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)∫_a^∞ e^-μ sψ(-s) ds da.

 Looking for the exponential eigenfunctions ψ(s) = exp(λ s) of the linearized equation ψ(t) =  F'(b)ψ_t, we obtain the characteristic equation for λ, cf.  <cit.>: 
 
    λ +(μ - β(x_m))=  g( b/μ) ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)e^-(μ+λ)ada =: χ_b(λ).

 

 
 
 First, we consider the simpler case when β is a non-decreasing function and  μ> β(x_m) (so that F'(b) ≤ 1 at the unique positive equilibrium): 
   Assume that β'(x_m+s) ≥ 0, g(s) >0, for s ≥ 0 and μ - β(x_m) >0. Then the characteristic equation has exactly one real solution λ=λ_0. Moreover, λ_0≤ 0 and λ_j < λ_0 for each other (complex) solution λ_j. If g is strictly decreasing on _+, then λ_0<0. 
  
  The proof is obvious if β'(s) = 0 a.e. so that we consider the case when β'≠ 0. Clearly each eigenvalues satisfies the inequality λ_j +(μ - β(x_m)) ≥ 0. Taking into account  opposite types of monotonicity of λ +(μ - β(x_m)) and χ_b(λ) and the fact that χ_b(λ) >0 and χ_b(0)= μ F'(b)-  β(x_m)≤μ -  β(x_m), we deduce the existence of a unique real (non-positive) solution λ_0. 
 
 Assuming that there is another eigenvalue λ_j with λ_j ≥λ_0, we get a contradiction: 
 
    λ_0 +(μ - β(x_m)) = |λ_0 +(μ - β(x_m))| < |λ_j +(μ - β(x_m))|≤

 
    g( b/μ) ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)e^-(μ+λ_j)ada ≤


    g( b/μ) ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)e^-(μ+λ_0a)da =λ_0 +(μ - β(x_m)).

Finally, suppose that g is a strictly decreasing and λ_0=0. 
Then we get the following contradiction: 

    μ - β(x_m)=   ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)g( b/μ)e^-μ ada <


    ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)g( e^-μ ab/μ)e^-μ ada =


    - β(x_m) + μ∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ)e^-μ ada = - β(x_m) + μ.

 
 Lemma <ref>  simplifies and improves arguments in   <cit.> (apparently,  the uniqueness statement of this lemma  contradicts to the existence of two negative eigenvalues in the particular case presented in <cit.>; 
however,  one of these eigenvalues is false and should be excluded to allow the convergence of integrals in χ_b(λ). 

 Next, we analyse  the situation when β is not necessarily  monotone function and  μ> β(x_m) (we assume that the model possesses at least one positive equilibrium b): 
  Assume that  g>0 is strictly decreasing Lipschitz continuous function on _+, μ > β(x_m). Then in the half-plane { z ≥ - μ} the characteristic equation has exactly one real solution λ=λ_0 and λ_j < λ_0 for each other (complex) solution λ_j. Moreover, if F'(b) < 1 then λ_0< 0, λ_0 = 0 when F'(b)  = 1 and  λ_0>0   if F'(b) > 1. 
  
 Note that 

    λ +(μ - β(x_m))=  g( b/μ)  ∫_0^∞β'( x_m + ∫_0^a g( e^-μτb/μ) dτ)e^-(μ+λ)ada =


    g( b/μ)  ∫_0^∞e^-(μ+λ)a/g( e^-μ ab/μ) dβ( x_m + ∫_0^a g( e^-μτb/μ) dτ) =


    -β(x_m) -  g( b/μ)  ∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ) [ e^-(μ+λ)a/g( e^-μ ab/μ)]'da.

Since 

    [ e^-(μ+λ)a/g( e^-μ ab/μ)]'= e^-(μ+λ)a-(μ+λ)g( e^-μ ab/μ) + e^-μ ab g'( e^-μ ab/μ)/(g( e^-μ ab/μ))^2 <0,

we obtain the following equivalent form of the characteristic equation: 

    1 =  g( b/μ)  ∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ)  e^-(μ+λ)a/g( e^-μ ab/μ)[1- e^-μ ab g'( e^-μ ab/μ)/(μ+λ)g( e^-μ ab/μ)]da.

The right side of this equation (we will denote it by ξ_b(λ), clearly 
ξ_b(λ) = (χ_b(λ)+β(x_m))/(μ+λ) and ξ_b(0)= (χ_b(0) + β(x_m))/μ = F'(b)) is a strictly decreasing function of λ on the interval (-μ, + ∞), ξ_b(-μ^+)= +∞, ξ_b(+∞)=0,  so that the equation  has exactly one real root λ_0 on this interval.
Furthermore, we see that λ_0 is a dominating eigenvalue. Indeed if λ_j is another eigenvalue with λ_j ≥λ_0 then 
λ_0+μ < |λ_j+μ| and we get a contradiction: 

    1 ≤   g( b/μ)  ∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ)  e^-(μ+λ_j)a/g( e^-μ ab/μ)[1+e^-μ ab |g'( e^-μ ab/μ)|/|μ+λ_j|g( e^-μ ab/μ)]da <


    g( b/μ)  ∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ)  e^-(μ+λ_0)a/g( e^-μ ab/μ)[1+ e^-μ ab |g'( e^-μ ab/μ)|/|μ+λ_0|g( e^-μ ab/μ)]da =1

Finally, since ξ_b(0)=F'(b) we conclude that λ_0 < 0 when F'(b)  < 1,  λ_0 = 0 when F'(b)  = 1 and λ_0 >0 when F'(b) > 1.  
Lemma <ref> implies that, except for the critical case F'(b)= 1, the local stability/instability of an equilibrium b in (<ref>) is equivalent to the local stability/instability of the same equilibrium for one-dimensional system b → F(b):   

  Assume that functions g: →, β: _+ → have a bounded and globally Lipschitzian first derivative and ρ < μ/5.  Then an equilibrium b of continuous semiflow 
  is  locally asymptotically stable if  F'(b) <1 and is unstable when F'(b) >1.  

 As it was established in <cit.>, under  additional assumptions of this corollary,  the functional F: L^1_ρ(_-)
 →  is continuously differentiable with bounded derivative.  Then an application of  <cit.> (based on Theorem 3.15 in <cit.>) completes the proof.    
 Note that the inequalities  F'(b) <1 and F'(b) >1 are equivalent, respectively, to the inequalities R'(b) <0 and R'(b) >0, where 

    R'(b)=
    -1/b -  1/μ bg( b/μ)  ∫_0^∞β( x_m + ∫_0^a g( e^-μτb/μ) dτ) [ e^-μ a/g( e^-μ ab/μ)]'da.




§ NUMERICAL SIMULATIONS AND OPEN PROBLEMS
  
Stability analysis and discussion presented in the previous sections suggest 
the existence of a 'large'  (we believe, open dense)  set 𝒪 in the phase space K
 such that solutions with initial data from 𝒪  converge to one of stable equilibria. Thus we can expect that a typical numerical solution of equation (<ref>) is asymptotically constant independently on the monotonicity properties of β. To check this conclusion and to find out possible configurations for the set of equilibria for (<ref>) with the unimodal fertility function β, we took β(x) = α xe^-x (which is used in the  Nicholson blowflies equation and other models, cf.  <cit.>);  x_m=0, μ=1, and g(x)=pe^-x, for some specific parameters α, p>0.  With these choices, our computations exhibit certain diversity in  dynamical patterns for (<ref>).


First, in full accordance with Theorem <ref>, Figure <ref> shows that  the function F is strictly increasing even in lack of monotonicity for β. Nevertheless, the function R is not necessarily monotone and we can observe that for different values of p,  R can be a decreasing function as well as  a unimodal one, see Figure <ref>. 
 It is important to mention that these shapes for R  do not change essentially even if instead of g(x)=pe^-x we take g with other orders of decay (e.g. polynomial) at +∞. 

Interestingly,  regardless the fact R(0) <1, equation (<ref>)  with α=6 and p=5 has three equilibria, b_1=0 < b_2≈0.47<b_3≈ 3.2. Since F'(0)=R(0)<1, F'(b_2) = 1+ b_2R'(b_2) >1, F'(b_3) = 1+ b_3R'(b_3) <1, we deduce that equilibria b_1=0 and b_3 are locally asymptotically stable and b_2 is unstable (see Corollary <ref>).  Our computations show that the linearization of  (<ref>)  at b_2 has only one (real) eigenvalue with non-negative real part for  α=6, p=5, therefore we cannot invoke the Hopf bifurcation approach to prove the existence of a periodic orbit for (<ref>).  

Simulating solutions b(t) of (<ref>) for different initial data ϕ, we took constant and periodic initial functions  close to the intermediate equilibrium b_2. Figures 3 and 4 show that these solutions are asymptotically constant. Clearly, the nonlinearity F(b) is bistable here and we observe a kind of the Allee effect for the forest dynamics. 
























§ ACKNOWLEDGMENTS
    This research was supported in part by the projects FONDECYT 1231169 and AMSUD220002 (ANID, Chile). The first author was also supported by  ANID-Subdirección de Capital Humano/Doctorado Nacional/2024–21240616. 
 The second author gratefully acknowledges the hospitality and the support of the MATRIX Institute allowing him to participate in the Research Program "Delay Differential Equations and Their Applications" (Australia, Creswick,  12 –20 Dec 2023). 






99
FBCD C. Barril, À. Calsina, O. Diekmann, J. Z. Farkas,  On competition through growth reduction, e-print 	arXiv:2303.02981, 
https://doi.org/10.48550/arXiv.2303.02981

FBCD1
C. Barril, A. Calsina, O. Diekmann,  J. Z. Farkas, On
  hierarchical competition through reduction of individual growth, Journal of
  Mathematical Biology, (2024) 88:66 https://doi.org/10.1007/s00285-024-02084-x

CO M. A. Caraballo-Ortiz, E. Santiago-Valentín, T. A. Carlo, 
Flower number and distance to neighbours affect the fecundity of Goetzea elegans (Solanaceae)
Journal of Tropical Ecology, 27 (2011) 521–528. 



DG  O. Diekmann, M. Gyllenberg, Equations with infinite delay: blending the abstract and the concrete, J.  Diff. Equations, 252 (2012), 819–851.




 


HT  F. Herrera, S. Trofimchuk, Dynamics of one-dimensional maps and Gurtin-MacCamy's population model. Part I: asymptotically constant solutions, Ukrainian Math. J.  75 (2023), 1635–1651 https://doi.org/10.3842/umzh.v75i12.7678

HTa  F. Herrera, S. Trofimchuk, Global dynamics of a size-structured forest model, e-print 	arXiv:2401.08618, 
 https://doi.org/10.48550/arXiv.2401.08618

KWW T. Krisztin, H.-O. Walther and J. Wu,  Shape, smoothness and invariant stratification of an attracting set for delayed monotone positive feedback, Fields Institute Monograph Series, Vol. 11, AMS, Providence, RI, 1999. 

KW T. Krisztin, H.-O. Walther,  Unique periodic orbits for delayed positive feedback and the global attractor.
J. Dyn. Differ. Equ. 13, 1–57 (2001).

KV T. Krisztin,  G. Vas,  The unstable set of a periodic orbit for delayed positive feedback, J. Dyn. Differ. Equ.  28,  (2016) 805–855.

MM  Z. Ma, P. Magal, Global asymptotic stability for Gurtin-MacCamy's population dynamics model, Proceedings of the AMS,152 (2024), 765–780.



MD J.A.J Metz,  O. Diekmann, The dynamics of physiologically structured populations, Lecture Notes in Biomathematics, 68, 1986.


Smith  H.L. Smith,   Monotone dynamical systems, Amer. Math. Soc.,  Providence, 1995.









