Abstract

We consider a class of evolution equations describing population dynamics
in the presence of a carrying capacity depending on the population with
delay. In an earlier work, we presented an exhaustive classification of the
logistic equation where the carrying capacity is linearly dependent on the
population with a time delay, which we refer to as the
“linear delayed carrying capacity” model. Here, we generalize it to the
case of a nonlinear delayed carrying capacity. The nonlinear functional form
of the carrying capacity characterizes the delayed feedback of the evolving
population on the capacity of their surrounding by either creating additional
means for survival or destroying the available resources. The previously
studied linear approximation for the capacity assumed weak feedback, while
the nonlinear form is applicable to arbitrarily strong feedback. The
nonlinearity essentially changes the behavior of solutions to the evolution
equation, as compared to the linear case. All admissible dynamical regimes are
analyzed, which can be of the following types: punctuated unbounded growth,
punctuated increase or punctuated degradation to a stationary state, convergence
to a stationary state with sharp reversals of plateaus, oscillatory attenuation,
everlasting fluctuations, everlasting up-down plateau reversals, and divergence
in finite time. The theorem is proved that, for the case characterizing the
evolution under gain and competition, solutions are always bounded, if the
feedback is destructive. We find that even a small noise level profoundly
affects the position of the finite-time singularities. Finally, we
demonstrate the feasibility of predicting the critical time of solutions
having finite-time singularities from the knowledge of a simple quadratic
approximation of the early time dynamics.

There exists a number of evolution equations characterizing population dynamics,
which have been applied to numerous concrete systems, ranging from populations
of human and biological species to the development of firms and banks (see review
articles [Kapitza(1996), Hern(1999), Korotayev(2007), Yukalov et al.(2012a)]). The mathematical
structure of these equations usually represents some generalizations of the
logistic equation. Such equations can be classified into three main classes,
depending on the nature of the dynamics of the carrying capacity. The first
class, independently of whether the growth rate is a nonlinear function without
or with delay of the population, assumes that the carrying capacity is a constant
quantity given once for all that describes the total resources available to the
population, in agreement with the initial understanding of the carrying capacity
(e.g., [Haberl & Aubauer(1992), Varfolomeyev & Gurevich(2001), Hui & Chen(2005), Gabriel et al.(2005), Berezowski & Fudala(2006), Arino et al.(2006)]).
The second class allows the carrying capacity to change as a function of time, but
for exogenous reasons, either by explicitly prescribing its evolution or by specifying
its own independent dynamics for instance also given by a logistic equation
[Dolgonosov & Naidenov(2006), Pongvuthithum & Likasiri(2010)].

The third class of equations interprets the carrying capacity as a functional
of the population itself, implying that the population does influence the carrying
capacity, either by producing additional means for survival or by destroying the
available resources [Yukalov et al.(2009), Yukalov et al.(2012a), Yukalov et al.(2012b)]. This feedback makes
it possible to describe the regime of punctuated evolution, which is often observed
in a variety of biological, social, economic, and financial systems
[Yukalov et al.(2009), Yukalov et al.(2012a)].

In our previous articles [Yukalov et al.(2009), Yukalov et al.(2012a)], the carrying capacity was
approximated by a linear dependence on the population variable, which, strictly
speaking, assumes that the population influence should be smaller than the
initial given capacity. In the present paper, we generalize the approach by
accepting a nonlinear carrying capacity that allows us to consider a population
influence of arbitrary strength. Moreover, the population variable enters the
carrying capacity with a time delay, since to either create or destroy resources
requires time.

The outline of the present paper is as follows. In Sec. 2, we explain the
deficiency that is typical of the linear capacity specification and suggest the
way of generalizing its form to a nonlinear expression, by the use of the theory
of self-similar approximations. The existence and general conditions for the
stability of evolutionary stationary states are formulated in Sec. 3. The
temporal behavior of solutions essentially depends on the system parameters
characterizing different prevailing situations, when the main features are
described as gain and competition (Sec. 4), loss and cooperation (Sec. 5), loss
and competition or gain and cooperation (Sec. 6). The equations display a rich
variety of dynamical regimes, including punctuated unbounded growth, punctuated
increase or punctuated degradation to a stationary state, convergence to a
stationary state with sharp reversals of plateaus, oscillatory convergence,
everlasting fluctuations, everlasting up-down plateau reversals, and divergence
in finite time. All admissible dynamical regimes are studied and illustrated. The
role of noise on the dynamics is investigated in Sec. 7. The possibility of
predicting finite-time singularities by observing only the initial stage of
motion is discussed in Sec. 8. Section 9 concludes.

Here and in what follows, we use the dimensionless variable for the population
x(t) as a function of time t. The reduction of the dimensional equation to
the dimensionless form has been explained in full details in our previous papers
[Yukalov et al.(2009), Yukalov et al.(2012a)] and we do not repeat it here.

2.1 Singular solutions for linear delayed carrying capacity

The general expression for the evolution equation with functional carrying
capacity, in dimensionless notation, reads as

dx(t)dt=σ1x(t)−σ2x2(t)y(x),

(1)

where the carrying capacity functional

y(x)=y[x(t−τ)]

(2)

depends on the population at an earlier time, with a constant delay time τ,
which embodies that any influence of the population on the capacity requires
time in order to either create additional means or to destroy the given resources.
Here and in what follows, we use the term population, although the variable
x can characterize either population, or firm assets, or other financial and
economic indices [Yukalov et al.(2009), Yukalov et al.(2012a)].

The coefficients σi describe the prevailing features in the balance between
gain (birth) or loss (death) and competition versus cooperation. There exist four
situations characterized by these coefficients:

If we fix σ1=σ2=1 and take a constant capacity y0=y(0),
we come back to the logistic equation. In our previous papers
[Yukalov et al.(2009), Yukalov et al.(2012a)], we modelled the influence of population on the
carrying capacity by the linear approximation

y1(x)=1+b1x(t−τ),

(4)

with the parameter b1 describing either destructing action of population on
the resources, when b1<0, or creative population activity, if b1>0.

In the case of destructive action, it may happen that the capacity (4) reaches
zero and becomes even negative, with the effect of the growing population. It
crosses zero at a critical time tc defined by the equation

1+b1x(tc−τ)=0.

Being in the denominator of the second term of (1), the vanishing capacity leads
to the appearance of divergent or non-smooth solutions. In some cases, having to
do with financial and economic applications, the arising negative capacity can
be associated with the leverage effect [Yukalov et al.(2012a)]. However, in the usual
situation, the solution divergencies, caused by the zero denominator, look
rather unrealistic, reminding of mathematical artifacts. Therefore, it would
be desirable to define a carrying capacity that would not cross zero in finite
time.

2.2 Evolution equation with nonlinear delayed carrying capacity

It would be possible to replace the linear form (4) by some nonlinear
function. This, however, would be a too arbitrary and ambiguous procedure, since
it would be always unclear why this or that particular function has been chosen.
In order to justify the choice of a nonlinear function, we propose
the following procedure to select the form of the nonlinear carrying capacity.

Strictly speaking, the linear approximation presupposes that the second term
in the r.h.s of expression (4) is smaller than one. Generally, a function
y(x) can be expanded in powers of its variable x according to the series

y(x)≃1+b1x+b2x2+…,

(5)

where different terms describe the influence with different action intensity.
Since expression (4) can be interpreted as the first-order term in the
general series expansion (5) with, in principle, infinite many terms,
it is convenient to think of it as the general expansion of some nonlinear
function to be determined by a suitable summation. We thus propose to construct
a nonlinear extension of (4) by defining an effective sum of these
series (5). A standard way to realize the summation (5) is via
Padé approximants [Baker & Graves-Moris(1996)]. However, as is well known, the Padé
approximants are very often plagued by the occurrence of artificial zeroes and
divergencies, which makes them inappropriate for the summation of a quantity that
is required to be finite and positive. For this purpose, it is more appropriate
to resort to the method of self-similar approximations
[Yukalov(1990a), Yukalov(1990b), Yukalov(1991), Yukalov(1992), Yukalov & Yukalova(1996)]. This is a
mathematical method allowing for the construction of effective sums of power-law
series, even including divergent series. According to this method, a series
(5) is treated as a trajectory of a dynamical system, whose fixed point
represents the sought effective sum of the series. In the vicinity of a fixed
point, the trajectory becomes self-similar, which gives its name to the method
of self-similar approximations.

Employing this method, under the restriction of getting a positive effective
sum of the series, which is done by using the self-similar exponential
approximants [Yukalov & Gluzman(1997), Yukalov & Gluzman(1998)], we obtain the effective sum

y(x)=exp{bx(t−τ)}.

(6)

Here the production parameterb characterizes the type of the influence
of the population on the carrying capacity. In the case of creative activity of
population, producing additional means for survival, the creation parameter is
positive, b>0. And if the population destroys the given carrying capacity,
then the destruction parameter is negative (b<0).

In this way, we come to the evolution equation

dx(t)dt=σ1x(t)−σ2x2(t)exp{−bx(t−τ)},

(7)

with the nonlinear carrying capacity that generalizes the delayed logistic
equation (1) with the linear capacity. The proposed generalization (7)
allows us to consider population influences of arbitrarily strong intensity.
The initial condition to the equation defines the history

x(t)=x0(t≤0).

(8)

By definition, the population is described by a positive variable, so that we
will be looking for only positive solutions x(t)>0 for t>0.

This section gives the general conditions for the existence and stability of
stationary states of the delayed Eq. (7). The details of such conditions depend
on the type of the system characterized by the values of σi. More
specific investigations of the evolutionary stable states as well as the overall
dynamical regimes will be analyzed in the following sections.

3.1 Existence of stationary states

The stationary states of Eq. (7) are defined by the solutions to the equation

σ1x∗−σ2(x∗)2exp(−bx∗)=0.

(9)

There always exists the trivial solution

x∗1=0(−∞<b<∞,σ1=±1,σ2=±1).

(10)

But the nontrivial solutions require the validity of the relation

σ1σ2=x∗exp(−bx∗),

(11)

which tells us that they may happen only for

σ1=σ2=±1(x∗>0).

(12)

Two nontrivial states may exist depending on the value of the production
parameter. One of the states exists for negative values of the latter in the
range

0<x∗2≤1(b≤0).

(13)

Below we show that varying the parameters of Eq. (7) generates a number of
bifurcations and provides a rich variety of qualitatively different solutions.
We give a complete classification of all possible solutions demonstrating how
the bifurcation control [Chen et al.(2003)] can be realized for this equation.

For positive b, there can occur two states, such that

1<x∗2<e,x∗3>e(0<b<1e).

(14)

At the bifurcation point b=1/e, the states coincide: x∗2=x∗3=e.
There are no stationary nontrivial states for b>1/e. The existence of the
nontrivial states is illustrated in Fig. 1.

Figure 1: Existence of nontrivial stationary states depending on the value
of the production parameter b.

3.2 Stability of stationary states

Studying the stability of the stationary states for the delay differential
equation, we follow the Lyapunov-type procedure developed for delay equations
in the book [Kolmanovskii & Myshkis(1999)] and described in detail in Refs.
[Yukalov et al.(2009), Yukalov et al.(2012a)]. The steps of this procedure are as follows. We
consider a small deviation from the solution of equation (7) by writing
x=x∗+δx. Substituting this in (7) yields the linearized equation

ddtδx(t)=(σ1−2σ2x∗e−bx∗)δx(t)+bσ1x∗δx(t−τ),

(15)

which is the basis for the stability analysis. The solution of equation (7)
is Lyapunov stable, when the solution to (15) is bounded. A fixed point is
asymptotically stable, when the solution to (15) converges to zero, as time
tends to infinity. Thus, the conditions of fixed-point stability are prescribed
by the convergence to zero of the small deviation described by equation (15).

We find that the trivial state x∗1=0 is stable when

σ1=−1,σ2=±1,−∞<b<∞,τ≥0,

(16)

while for σ1=1 it is always unstable.

In the case of the nontrivial states, we use relation (11) and reduce Eq. (15)
to

ddtδx(t)=−σ1δx(t)+bσ1x∗δx(t−τ).

(17)

According to the existence condition (12), we need to study the stability
of the nontrivial states only for coinciding σi.

The following analysis of the stationary states and the solution of the full
evolution equation (7) requires to specify the values of σi.

4.1 Evolutionary stable states

The stability analysis shows that the nontrivial state x∗2 is stable either
in the range

0<x∗2≤1e(b≤−e,τ<τ∗2)

(18)

or in the range

1e<x∗2≤e(−e<b≤1e,τ>0),

(19)

where

τ∗2=1√(bx∗2)2−1arccos(1bx∗2).

(20)

Close to the boundary b→−e one has

x∗2≃1e(1+b+e2e),τ∗2≃π√e|b|−e(b→−e).

The state x∗3>e is always unstable under σ1=σ2=1.

Thus, there can exist just one stable stationary state x∗2∈(0,e),
whose region of stability is given by Eqs. (18) to (20) and shown
in Fig. 2.

Figure 2: Stability region (shadowed) for the stationary state x∗2 under
σ1=σ2=1.

For some values of the production parameter b, the basin of attraction of
x∗2 is not the whole positive semiline of x0, but a limited interval.
This happens for b∈(0,1/e), when the basin of attraction is defined by
the inequalities

0<x0<x∗3(0<b≤1e).

(21)

For these b values, the solution x tends to infinity, as t→∞,
for the history x0>x∗3, even for the parameters b and τ in the
region of stability of x∗2.

4.2 Boundedness of solutions for semi-negative production parameters

When σ1=σ2=1 and when the production parameter b is non-positive,
this means that the population does not produce its carrying capacity but rather
destroys it or, in the best case, retains the given capacity value. In this case,
an important result for the overall temporal behavior of solutions can be derived
rigorously: the population growth has to be limited.

{proposition}

The solution x(t) to the evolution equation (7), under the condition
σ1=σ2=1, for b≤0, any finite τ≥0, and any history
x0≥0, is bounded for all times t≥0, and, for b<0, there exists
a time
t0=t0(x0,τ) such that

0≤x(t)≤1(t≥t0).

(22)

{proof}

When b=0, the explicit solution is

x(t)=x0x0+(1−x0)e−t(b=0).

If x0<1, then x→1 from below as t→∞. If x0=1,
then x=1 for all t>0. If x0>1, then x→1 from above, as
t→∞. So, the solution is always bounded.

When b is arbitrary and t≤τ, then the explicit solution

x(t)=x0ebx0+tebx0+x0(et−1)

is evidently bounded.

For b<0, any finite τ≥0 and all t≥0, the evolution
equation (7) reads as

dxdt=x−x2exp{|b|x(t−τ)}.

(23)

If, at some moment of time t>0, it happens that x≥1, then the above
equation defines a semi-negative derivative

dxdt≤x(1−exp{|b|x(t−τ)})≤0,

implying that x decreases or does not grow.

If the history is such that x0<1, then either x stays always below one,
or it grows and reaches one at some moment of time t>0. But x cannot cross
the line x=1, since, as is shown above, at the time when x would become
≥1, it has to either stay on this line x=1 or has to decrease. The
solution cannot stay forever on the unity line, as far as x=1 is not a stable
stationary state, which is x∗2<1 for b<0. This means that there is a
moment of time t0<+∞ such that the solution x has to go down for
t>t0.

For x0=1, again the solution cannot rise, having a non-positive derivative,
and cannot stay forever at this value, which is not a stable fixed point.
The sole possibility is that x starts diminishing beyond a finite t0.

When x0>1, then its derivative is non-positive. The solution cannot grow and
cannot stay forever at the value that is larger than the fixed point x∗2<1.
Hence x must decrease. When diminishing, it reaches the value x=1, where
it cannot stay forever, but has to go down beyond a finite time t0.

Thus, for b<0, there always exists such a moment of time, when the solution
goes below the value x=1 and can never cross this line from below.

4.3 Punctuated unbounded growth

A different situation happens for positive production parameters, when the
solutions may be unbounded. For example, when b is outside of the stability
region of x∗2, so that

b>1e,τ≥0,x0>0,

(24)

then x tends by steps to infinity, as t→∞. A similar unbounded
punctuated growth occurs when b is inside the stability region, but x0
is outside of the attraction basin of x∗2, which takes place if

The meaning of such a behavior is clear: for either sufficiently high production
parameter or for sufficiently high startup and creative activity, the population
(or firm development) can demonstrate unlimited punctuated growth with time.

4.4 Punctuated convergence to stationary state

For positive production parameters, the solutions can also tend to the
stationary state x∗2 by punctuated steps. They tend to x∗2 from below
if x0<x∗2, and from above if x0>x∗2. This happens when

0<b<1e,τ≥0,x0<x∗3.

(26)

When the production parameter is negative, the approach to the stationary
state becomes quite nonmonotonic, with sharp reversals after almost horizontal
plateaus. This regime arises when

−e<b<0,τ≥0,x0>0.

(27)

With decreasing further the negative destruction parameter, the plateaus shorten
and the dynamics reduce to strongly fluctuating convergence to the focus x∗2,
which occurs for

Different regimes of punctuated convergence to the stationary state x∗2
are shown in Fig. 4, which include punctuated growth, punctuated decay, and
convergence with plateau reversals. Strongly oscillatory convergence to x∗2
is demonstrated in Fig. 5.

Figure 5: Strongly oscillatory convergence to the stationary state x∗2=0.383,
under σ1=σ2=1, with b=−2.5<−e and x0=0.25, for different
time lags: (a) τ=20; (b) τ=120.

4.5 Everlasting oscillations

When the destructive action is rather strong, increasing the time delay leads
to the switch from the oscillatory convergence to a stationary state to the
regime of everlasting oscillations, which happens for

b<−e,τ≥τ∗2,x0>0.

(29)

Figure 6 illustrates this effect of changing the dynamical regime with increasing
the time lag above τ∗2 that plays the role of a critical point.

A rather exotic regime develops when the feedback action on the carrying
capacity is strongly destructive and the time lag is very long. Then, there
appears the regime of everlasting up-down reversals of plateaus located at zero
and one, as is shown in Fig. 7.

Figure 7: Everlasting up-down plateau reversals, σ1=σ2=1,
for the destruction parameter b=−10≪−e, history x0=0.1, and
the time lag τ=100≫τ∗2=1.524.

5.2 Evolutionary stable states

Under loss and cooperation, the state x∗2 is always unstable. But
the trivial stationary state x∗1=0 is also stable for all b and τ.
The nontrivial state x∗3>e is stable for

0<b<1e,τ<τ∗3,

(31)

where

τ∗3=1√(bx∗3)2−1arccos(1bx∗3).

(32)

On the boundary of stability, when b approaches 1/e, then

x∗3≃e[1+√2(1−be)],τ∗3≃1−23√2(1−be)(b→1e).

The stability region of x∗3 is presented in Fig. 9. Because the state
x∗1 is stable everywhere, hence the stability region of x∗3 corresponds
to the region of bistability.

Figure 9: Stability region of the stationary state x∗3 (shadowed),
under σ1=σ2=−1. Since x1=0 is always stable,
the shadowed region is also the region of bistability.

5.3 Dynamical regimes of evolution

In the case of loss and cooperation, depending on the parameters b, τ,
and the history x0, there can occur the following dynamical regimes: monotonic
convergence to a stationary state, convergence with oscillations, everlasting
oscillations, and finite-time singularity. All these regimes are demonstrated
in Figs. 10 to 13. The summary of all possible solution types is illustrated
in the scheme of Fig. 14.

Figure 11: Change of dynamical regimes, under σ1=σ2=−1,
for the history x0>x∗3=5.938, with varying time lags:
(a) convergence to the stationary state x∗1=0 for b=0.3,
τ=20, x0=5.94 (solid line), x0=7 (dashed line), and
x0=10 (dashed-dotted line); (b) finite-time singularity at
tc=18.09 for b=0.3, x0=7, τ=0.8 in the interval
τ1<τ<τ2, with τ1=0.777, τ2=1.398;
(c) everlasting oscillations for b=0.3, x0=7, τ=0.68>τ∗3=0.661;
(d) transformation of the oscillatory convergence to the state x∗3=5.938
for b=0.3, x0=7 and τ=0.6>τ3=0.31 (solid line) to the monotonic
decay to the same state for τ=0.1<τ3=0.31 (dashed line).

Figure 12: Dynamical regimes, under σ1=σ2=−1, for the
history x∗2<x0=2<x∗3, with x∗2=1.631, x∗3=5.938,
the fixed production parameter b=0.3, varying only the time lag:
(a) finite-time singularity at tc=2.42 for τ=10>τ2=0.713277
(dashed line) and τ=0.713278>τ2, with the singularity at
tc=3.967 (solid line);
(b) everlasting oscillations for τ∗3<τ=0.68<τ2,
with τ∗3=0.661, τ2=0.713277; (c) convergence to the
stationary state x∗3 for τ1<τ=0.65<τ∗3, with
τ1=0.285 (solid line) and τ=0.1<τ1 (dashed line).

Figure 13: Dependence of dynamical regimes, under σ1=σ2=−1,
in the case of the destruction parameter b=−0.25, on the history and time
lags: (a) finite-time singularity at tc=1.509 for x0=1>x∗2=0.816,
τ=10>tc (solid line) and the singularity at tc=1.223 for
τ=0.1<tc ( dashed line); (b) monotonic degradation to x∗1=0 for the
same b=−0.25, but for x0=0.81<x∗2, τ=5 (solid line), for
x0=0.5<x∗2, τ=0.1 (dashed line), and x0=0.81<x∗2,
τ=0.1 (dashed-dotted line). {sidewaysfigure}

Classification of all possible dynamical regimes, under
σ1=σ2=−1, depending on the production parameter b,
time lag τ, and history x0.

6 Dynamics under loss and competition or gain and cooperation
(σ1σ2=−1)

Contrary to the previous cases, where σ1 and σ2 were equal, now they
are of opposite signs, so that σ1σ2=−1.

6.1 Decay under loss and competition

In such an unfortunate situation, when

σ1=−1,σ2=1,

(33)

there is only the stationary state x∗1=0 that is stable for all production
parameters b∈(−∞,∞), any time lag τ≥0, and arbitrary
history x0≥0. The solutions always monotonically decay to zero, as shown
in Fig. 15.

Figure 14: Monotonic decay to zero, under loss and competition
(σ1=−1 and σ2=1), with the history x0=5 and time
lag τ=10, for b=1 (solid line) and b=−0.1 (dashed line).

When σ1=1,σ2=−1, the solutions always grow with time, exhibiting
either unbound increase as t→∞, or a finite-time singularity at a
critical time t∗c.

If b<0, the solutions exhibit finite-time singularities. For τ≥tc,
where

tc=ln(1+ebx0x0),

(34)

the point of singularity is t∗c=tc(x0,b). But if τ<tc, then
the singularity point, defined numerically, is t∗c(x0,b,τ)≥tc,
such that t∗c→tc+0, as τ→tc−0. The corresponding
finite-time singularities are illustrated in Fig. 16a.

When b>0, there can occur either unbounded growth as t→∞ or a
finite-time singularity. If 0≤τ<τ∗, where τ∗=τ∗(x0,b)
is defined numerically, the solution unboundedly grows as t→∞. For
τ∗≤τ<tc, there exists a critical time t∗c=t∗c(τ)>tc
at which the solution diverges. And when τ≥tc, the divergence happens
at t∗c=tc, where tc is given by expression (38). The change of
the regime for the same production parameter b and history x0, but a varying
time lag τ is illustrated in Fig. 16b.

Figure 15: Fig. 16. Behavior of the solutions in logarithmic scale, under σ1=1,
σ2=−1, exhibiting either finite-time singularities or unbounded
growth at t→∞: (a) finite-time singularity, with the divergence
at tc=0.343 for b=−0.1, x0=2, τ=10>tc (solid line)
and with the divergence at t∗c=0.256 for b=−0.1, x0=2,
τ=0.01<tc (dashed line); (b) the change of regime for b=5,
x0=2 and varying time lag, when the finite-time singularity at tc=9.307
for τ=20>tc (dashed line), or at t∗c=9.655>tc for
τ∗<τ=9.212<tc, where τ∗=9.2117 (dashed-dotted line),
changes to the unbounded growth as t→∞ for τ=9.211<τ∗
(solid line).

7.1 Stochastic differential equation

In the presence of noise, the evolution equation (7) becomes the stochastic
differential equation

dx=g(x,t)dt+αdWt,

(35)

where α=√2D characterizes the noise strength, with D being
the diffusion coefficient, and where

g(x,t)=σ1x−σ2x2exp{−bx(t−τ)}.

(36)

We consider Eq. (35) in the sense of Ito, where Wt is the standard
Wiener process.

The addition of the noise does not influence much those solutions that do not
exhibit finite-time singularities, while the latter can be strongly influenced
by even weak noise. Therefore we concentrate our attention on the most
interesting case of the noise influence on the solutions with finite-time
singularities.

7.2 Influence of noise on finite-time singularities

Recall that a finite-time singularity can occur only in the case of cooperation,
when σ2=−1. Even rather weak noise can essentially shift the
singularity point. Moreover the same noise strength, in different stochastic
realizations, shifts the singularity point in a random way. Under the occurrence
of a finite-time singularity, the influence of noise turns out to be more
important than the variation of the time lag. This is in agreement with the Mao
theorem [Mao(1996)], according to which there can exist a finite range of time
lags for which the solution to the differential delay equation is close to that
of the related ordinary differential equation. Figure 17 illustrates the
influence of noise on the singularity point.

Figure 16: Influence of noise, in the case of b=−0.25 and x0=1, on the
point of the finite-time singularity: (a) several realizations of stochastic
trajectories, under σ1=σ2=−1 and τ>tc=1.509 for
the same noise strength α=0.25; (b) stochastic trajectories, under
σ1=1, σ2=−1, with τ≪tc=0.371, for the same
noise strength α=0.5; (c) stochastic trajectories, with the parameters
as in (b), but for the larger noise strength α=1; (d) singular solutions
for the parameters as in (b), but with α=0, for different time lags,
τ=0 (solid line), τ=0.01 (dashed line), τ=0.1 (dotted line),
and τ=1 (dashed-dotted line).

7.3 Fokker-Planck equation and condition for existence of a stationary
probability distribution

The stochastic differential equation (35) corresponds to the Fokker-Planck
equation

∂∂tP(x,t)=−∂∂x[g(x,t)P(x,t)]+D∂2∂x2P(x,t)

(37)

for the distribution function P(x,t), with D=α2/2. Looking for a
solution that would be defined for all x∈[0,∞) and all t≥0
requires the existence of a stationary distribution

P(x)≡limt→∞P(x,t).

(38)

The latter is defined by the equation

D∂2∂x2P(x)−∂∂x[g(x)P(x)]=0,

(39)

in which

g(x)=σ1x−σ2x2exp(−bx).

(40)

Equation (39) has to be complemented by boundary conditions, which are
usually taken in the form of absorbing boundary conditions

limx→∞P(x)=0,limx→∞g(x)P(x)=0,

(41)

since these conditions allow for the normalization of the distribution as

If b=0, then Eqs. (39) and (40) show that the limit (45) requires the same
conditions (48) as for b<0.

The stationary distribution (43) possesses maxima when the effective
potential (44) displays minima. The latter correspond to the stable fixed
points x∗ of the differential equation for x(t). Combining the analysis of
the conditions for the existence of the distribution P(x) with the conditions
for the existence of stable fixed points of the delay differential equation (7),
we find the following: When the solution x(t), for any history x0 and
τ→0, converges to a fixed point, then P(x) exists. Conversely, when
there is, at least for some history x0 and τ→0, an unbound solution
x(t), diverging either at a finite-time singularity or at increasing time
t→∞, then P(x) does not exist. Summarizing, we come to the
conclusion.

Statement. The necessary and sufficient condition for the existence
of the distribution P(x) is the convergence, for any history x0≥0 and
τ→0, of the solution x(t) to a fixed point.

Remark. Note the importance of the limit τ→0. As follows from
the analysis of the previous sections, P(x) can exist, though x(t) diverges
for some finite τ>0.

Different shapes of the distribution P(x), as a function of x, are shown
for the case of a single fixed point x∗1=0 in Fig. 18 and for either the
occurrence of the bistability region, with two fixed points x∗1=0 and
x∗3>0, or for the case of one nontrivial fixed point x∗2>0 in Fig. 19.

Figure 18: Distribution P(x), as a function of x, in the case of
either the bistability region with two fixed points x∗1=0 and x∗3>0
or for one nontrivial fixed point x∗2>0. The parameters are:
(a) σ1=σ2=−1, b=0.36, with x∗3=3.397, for α=0.5
(solid line), α=0.75 (dashed line), and α=1 (dashed-dotted line);
(b) σ1=σ2=−1, b=0.34, with x∗3=4.268, for the
same noise strengths as in (a); (c) σ1=σ2=−1, b=0.33,
with x∗3=4.671, for the same noise strengths as in (a);
(d) σ1=σ2=1, b=−1, with x∗2=0.567, for the
noise strengths α=0.125 (solid line), α=0.25 (dashed line), and
α=0.5 (dashed-dotted line).

Suppose, we make observations for the behavior of the function of interest x(t),
starting from an initial condition x0, during a short period of time, when this
function can be well characterized by the simple asymptotic form

x(t)≃x0+c1t+c2t2(t→0).

(49)

In real life, the coefficients ci can be defined from the observed data. And
for the considered equation, they are found by substituting expression (49)
into the evolution equation. Thus for the case τ>tc, we get

In the framework of self-similar factor approximants, the second-order factor
approximant reads

x∗(t)=x0(1+At)n.

(50)

The parameters A and n are defined by expanding (50) in powers of t
and comparing the expansion with the asymptotic form (49), which yields

A=σ2x0(1−bx0)exp(−bx0),n=σ1σ2−x0exp(−bx0)x0(1−bx0)exp(−bx0).

Let us recall that, as the analysis of the previous sections shows, for the
considered case of τ>tc, the finite-time singularity happens under one
of the following conditions: either for σ1=1,σ2=−1,b<0 and
for any history x0, or for σ1=σ2=−1,b<0 and x0>x∗2<e.
In both these cases, we find that A<0 and n<0, which allows us to
rewrite Eq. (50) as

x∗(t)=x0(1−|A|t)|n|.

(51)

The latter expression shows that the point of singularity is given by

tappc=1|A|=exp(bx0)σ2x0(1−bx0).

(52)

We have investigated the behavior of formula (52) for the different
situations studied in the previous sections. We find that, when the evolution
equation gives a solution x(t) diverging at a finite time, then the predicted
value (52) does approximate the real divergence point tc. When the
solution x(t) is bounded, approaching a stationary state, then either A
or n is positive, so that the factor approximant (50) does not predict
singularities. And, if the solution x(t) tends to infinity for t→∞,
then the factor approximant (50) either does not show a finite-time
singularity or, in some cases exhibits its appearance. Such artificial
singularities can be removed by constructing the factor approximants of higher
orders. In order not to complicate the consideration, here we limit ourselves
to the second-order factor approximant that does predict the singularity when
it really happens for x(t).

To estimate the accuracy of the prediction, we calculate the predicted tappc
for different values of the parameters and compare it with the tc given by the
evolution equation. The results in Fig. 20 demonstrate that the predicted singularity
point tappc is close to the real tc. The accuracy can be more precisely
characterized by the absolute and relative errors

Δ≡|tappc−tc|,ε≡Δtc×100%.

In Table 1, we present such error metrics, with fixed x0=1, for varying
parameters b.

Figure 19: Predicted singularity times tappc (red line with squares)
and real tc (blue line with circles), as functions of the parameter b=−|b|
and history x0, under σ1=1, σ2=−1 for: (a) varying |b|,
with fixed x0=1; (b) varying x0, with fixed b=−1.
\tbl

The singularity time tappc predicted by the self-similar factor
approximant (52), as compared to the exact singularity time tc following from
the evolution equation
btctappcε%Δ-0.00010.6921.0045%0.308-0.010.6540.98050%0.326-0.10.5010.82364%0.322-0.50.2490.40462%0.155-1.00.1240.18448%0.060-1.50.0650.08937%0.024-2.00.0350.04529%0.010-2.50.0190.02426%0.005-3.00.0100.01220%0.0020-3.50.00570.006718%0.0010-4.00.00320.003714%0.0005

Note that tappc is systematically larger than the true singularity
time tc, as can be expected from the fact that the sole information used
in the prediction is the quadratic asymptotic representation (49), which
necessarily underestimates the full strength of the nonlinear feedback leading
to the singularity. Taking third and higher order terms into account would
lead to significant improvement of the prediction accuracy. But we believe that
using the quadratic asymptotic form (49) is a realistic proxy for the
capture of early time dynamics in real life situations. While the relative
errors are significant (from 14% to 64% in the investigated cases), we believe
that these predictions are useful to provide an approximate estimation of the
critical time of the singularity. Taking into account more terms in the
asymptotic expansion for x(t) would improve the accuracy of the prediction,
however involving more complicated expressions for the critical time.

We have considered the evolution equation describing the population dynamics
with functional delayed carrying capacity. The linear delayed carrying capacity,
advanced earlier by the authors, has been generalized to the case of a nonlinear
delayed carrying capacity. This allowed us to treat the delayed feedback of
the evolving population on the capacity of their surrounding, by either creating
additional means for survival or destroying the available resources, when the
feedback can be of arbitrary strength. This is contrary to the linear approximation
for the capacity, which assumes weak feedback. The nonlinearity essentially
changes the behavior of solutions to the evolution equation, as compared to the
linear case.

The justification for the exponential form of the nonlinearity is based on the
derivation of an effective limit of expansion (5) for the carrying capacity by
invoking the self-similar approximation theory.

All admissible dynamical regimes have been analyzed, which happen to be of the
following types: punctuated unbounded growth, punctuated increase or punctuated
degradation to a stationary state, convergence to a stationary state with sharp
reversals of plateaus, oscillatory attenuation, everlasting fluctuations,
everlasting up-down plateau reversals, and divergence in finite time. The theorem
has been proved that, for the case of gain and competition, the solutions are
always bounded, when the feedback is destructive.

We have studied the influence of additive noise in two cases: (i) on the solutions
exhibiting finite-time singularities and (ii) in the presence of stationary solutions.
For the former case, we found that even a small noise level profoundly affects
the position of the finite-time singularities. For the later case, we have used
the Fokker-Planck equation and derived the general condition for the existence of
a stationary distribution function.

Finally, we showed that the knowledge of a simple quadratic asymptotic behavior
of the early time dynamics of a solution exhibiting a finite-time singularity
provides already sufficient information to predict the existence of a critical
time, where the solution diverges.

It is necessary to stress that taking into account the nonlinear delayed
carrying capacity not merely changes quantitatively the behavior of the solutions
to the evolution equation, but also removes artificial finite-time divergence
and finite-time death that exist in the equation with the linear form of the
carrying capacity. For example, the linear carrying capacity can lead to the
appearance of finite-time singularity or finite-time death even in the case of
prevailing competition (σ2=1), as is found in
[Yukalov et al.(2009), Yukalov et al.(2012a)]. But with the nonlinear carrying capacity, as used in
the present paper, these finite-time critical phenomena are excluded. Now,
finite-time singularity can occur only in the logically clear case of cooperation
(σ2=−1). The reason why the linear approximation for the carrying
capacity leads to such artificial singularities and deaths has been explained in
Sec. 2.1 of the present paper.

\nonumsection

Acknowledgments

Financial support from the ETH Competence Center ”Coping with Crises in Complex
Socio-Economic Systems” is appreciated.