**1. Introduction**

Currently, special attention is paid to important objects such as chains of interacting oscillators. Such chains arise when modeling many applied problems in radiophysics (see [1–3]), laser optics (see [4–6]), mechanics (see [7,8]), neural network theory (see [9–13]), biophysics (see [14]), mathematical ecology (see [15–18]), etc. This paper investigates the chains of coupled logistic equations with delay that are relevant for biophysics and mathematical ecology.

As a basic example describing certain population size changes, the well-known logistic equation with delay

$$
\dot{u} = r[1 - u(t - T)]u \tag{1}
$$

is considered. Here, *u*(*t*) > 0 is the normalized population or the population density, *r* > 0 is the Malthusian parameter, *T* > 0 is the delay associated with the reproductive age of individuals.

Let us recall some well-known facts (see, for example, [19]). In (1), the equilibrium state *u*<sup>0</sup> ≡ 1 is asymptotically stable as *rT π*/2, but this equilibrium state is unstable and there is a stable cycle *u*0(*t*) as *rT* > *π*/2. Under the condition

$$0 < rT - \frac{\pi}{2} \ll 1$$

**Citation:** Kashchenko, S. Quasinormal Forms for Chains of Coupled Logistic Equations with Delay. *Mathematics* **2022**, *10*, 2648. https://doi.org/10.3390/ math10152648

Academic Editor: Ferenc Hartung

Received: 3 July 2022 Accepted: 26 July 2022 Published: 28 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

its asymptotic behavior is:

$$u\_0(t) = 1 + \sqrt{rT - \frac{\pi}{2}} \left[ \xi\_0(\tau(1 + O(\varepsilon))) \exp\left(i\frac{2\pi}{T\_0}t\right) + 1 \right]$$

$$\bar{\xi}\_0(\tau(1 + O(\varepsilon))) \exp\left(-i\frac{2\pi}{T\_0}t\right) \Big| + \dots \dots$$

This cycle is relaxational as *rT* 1. Its asymptotic behavior is presented in [20]. We examine two types of chains of the coupled logistic equations with delay. Type 1: The chains of the 2*N* + 1 coupled equations

$$\dot{v}\_{\dot{\gamma}}(t) = r[1 - v\_{\dot{\gamma}}(t - T)]v\_{\dot{\gamma}}(t) + r \sum\_{m = -N}^{N} a\_{m\dot{\gamma}} \left( v\_m(t - h) - v\_{\dot{\gamma}}(t - h) \right). \tag{2}$$

Here, *j*, *m* = 0, ±1, ±2, ... , ±*N*, *amm* = 0, and *h* ≥ 0. Certain constraints may be imposed on the coupling coefficients *amj*. For example, the biologically derived conditions *uj*(*t*) ≥ 0 must be satisfied for *h* = 0. Therefore, the solutions must keep the quadrant *uj* ≥ <sup>0</sup> (*<sup>j</sup>* = 0, ±1, ±2, ... , ±*N*) invariant for nonnegative initial conditions *uj*(*s*) ∈ *<sup>C</sup>*[−*T*,0] (*<sup>j</sup>* = 0, ±1, ±2, . . . , ±*N*). Hence,

$$a\_{m\bar{j}} \ge 0.\tag{3}$$

System (2) has the equilibrium state *uj*<sup>0</sup> = *u*<sup>0</sup> = 1 (*j* = 0, ±1, ±2, ... , ±*N*). The substitution

$$
v\_{\rangle} = 1 + \mu\_{\rangle} \tag{4}$$

leads to the system of equations

$$\dot{u}\_{\dot{\gamma}}(t) = -r u\_{\dot{\gamma}}(t - T) + r \sum\_{m=-N}^{N} a\_{m\dot{\gamma}} \left( u\_m(t - h) - u\_{\dot{\gamma}}(t - h) \right) - r u\_{\dot{\gamma}}(t - T) u\_{\dot{\gamma}}.\tag{5}$$

Type 2: Here, the site of couplings is the only difference between the chains of coupled equations and (2):

$$
\psi\_j(t) = r \left[ 1 - \upsilon\_j(t - T) + \sum\_{m = -N}^{N} a\_{mj} \left( \upsilon\_m(t - h) - \upsilon\_j(t - h) \right) \right] \upsilon\_j(t). \tag{6}
$$

Substitution (4) converts this system to the system of equations

$$\begin{split} \dot{u}\_{j}(t) &= -r u\_{j}(t - T) + r \sum\_{m = -N}^{N} a\_{mj} \left( u\_{m}(t - h) - u\_{j}(t - h) \right) - \\ r u\_{j}(t - h) u\_{j} &+ r u\_{j}(t) \sum\_{m = -N}^{N} \left( u\_{m}(t - h) - u\_{j}(t - h) \right). \end{split} \tag{7}$$

The last term of System (7) distinguishes it from System (5).

It is convenient to associate the element *uj*(*t*) with the value of the two-variable function *u*(*t*, *xj*). Here, *xj* is the point of some circle with the angular coordinate *xj* = <sup>2</sup>*π*(2*<sup>N</sup>* + <sup>1</sup>)−<sup>1</sup> · *<sup>j</sup>*. Further, the periodicity condition *uj*+2*N*+1(*t*) = *<sup>u</sup>*(*t*, *xj* + <sup>2</sup>*π*) = *uj*(*t*) = *u*(*t*, *xj*) and *am*+2*N*+1,*<sup>j</sup>* = *am*,*j*+2*N*+<sup>1</sup> = *amj* (*m*, *j* = 0, ±1, ±2, . . .) holds.

We assume the coupling between elements to be homogeneous, i.e.,

$$a\_{m\circ} = a\_{m-j}.\tag{8}$$

The basic assumption is that the quantity (2*N* + 1) of elements is large enough:

$$0 < \varepsilon = 2\pi (2N + 1)^{-1} \ll 1. \tag{9}$$

The above condition gives ground for moving from the discrete spatial variable *x* to the continuous variable *x*.

In this paper, the three most important cases in terms of applications will be considered. In the first case, it is assumed that the system is fully connected and all coupling coefficients are the same: the equalities

$$a\_{\rm mj} = \gamma (2N + 1)^{-1} (2\pi)^{-1} \tag{10}$$

hold for some coefficient *γ*. The expression *γ* 2*π*(2*N* + 1) <sup>−</sup><sup>1</sup> *<sup>N</sup>* ∑ *j*=−*N f*(*xj*) is the partial

Darboux sum for the expression *γ*(2*π*)−<sup>1</sup> 2*π* 0 *f*(*x*)*dx*. Thus, under condition (10), the boundary value problem

$$\begin{split} \frac{\partial u}{\partial t} &= -ru(t - T, \mathbf{x}) + r\gamma \left[ (2\pi)^{-1} \int\_0^{2\pi} u(t - h, s) ds - u(t - h, \mathbf{x}) \right] - \\ &\quad ru(t - T, \mathbf{x})u + r\delta\gamma u \cdot \left[ (2\pi)^{-1} \int\_0^{2\pi} u(t - h, s) ds - u(t - h, \mathbf{x}) \right] \end{split} \tag{11}$$

is the asymptotic approximation for Systems (5) and (7) as *ε* → 0. Here, *δ* = 0 in the case of (5), and *δ* = 1 in the case of (7).

For Equation (11), the periodic boundary conditions

$$
\mu(t, x + 2\pi) \equiv \mu(t, x) \tag{12}
$$

hold. Note that in the case when coefficients *ai*−*<sup>j</sup>* are not identical, the integral term becomes more complicated. We obtain the expression <sup>1</sup> 2*π* <sup>2</sup>*<sup>π</sup>* <sup>0</sup> *a*(*s*)*u*(*t* − *h*, *x* + *s*)*ds* instead of the expression <sup>1</sup> 2*π* <sup>2</sup>*<sup>π</sup>* <sup>0</sup> *<sup>u</sup>*(*<sup>t</sup>* <sup>−</sup> *<sup>h</sup>*,*s*)*ds*, and <sup>1</sup> 2*π* <sup>2</sup>*<sup>π</sup>* <sup>0</sup> *a*(*s*)*ds* = 1.

The local (i.e., in the zero equilibrium state neighborhood) dynamics of the boundary value problem (11) and (12) are studied in Section 2.

Then, we consider the case of so-called diffusional coupling, where

$$a\_1 = a\_{-1} = \gamma \text{ and } a\_k = 0 \text{ for } k \neq \pm 1. \tag{13}$$

Now, from systems (5) and (7), we arrive at the boundary value problem

$$\begin{split} \frac{\partial u}{\partial t} &= -ru(t - T, \mathbf{x}) + r\gamma \left[ u(t - h, \mathbf{x} + \varepsilon) - 2u(t - h, \mathbf{x}) + u(t - h, \mathbf{x} - \varepsilon) \right] - \\ &\quad ru(t - T, \mathbf{x})u(t, \mathbf{x}) + r\delta\gamma \left[ u(t - h, \mathbf{x} + \varepsilon) - \\ &\quad 2u(t - h, \mathbf{x}) + u(t - h, \mathbf{x} - \varepsilon) \right] u(t, \mathbf{x}), \quad u(t, \mathbf{x} + 2\pi) \equiv u(t, \mathbf{x}). \end{split} \tag{14}$$

If the coefficients *ak* (*k* = ±1) do not differ much from zero, the dynamic behavior of (14) could change significantly. Thus, we examine a more general (compared to (14)) boundary value problem

$$\begin{split} \frac{\partial u}{\partial t} &= -r u(t - T, x) - r u(t - T, x) u(t, x) + \\ &\quad \qquad r \gamma \left( \int\_{-\infty}^{+\infty} \left( F\_t(s) - 2F\_0(s) + F\_{-t}(s) \right) u(t - h, x + s) ds \right) \times \\ &\quad \quad \left[ 1 + \delta u(t, x) \right], \quad u(t, x + 2\pi) \equiv u(t, x). \end{split} \tag{15}$$

,

Here,

$$\begin{array}{rcl} F\_{\pm\varepsilon}(s) & = & \frac{1}{\sigma\sqrt{2\pi}} \exp\left[ - (2\sigma^2)^{-1} (s \pm \varepsilon)^2 \right], \\ F\_0(s) & = & \frac{1}{\sigma\sqrt{2\pi}} \exp\left[ - (2\sigma^2)^{-1} s^2 \right]. \end{array}$$

Undoubtedly, the integral expression in (15) can be presented in the form of integrals from 0 to 2*π* of some function *F*˜(*s*). However, the form of (15) is preferable. Firstly, the asymptotic equality

$$\int\_{-\infty}^{+\infty} \left( F\_{\varepsilon}(s) - 2F\_{0}(s) + F\_{-\varepsilon}(s) \right) w(\mathbf{x} + s) ds = w(\mathbf{x} + \varepsilon) - 2w(\mathbf{x}) + w(\mathbf{x} - \varepsilon) + o(1)$$

holds for the fixed function *w*(*x*) as *δ* → 0. Secondly, it is technically convenient to calculate the importance of further usage of integrals explicitly. For example,

$$\int\_{-\infty}^{+\infty} F\_{\varepsilon}(s) \exp(iks) ds = \exp\left(ik\varepsilon - \frac{1}{2}\delta^2(\varepsilon k)^2\right). \tag{16}$$

In Section 3, the boundary value problem (15) is examined.

According to the equalities in (13), interactions with one 'neighbor' on the right and one 'neighbor' on the left make a dominating contribution to the couplings between elements as *ε* → 0. In the next, fourth, section, the case of coupling coefficient interactions for unidirectional couplings is studied. For *a*<sup>1</sup> = *γ*, *aj* = 0 as *j* = 1, the resulting boundary value problem assumes the form

$$\begin{array}{lcl}\frac{\partial u}{\partial t} &=& -ru(t-T,\mathbf{x}) - ru(t-T,\mathbf{x})u(t,\mathbf{x}) + \\ & r\gamma \Big(\int\_{-\infty}^{+\infty} \left(F\_{\mathbf{t}}(\mathbf{s}) - F\_{\mathbf{0}}(\mathbf{s})\right) u(t-h,\mathbf{x}+\mathbf{s})d\mathbf{s}\Big) \cdot \Big[1+\delta u(t,\mathbf{x})\Big]\_{\mathbf{y}} \\ & u(t,\mathbf{x}+2\pi) \equiv u(t,\mathbf{x}). \end{array} \tag{17}$$

The paper is devoted to the investigation of the behavior of the solutions to the boundary value problems (11), (12); Equations (15) and (17) with initial conditions from some sufficiently small metric *<sup>C</sup>*[−*T*,0] × *<sup>W</sup>*[0,2*π*] (*T*<sup>0</sup> = max(*T*0, *<sup>h</sup>*)) *<sup>ε</sup>*-independent zero equilibrium state neighborhood for small *ε*. In each of the problems, critical cases are identified in the study of the stability of the zero solution. In critical cases, the characteristic equations of the linearized-at-zero problems have infinitely many roots, with real parts tending to zero as *ε* → 0. This is a distinctive feature of the problems under consideration. Thereby, we may speak of the infinite-dimensional critical cases implementation. The standard methods of integral manifolds (see, for example, [21–23]) and normal forms (see, for example, [24]) are not directly applicable for their study. Thus, we employ the methods developed by the author in [25–28]. Their essence lies in the construction of special first-approximation equations called quasinormal forms (QNFs), whose nonlocal dynamics determine the local structure of the initial boundary value problems' solutions. These QNFs are nonlinear boundary value problems of neutral or parabolic types with either one or two spatial variables. The solutions of these QNFs make it possible to determine the principal terms of the asymptotic representations of the considered boundary value problems' solutions.

Section 5 is a natural extension of Section 4. It studies the model of a chain with a unidirectional coupling under the condition that the coupling coefficient between elements takes sufficiently large values. It is about the boundary value problem

$$\dot{N} = r[1 - N(t - T, \mathbf{x})]N + \gamma \left[ \int\_{-\infty}^{\infty} F(s)N(t, \mathbf{x} + s)ds - N \right],\tag{18}$$

$$N(t, \ge +2\pi) \equiv N(t, \ge),\tag{19}$$

where *γ* 1. We note that (18) is interpreted as a logistic delay equation with spatially distributed control. The function *F*(*s*) describing spatial interactions is given by

$$F(s) = \frac{1}{\sqrt{\mu \pi t}} \exp[-\mu^{-1}(s+h)^2], \quad \mu > 0.$$

We demonstrate that the dynamic properties of the problem under consideration change significantly depending on the various relations between *γ*, *μ*, and *h*.

We present without proof two analogs of the classical Lyapunov stability theorems in the first approximation statements.

For the linearized-at-zero Equations (11), (15), (17) and (18) with periodic boundary conditions (12), the characteristic equation has the form

$$\lambda + r \exp(-\lambda) = d \left( \int\_{-\infty}^{\infty} F(s, \varepsilon) \exp(iks) ds - 1 \right), \quad k = 0, \pm 1, \pm 2, \ldots \tag{20}$$

**Statement 1.** *Let the roots of Equation* (20) *have a negative real part and be separated from the imaginary axis as ε* → 0*. Then, an ε*<sup>0</sup> > 0 *can be found such that, for ε* ∈ (0,*ε*0)*, the solutions to the problem under consideration from the sufficiently small ε-independent zero equilibrium state neighborhood tend to zero as t* → ∞*.*

**Statement 2.** *Let Equation* (20) *have a root with a positive real part separated from zero as ε* → 0*. Then, the zero equilibrium state of the boundary value problem is unstable for small ε, and there is no attractor of this boundary value problem in some sufficiently small ε-independent zero neighborhood.*

Thus, in the posed problem about local conditions in the zero equilibrium state neighborhood, the dynamics are trivial under the condition of Statement 1, and it cannot be studied by local analysis methods under the condition of Statement 2.
