The Lynx and Hare Data of 200 Years as the Nonlinear Conserving Interaction Based on Noether’s Conservation Laws and Stability

Hiroshi Uechi¹, Lisa Uechi², Schun T. Uechi³

¹Osaka Gakuin University, Kishibe-minami, Osaka, Japan.
²University of California, Los Angeles, USA.
³KPMG Ignition Tokyo, Data Tech., Tokyo, Japan.

DOI: 10.4236/jamp.2021.911181                

Abstract

We applied n-variable conserving nonlinear differential equations (n-CNDEs) to the population data of the 10-year cycles of Canadian lynx (1821-2016) and the snowshoe hare (1845-1921). Modeling external effects as perturbations to population dynamics, recovering and restorations from disintegrations (or extinctions), stability and survival strategies are discussed in terms of the conservation law inherent to dynamical interactions among species. The 2-variable conserving nonlinear interaction (2CNIs) is extended to 3, 4, ... n-variable conserving nonlinear interactions (n-CNIs) of species by adjusting minimum unknown parameters. The population cycle of species is a manifestation of conservation laws existing in complicated ecosystems, which is suggested from the CNDE analysis as a standard rhythm of interactions. The ecosystem is a consequence of the long history of nonlinear interactions and evolutions among life-beings and the natural environment, and the population dynamics of an ecosystem are observed as approximate CNIs. Physical analyses of the conserving quantity in nonlinear interactions would help us understand why and how they have developed. The standard rhythm found in nonlinear interactions should be considered as a manifestation of the survival strategy and the survival of the fittest to the balance of biological systems. The CNDEs and nonlinear differential equations with time-dependent coefficients would help find useful physical information on the survival of the fittest and symbiosis in an ecosystem.

Keywords

The 10-Year Population Cycles of Canadian Lynx and Snowshoe Hare, The Standard Rhythm of Interactions, Noether’s Conservation Law, Conserving Nonlinear Interactions of Species

1. Introduction

Ecology is a field of study on complex many-body interactions among organisms, living beings and natural environment. The conservation status of a group of living animals is measured in many countries, whether certain species still exist or how likely to become extinct in the near future. In reality, human activities from political, social and cultural conduct to agricultural, engineering procedures caused serious problems for the conservation of species, protection of biological diversity in the scale of the earth. The development of human societies means to thrive against the severe natural environment by utilizing natural resources for wellness of human beings and society, but environmental pollutions and contaminations caused by disposed chemical substances affect living beings. The food chain among life is extremely complex that human activities affect the food chain, natural environment and biodiversity. It is important to understand interactions of vast fields of physical world among chemical substances, natural environment and ecology in terms of biodiversity of plants and animals, dynamics among microbes, pathogens and viruses. The carbon, nitrogen and water cycles on the level of global environment are essential, and the balance of cycles will comprise healthy environment to sustain life on the earth.

The Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES, 2021) [1] is a multilateral international treaty to protect endangered species, which has goals for conservation of biodiversity, environment and natural resources, including protection and management of living habitat. In order to understand the importance of biodiversity and ecology, it is essential to collect specific population data and the natural habitat of endangered species. Any population data would reveal interactions of species and natural environment, and nonlinear differential equations have been applied to population changes among 2 - 4-variable factors of microbes, food chains and interacting species [2] [3] [4]. Though population changes can be reasonably expressed by nonlinear differential equations with data accumulations and computer simulations, it is imperative for us to extract physical reasons and dynamics of nonlinear interactions in ecological systems.

Nonlinear differential equations provide us with a possible account for complex change of populations by adjusting coupling constants, but they do not explain physical principles and reasons why such nonlinear interactions and coupling constants are physically chosen by nature, and in this sense, dynamical change of ecological system is not simply solved by computer simulations and big-data accumulations. Hence, we have proposed nonlinear conserving equations for n-body interactions [5] [6] [7] and applied them to population data of species and microbes in order to understand and extract useful information of what the population growth and survival of species would signify. The population cycles in ecosystems and those found in self-organization of biological systems have similar periodic oscillations, and we found a characteristic dynamical rhythm [6] [7] expressed by conserving nonlinear interactions (CNIs). Useful concepts such as the survival of the fittest and relations to conservation laws manifest in nonlinear interactions are gradually recognized when a rhythm of big data is analyzed by the model of conserving nonlinear differential equations (CNDEs).

The conserving nonlinear equation is applied dynamical interactions of species such as the population regulation in Canadian lynx and snowshoe hare [8] [9] [10] [11], the food-web of Microbes in Okanagan Lake [12]. Characteristic changes in population data, the relations between initial conditions and nonlinear interactions are observed in CNDEs, which restricts physically meaningful nonlinear forms of interactions. One should notice that the physically meaningless nonlinear interactions known as atto-fox problems [13] [14] [15] [16] [17] are strictly restricted in the conserving nonlinear equations. The CNDE denotes the strenuous correlation between an initial condition and a conservation law, which may designate symbiosis and the survival of the fittest of a group of species. A dynamical mechanism of biome on the earth demands knowledge of vast fields of science and ecology, such as microbiological dynamics through cycles of organic chemical substances [18] [19] [20], the food chain [21] [22], environmental changes and natural resources such as air, water and soil, resulting in increasing cooperative studies of intersecting fields [23] [24] [25] [26]. The interactions of research in many fields would develop and enhance our understandings of concepts of life, diversity and ecology.

Problems of population dynamics such as termination, extinction [27] [28] or restoration of species are important for human society to thrive on the earth. Though interactions among species could be expressed by nonlinear differential equations, scientific and physical meanings must be yet investigated. The functional form of nonlinearity is generally discussed in [5] [6] [7], and in the current analysis, the nonlinear polynomial of degree two is applied with the variational method (Noether’s theorem) and conservation law to examine dynamics from big-data. There are two types of nonlinear differential equations: nonconserving, dissipative nonlinear equations, and conserving nonlinear equations which have the conserved quantity, CNIs and equations are discussed, and symbiosis, the food-chain and the patterns of interaction, the survival of the fittest of consisting species are discussed as much as possible by comparing conserving nonlinear equations with conventional nonlinear differential equations.

The Canadian lynx is a lynx species native to North America, inhabiting from Canada to Alaska and the northern United States. The Canadian lynx is reported to be dependent heavily on snowshoe hares for food, which has been studied by many researchers for more than 200 years. These hares comprise 35% - 97% of lynx food, but the diet varies by abundance of hares and the season. Lynxes consume other animals, such as ducks, grouse, ptarmigan, red squirrels, voles and moles, young sheep, mule deer and reindeer, but snowshoe hares are still the primary chief component, 600 g - 1200 g of food every day [29]. The pray-predator cycles, rises and falls of lynx populations are strongly correlated o those of hare population cycles, consequently known as the 10-year cycle of population density of Canadian lynx and snowshoe hare; a period of hare scarcity occurs every 8 to 11 years [30]. The abundance of lynx has been estimated from the data of Canada lynx fur-trades return of the Hudson’s Bay Company, since the 1730s [31] [32]. A cycle of its abundance is characterized by huge rises and falls of lynx populations with the peaks typically ten times higher than the troughs, resulting in population changes in other related species. We examined the lynx-hare data for about 200 years for 1821-2016 supplied by courtesy of Government of Northwest Territories [33], which is numerically simulated and analyzed by the model of conserving nonlinear differential equations.

The CNI generated a characteristic stable population cycle, and we termed it as the standard rhythm, which is a consequence of conservation laws maintained by interacting species. The standard rhythm may be related to the concept of the survival of the fittest and symbiosis among species, which should be investigated with big data in the various fields of ecology. The conserving nonlinear model suggests that biomechanical and microbiophysical systems show symbiosis among cells and organisms, characteristic stability and balance [34] [35] and restoration phenomena. The standard and stable rhythm of population density could be interpreted as a manifestation of the survival of the fittest to the balance of nonlinear dynamical interactions. The conservation laws and their physical meanings are neither clearly understood nor investigated yet in many fields of macro- and micro-biophysics.

The basic property of 2-variable conserving dynamical interactions, Noether’s theorem and the conserved quantity, Ψ-function, is reviewed in Section 2 and properties of Ψ-function are explained in detail in Section 3. Sections 2 and 3 are brief reviews to understand conserving differential equations. The new approach based on conserving nonlinear interactions and extensions to (n + m)-variable interactions are discussed in Section 4. The indications and inferences of CNIs for 10-year population cycles in 200 years (1821-2016 years) of Canadian Lynx-Snowshoe hare are discussed in Section 5. We applied the CNI method to Lynx-Hare data supplied from the Environment and Natural Resources 2017 Tundra Ecological Research Station Small Mammal Database. The difference between Ψ-function and Hamiltonian interpretation, the relation between Noether’s theorem and Lyapunov function are discussed in Section 6. Hence, the new results and approaches, applications to Lynx-Hare survival strategies are discussed from Section 4 to Section 6 in detail. The erroneous notion of the atto-fox problems is strictly criticized and debunked. Comments on the survival of the fittest and symbiosis based on conserving approximations are concluded in Section 7. The conserving nonlinear interactions compared with nonconserving interactions would help understand a fundamental concept of survivals and evolutions of ecological systems.

2. The 2-Variable CNIs, Restoration from External Perturbations

Nonlinear equations are classified into non-conserving dissipative equations suitable for transient behaviors and phenomena short in time and conserving nonlinear equations suitable for life living-long and surviving against external destructive perturbations, which corresponds to irreversible dissipative phenomena and thermomechanical equilibrium [36]. Though life-involving phenomena are very different from those of physical materials and particles, viruses and cells, plants and animals are not against physical laws and principles; in other words, they are making good use of natural laws for their existence.

The Lotka-Volterra types of nonlinear differential equations [37] [38] are primitive nonlinear equations and have been applied to study competitive and predator-pray interactions, mathematical models of disease spreading among species, pest-control, evolving ecosystem networks, population and epidemiology and so forth [39] [40] [41]. The Lotka-Volterra types of nonlinear differential equations are easy to use, but they have intrinsic blunder known as the atto-fox problem, which is prevented in the conserving nonlinear equations [5] [6] [7]. Since the Lynx-hare cycle is an approximate nonlinear conserving phenomena, it has a conserved quantity, Ψ-function which restricts unphysical atto-fox problem.

The 2n-variable conserving nonlinear differential equations give the following useful results.

1) The form of differential equations and coefficients of nonlinear interactions are strictly confined with initial conditions when a population of species is stable for a long period of time, designating a conservation law.

2) The conserved quantity Ψ-function produces a Lyapunov function usually employed to study solutions to nonlinear differential equations, which is helpful to study non-conserving dissipative interactions.

3) A complex interacting system can be approximately decomposed into an assembly of binary-coupled forms (BCF).

4) The binary coupled system with the conservation law indicates an addition law interpreted as the restoration or rehabilitation phenomena, such as a small damaged device replaced by a new one.

5) The conservation law is also useful to check accuracy of numerical solutions to conserving nonlinear differential equations.

The purpose of conserving or dissipative nonlinear analyses among species is to study stability and restoration, prosperity and degeneration mechanism of ecological systems, effects on natural or artificial environmental changes in dynamics of ecology. It should help superintend proper procedures and methods to sustain natural and ecological systems in the future based on data in the past.

2.1. A Brief Review of Noether’s Theorem and Conserved Quantities

The necessary condition for extrema in Lagrangian formulation, $\mathcal{L}\left(t\mathrm{,}{x}^k\left(t\right)\mathrm{,}{\stackrel{\dot{}}{x}}^k\left(t\right)\right)$, is given by

$\delta J=\delta {\displaystyle \int \mathcal{L}\left(t\mathrm{,}{x}^k\left(t\right)\mathrm{,}{\stackrel{\dot{}}{x}}^k\left(t\right)\right)\text{d}t}=0\left(k=\mathrm{1,}\cdots \mathrm{,}n\right)\mathrm{,}$ (2.1)

and all functions, $x\left(t\right)=\left(x^1\left(t\right),\cdots ,x^n\left(t\right)\right)$, $t\in \left[a,b\right]$, belong to $C^2\left[a\mathrm{,}b\right]$, which denotes the set of all continuous functions on the interval $\left[a\mathrm{,}b\right]$ and the second derivatives of all functions are continuous. If $x\left(t\right)$ is a relative minimum of the functional J, the condition (2.1) generates,

$E_k\equiv \frac{\partial L}{\partial x^k}-\frac{\text{d}}{\text{d}t}\left(\frac{\partial \mathcal{L}}{\partial {\stackrel{\dot{}}{x}}^k}\right)=0.$ (2.2)

This is the Euler-Lagrange equation which determines the equations of motion of a system.

Noether’s theorem describes that the Lagrangian with Euler-Lagrange equation is invariant under certain space-time transformations, and invariance of Lagrangian generates respective conservation laws [42] [43]. For example, the time-translation invariance of Lagrangian corresponds to the conservation of energy of a system. Let us consider r-parameter transformations in general that will be regarded as transformations of configuration space, $\left(t\mathrm{,}{x}^1\mathrm{,}\cdots \mathrm{,}{x}^n\right)$ -space, depending upon r real, independent parameters ${\epsilon }^1\mathrm{,}\cdots \mathrm{,}{\epsilon }^r$. The transformations are defined by

$\begin{array}{l}\bar{t}=t+{\tau }_s\left(t\mathrm{,}x\right){\epsilon }^s+o\left(\epsilon \right)\mathrm{,}\\ {\bar{x}}^k=x^k+{\xi }_s^k\left(t\mathrm{,}x\right){\epsilon }^s+o\left(\epsilon \right)\mathrm{,}\end{array}$ (2.3)

where s ranges over $\mathrm{1,}\cdots \mathrm{,}r$, and $o\left(\epsilon \right)$ denotes the terms which go to zero faster than $\left|\epsilon \right|$, ${\lim }_{\left|\epsilon \right|\to 0}o\left(\epsilon \right)/\left|\epsilon \right|=0$. The functions ${\tau }_s\left(t\mathrm{,}x\right)$ and ${\xi }_s^k\left(t\mathrm{,}x\right)$ for linear parts of $\bar{t}$ and ${\bar{x}}^k$ with respect to $\epsilon$ are commonly called the infinitesimal generators of the transformations. Classical theorem of Emmy Noether on invariant variational problems can be obtained under the hypotheses of extrema (2.1), and transformation (2.3). The result is the r identities produced by the transformation (2.3),

$-E_k\left({\xi }_s^k-{\stackrel{\dot{}}{x}}^k{\tau }_s\right)=\frac{\text{d}}{\text{d}t}\left(\left(\mathcal{L}-{\stackrel{\dot{}}{x}}^k\frac{\partial \mathcal{L}}{\partial {\stackrel{\dot{}}{x}}^k}\right){\tau }_s+\frac{\partial \mathcal{L}}{\partial {\stackrel{\dot{}}{x}}^k}{\xi }_s^k-{\Phi }_s\right),$ (2.4)

where $k=1,\cdots ,n$ is summed, and ${\Phi }_s$ is an arbitrary function defined as a gauge function of the transformation. Note that the arbitrary choice of a gauge function will not change equations of motion of a system, which can be usually used for a convenient expression of conservation laws. If the fundamental integral (2.1) is divergence-invariant under the r parameter group of transformation (2.3), and if $E_k=0$ for $k=1,\cdots ,n$, then following r expressions are obtained:

${\Psi }_s\equiv \left(\mathcal{L}-{\stackrel{\dot{}}{x}}^k\frac{\partial \mathcal{L}}{\partial {\stackrel{\dot{}}{x}}^k}\right){\tau }_s+\frac{\partial \mathcal{L}}{\partial {\stackrel{\dot{}}{x}}^k}{\xi }_s^k-{\Phi }_s=\text{constant}.$ (2.5)

This defines the conserved quantities of a system. Since the expressions Ψ defined in (2.5) are constant with the condition: $E_k=0\left(k=1,\cdots ,n\right)$, they are the first integrals of the differential equations of motion. In physical applications, the first integral (2.5) is interpreted as the energy of the system, whose governing equations are $E_k=0$. In general, the function Ψ is constant with respect to time. By employing the formalism, the conservation law and symmetry of nonlinear differential equations for coupled conserving systems are discussed.

The meaning of a conservation law in biological complex systems may be very different from conservation laws in physics described by way of a Lagrangian or a Hamiltonian. It is difficult to directly use mechanical concepts and analogies such as potential and kinetic energy, but it should be useful to employ the concept of conservation law and symmetry in terms of Noether’s theorem in order to investigate conservation laws corresponding to dynamics of biodiversity.

2.2. Nonlinear Conserving 2-Variable Interactions

It may seem natural to assume that a nonlinear system eventually dissipates energy and arrives at an equilibrium state, and the equilibrium state would be considered as a maximum entropy state. If the state is to be activated, it needs some external inputs of energy, designating external inputs or perturbations. They are not usually included in the analysis of Lotka-Volterra and Kolmogorov type nonlinear equations. There are numerous examples of stable nonlinear dynamical systems from microbiological systems such as DNA and RNA, amino acids and proteins, cells and organs, to macroscopic systems such as the 10-year cycle of Canadian lynx and snowshoe hare, wolves and caribous, creatures in marine coral reefs. The stability suggests that some conservation laws and symmetries be inherent from complex microbiological systems to macroscopic ecological systems. Hence, it could be possible to investigate conservation laws corresponding to biological stability based on Noether’s theorem in dynamical systems [5] [6] [42] [43] [44] [45] [46]. The 2-variable CNI is briefly reviewed in order to extend and construct solutions to higher-variable conserving nonlinear interactions.

The Lagrangian of 2-variable CNI with external perturbations is described as [6]:

$\begin{array}{c}\mathcal{L}={\alpha }_1{\stackrel{\dot{}}{x}}_1{x}_2+{\alpha }_2{x}_1{\stackrel{\dot{}}{x}}_2+{\alpha }_3{x}_1^2+\left({\alpha }_4+{\alpha }_5\right)x_1{x}_2\\ +{\alpha }_6{x}_2^2+{\alpha }_7{x}_1^2{x}_2+{\alpha }_8{x}_1{x}_2^2+c_2{x}_1+c_1{x}_2.\end{array}$ (2.6)

From (2.6), we get the following nonlinear differential equation,

${\stackrel{\dot{}}{x}}_1=\frac{1}{d_{21}}\left\{\left({\alpha }_4+{\alpha }_5\right)x_1+2{\alpha }_6{x}_2+2{\alpha }_8{x}_1{x}_2+{\alpha }_7{x}_1^2\right\}+\frac{c_1}{d_{21}},$ (2.7)

${\stackrel{\dot{}}{x}}_2=\frac{1}{d_{12}}\left\{2{\alpha }_3{x}_1+\left({\alpha }_4+{\alpha }_5\right)x_2+2{\alpha }_7{x}_1{x}_2+{\alpha }_8{x}_2^2\right\}+\frac{c_2}{d_{12}},$ (2.8)

and the 2-variable conserving nonlinear model has the following constant function in time, calculated by (2.5):

$\Psi \left(x_1\mathrm{,}{x}_2\right)={\alpha }_3{x}_1^2+\left({\alpha }_4+{\alpha }_5\right)x_1{x}_2+{\alpha }_6{x}_2^2+{\alpha }_7{x}_1^2{x}_2+{\alpha }_8{x}_1{x}_2^2+c_2{x}_1+c_1{x}_2\mathrm{.}$ (2.9)

The parameter $d_{21}$ is given by $d_{21}={\alpha }_2-{\alpha }_1=-d_{12}$. The Ψ-function does not have any velocity (time-derivative) terms, and it looks like potential energy terms. The interpretation of the function is an open question, depending on nonlinear systems.

2.3. Exponential Types of 2 Variable Coupled LV System

The Lyapunov function discussed in nonlinear problems can be obtained by Ψ-function derivation, which is shown in a simple example. The classical type of LV equation is defined as:

$\begin{array}{l}{\stackrel{\dot{}}{x}}_1=a_{11}{x}_1+a_{12}{x}_1{x}_2,\\ {\stackrel{\dot{}}{x}}_2=a_{21}{x}_2+a_{22}{x}_1{x}_2,\end{array}$ (2.10)

and the Equation (2.10) can be reproduced by the following Lagrangian,

$\mathcal{L}=\exp \left\{{\alpha }_1{x}_1+{\alpha }_2{x}_2\right\}\left\{{\alpha }_3{\stackrel{\dot{}}{x}}_1+{\alpha }_4{\stackrel{\dot{}}{x}}_2+{\alpha }_5{x}_1{x}_2\right\}\mathrm{.}$ (2.11)

Conventionally, the Lotka-Volterra equations are usually written in the form,

$\begin{array}{l}{\stackrel{\dot{}}{x}}_1=\frac{1}{d_{12}}{\alpha }_5\left\{{\alpha }_2{x}_1{x}_2+x_1\right\},\\ {\stackrel{\dot{}}{x}}_2=\frac{1}{d_{21}}{\alpha }_5\left\{{\alpha }_1{x}_1{x}_2+x_2\right\},\end{array}$ (2.12)

where $d_{12}={\alpha }_1{\alpha }_4-{\alpha }_2{\alpha }_3=-d_{21}$. The conserved quantity Ψ of (2.12) is given by

$\Psi \equiv {\alpha }_5\exp \left\{{\alpha }_1{x}_1+{\alpha }_2{x}_2\right\}{x}_1{x}_2\mathrm{.}$ (2.13)

The logarithm of (2.13) is similar to a well known Lyapunov function $V\left(x_1\mathrm{,}{x}_2\right)$ which has the property $\text{d}V\left(x_1\mathrm{,}{x}_2\right)/\text{d}t=0$. It proves that the time derivative of conserved quantity vanishes, $\text{d}\Psi /\text{d}t=0$, and therefore, Lyapunov function can be derived from Noether’s theorem. The result will be helpful to understand a system of nonlinear differential equations in physical terms. There are two main types of Lyapunov functions that are strict Lyapunov and non-strict Lyapunov functions [47]. Our conserved quantity would correspond to the strict Lyapunov function, due to the global property of Lagrangian approach. The proof of equivalence would relate stability and dynamics with the conservation law of a system.

The solutions $x_1$ and $x_2$ are rapidly changing in time shown in Figure 1(a), whose parameters and initial starting values are explained in [5] [6]. One should be careful that solutions are sensitive to nonlinear coupling constants and initial starting values, which should be chosen carefully for numerical simulations with external perturbations. The closed line in the phase space solutions shows a definite periodic time. Though solutions ( $x_1\mathrm{,}{x}_2$ ) are rapidly changing in time, but the conserved function, Ψ, is constant in time. The Ψ-function diverges or cannot be constant with arbitrary parameters. The eight interaction terms on the right-hand sides of (2.7) and (2.8) require eight independent parameters ( $d_{12}$ and $d_{21}$ are common factors and irrelevant). Since one can regard ( ${\alpha }_4+{\alpha }_5$ ) as one parameter, we have only five independent parameters in the beginning in

(a)(b)(c)

Figure 1. The 2-variable conserving solution and the constant function in time, $\Psi \left(x_1\left(t\right)\mathrm{,}{x}_2\left(t\right)\right)$. (a) The solution to the 2-variable conserving nonlinear system. Solid and dashed lines represent $x_1$ (prey) and $x_2$ (predator), respectively. (b) Phase-space of the 2-variable conserving solution. (c) The conservation law, $\Psi \left(x_1\left(t\right)\mathrm{,}{x}_2\left(t\right)\right)$, of the 2-variable conserving solution. It is constant with respect to time.

order to solve (2.7) and (2.8), which is helpful for numerical simulations of conserving nonlinear differential equations.

2.4. The 3-Variable Conserving Nonlinear Interactions

The conserving nonlinear interaction of 3-variable system will be produced by the Lagrangian of the following type:

$\begin{array}{c}\mathcal{L}={\alpha }_1{\stackrel{\dot{}}{x}}_1{x}_2+{\alpha }_2{\stackrel{\dot{}}{x}}_2{x}_3+{\alpha }_3{\stackrel{\dot{}}{x}}_3{x}_1+{\alpha }_4{x}_1{x}_2+{\alpha }_5{x}_1{x}_3+{\alpha }_6{x}_2{x}_3+{\alpha }_7{x}_1^2\\ +{\alpha }_8{x}_2^2+{\alpha }_9{x}_3^2+{\alpha }_{10}{x}_1^2{x}_2+{\alpha }_{11}{x}_1{x}_2^2+{\alpha }_{12}{x}_1^2{x}_3+{\alpha }_{13}{x}_1{x}_3^2\\ +{\alpha }_{14}{x}_2^2{x}_3+{\alpha }_{15}{x}_2{x}_3^2+c_3{x}_1+c_1{x}_2+c_1{x}_3\mathrm{,}\end{array}$ (2.14)

and the conserved quantity of (2.14) is given by

$\begin{array}{c}\Psi ={\alpha }_4{x}_1{x}_2+{\alpha }_5{x}_1{x}_3+{\alpha }_6{x}_2{x}_3+{\alpha }_7{x}_1^2+{\alpha }_8{x}_2^2+{\alpha }_9{x}_3^2+{\alpha }_{10}{x}_1^2{x}_2+{\alpha }_{11}{x}_1{x}_2^2\\ +{\alpha }_{12}{x}_1^2{x}_3+{\alpha }_{13}{x}_1{x}_3^2+{\alpha }_{14}{x}_2^2{x}_3+{\alpha }_{15}{x}_2{x}_3^2+c_3{x}_1+c_1{x}_2+c_1{x}_3.\end{array}$ (2.15)

The Lagrangian (2.14) produces the following coupled nonlinear differential equations:

$\begin{array}{l}{\alpha }_1{\stackrel{\dot{}}{x}}_2-{\alpha }_3{\stackrel{\dot{}}{x}}_3=2{\alpha }_7{x}_1+{\alpha }_4{x}_2+{\alpha }_5{x}_3+2{\alpha }_{10}{x}_1{x}_2+{\alpha }_{11}{x}_2^2+2{\alpha }_{12}{x}_1{x}_3+{\alpha }_{13}{x}_3^2+c_3,\\ {\alpha }_2{\stackrel{\dot{}}{x}}_3-{\alpha }_1{\stackrel{\dot{}}{x}}_1={\alpha }_4{x}_1+2{\alpha }_8{x}_2+{\alpha }_6{x}_3+2{\alpha }_{11}{x}_1{x}_2+2{\alpha }_{14}{x}_2{x}_3+{\alpha }_{10}{x}_1^2+{\alpha }_{15}{x}_3^2+c_1,\\ {\alpha }_3{\stackrel{\dot{}}{x}}_1-{\alpha }_2{\stackrel{\dot{}}{x}}_2={\alpha }_5{x}_1+{\alpha }_6{x}_2+2{\alpha }_9{x}_3+{\alpha }_{12}{x}_1^2+2{\alpha }_{13}{x}_1{x}_3+{\alpha }_{14}{x}_2^2+2{\alpha }_{15}{x}_2{x}_3+c_2.\end{array}$ (2.16)

The coupled nonlinear Equation (2.16) can be solved for $x_3$, resulting in,

$\begin{array}{l}\left({\alpha }_2{\alpha }_{13}+{\alpha }_3{\alpha }_{15}\right)x_3^2+\{\left(2{\alpha }_1{\alpha }_{15}+2{\alpha }_3{\alpha }_{14}\right)x_2+\left(2{\alpha }_1{\alpha }_{13}+2{\alpha }_2{\alpha }_{12}\right)x_1\\ +\left(2{\alpha }_1{\alpha }_9+{\alpha }_2{\alpha }_5+{\alpha }_3{\alpha }_6\right)\}{x}_3+\left({\alpha }_1{\alpha }_5+2{\alpha }_2{\alpha }_7+{\alpha }_3{\alpha }_4\right)x_1\\ +\left({\alpha }_1{\alpha }_6+{\alpha }_2{\alpha }_4+2{\alpha }_3{\alpha }_8\right)x_2+\left({\alpha }_1{\alpha }_{12}+{\alpha }_3{\alpha }_{10}\right)x_1^2\\ +\left(2{\alpha }_3{\alpha }_{11}+2{\alpha }_2{\alpha }_{10}\right)x_1{x}_2+\left({\alpha }_1{\alpha }_{14}+{\alpha }_2{\alpha }_{11}\right)x_2^2+{\alpha }_3{c}_1+{\alpha }_1{c}_2+{\alpha }_2{c}_3=0.\end{array}$ (2.17)

This is another time-independent relation of the 3-variable CNI, but the time-independent relation is obtained from the structure of nonlinear differential equations, not from the Ψ-conservation law. There are two types of time-independent relations, determined by Lagrangian or Noether’s theorem and derived from a particular structure of nonlinear differential equations. The order of nonlinearity, nonlinear coupling terms and coupling strength are usually chosen to agree with global empirical behavior such as exponential decrease or logistic-type convergence which seem reasonable because of dissipative properties in natural world.

Nonlinear problems without classifications of conserving and nonconserving, dissipative phenomena would lead to physically incorrect and ambiguous results, such as the atto-fox problem, interactions among negative values of populations. Some researchers tend to think nonlinear equations have intrinsic drawbacks, which is not true in the nonlinear differential equation with conserving laws. The coupling constants and initial values are strictly restricted to generate physically meaningful solutions, and if conservation laws are neglected, it is easy to generate numerical solutions. Symmetries, conservation laws of nonlinear interactions are important to understand mechanism in microbiological and ecological systems.

3. The Conservation Law and Properties of $\Psi \left(x_1,x_2,\cdots \right)$ -Function

The Ψ-function may have similar physical meanings as Hamiltonian of a system of materials and particles. However, the Hamiltonian has definite meanings as the total energy of a system is a constant in time (the energy conservation law). The energy is convertible to heat, light, electric energy, other kind of work, and has the dimension of force × displacement [42] [43]. The conserved quantity Ψ is constant in time, but it is constructed from variations of populations or density of populations of 2n species, not from the force, kinetic energy and potential energy. The meaning of a conservation law in biological complex systems may be very different from conservation laws in physics, whereas the Ψ-function could be the conserved quantity corresponding to biodiversity.

A nonlinear system demands certain density-independent or density-dependent external perturbations which are not included in Lotka-Volterra and Kolmogorov type nonlinear Equations (6.1) and (6.2). The 2-variable conserving nonlinear differential Equations (2.6)-(2.9) produce generalized nonlinear equations of Lotka-Volterra, Kolmogorov, and Lyapunov function type nonlinear equations, and so, the conserved Ψ-function is related to stability, control and properties of a nonlinear system. The concept of stability is studied in terms of an addition law of Ψ-function [5] [6] [7]. The Ψ-function could be applied to recovery or restoration phenomena caused on a system, which is important for microbiological and ecological systems. The 10-year cycle of Canadian lynx and snowshoe hare is studied as an example of the standard rhythm for recovery or restoration. The standard rhythm would be a consequence of symbiosis [18] [19], living together to survive in severe nature.

3.1. The Classification of Restorations and Disintegrations (Extinctions) in Terms of Ψ-Function

It is essential for nonlinear theoretical and mathematical analyses to have reliable predictions on an ecological system if it thrives or extincts over a given period of time. Let us suppose that a binary-coupled n-variable nonlinear system with $\Psi \left(x_1\mathrm{,}\cdots \mathrm{,}{x}_{2n}\right)$ and another 2-variable interacting system written by $\left(n+1\right)$ -variable $\Psi \left(x_{2n+1}\mathrm{,}{x}_{2n+2}\right)$ interact with each other, resulting in constructing a new, stable and conserving coupled system. The addition property of $\Psi \left(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \right)$ function in the form:

$\Psi \left(x_1\mathrm{,}\cdots \mathrm{,}{x}_{2n}\right)+\Psi \left(x_{2n+1}\mathrm{,}{x}_{2n+2}\right)\to \Psi \left(x_1\mathrm{,}\cdots \mathrm{,}{x}_{2n+2}\right)\mathrm{.}$ (3.1)

may be interpreted as an approximation to the restoration or rehabilitation phenomena known in a large system of biological systems, neural networks, cells of organs, or computer networks replacing a small broken device with a normal device.

The recovery and degeneration in the current study are classified as: 1) Complete Recovery (CR), 2) Partial Recovery (PR), 3) Non-Recovery (NR), Degeneration or Extinction.

1) Complete Recovery (CR)

As shown in Figure 1, when a system maintains the form of $\Psi \left(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \right)$ -function, the interacting system is completely stable. The addition property of Ψ-function means that a defect or negative perturbations will be recovered when a form of binary coupled form is maintained [5]; in other word, a form of CNI to maintain the constant Ψ-function is regained.

2) Partial Recovery (PR)

The partial recovery is not a complete recovery, but a transition to a reasonably stable, approximate state which maintains a conserving stable solution. The Ψ-function changes to a different constant value (see for example, Figure 2).

3) Non-Recovery (NR)

Negative external perturbations or defects of an interacting system eventually lead to the degeneration of an entire interacting system or the extinction, and the Ψ-function becomes negative and is not constant or diverges (see for example, Figure 3 and Figure 4).

There are characteristic combinations of variables, coupling constants to be conservative or dissipative until it reaches an ecological equilibrium. The extinction of species appears as population number zero or negative. The concept of extinction should be categorized as the acute extinction (AEX) and the chronic extinction (CEX). The concept of AEX is well explained with poisons that cause death, which can be specifically determined. On the other hand, the chronic extinction (CEX) is regenerating extinction caused by multi-factored problems, whose causes and effects are difficult to determine. The conserving coupled system and addition law of Ψ-function can help us understand a stable complex system and the concept of recoveries, which is an important result for the conservation of nature and sustainability.

3.2. Recovering, Degeneration and Extinction of Species from External Perturbations

The dynamical interactions of species are affected by climate changes, temporal temperature fluctuations and changes of natural environment by human activities. These factors are mathematically introduced as external perturbations and expressed by piecewise continuous constants, $c_1$ and $c_2$, by using θ-functions such that

$c_i=f_i\left\{\theta \left(t-t_{start}\right)-\theta \left(t-t_{end}\right)\right\},\left(i=1,2\right),$ (3.2)

where $\theta \left(t-t^{\prime }\right)$ represents a step function:

$\theta \left(t-t^{\prime }\right)=(\begin{array}{cc}1,& \left(t\ge t^{\prime }\right),\\ 0,& \left(t<t^{\prime }\right),\end{array}$ (3.3)

and the strength of perturbation $f_i\left(i=1,2\right)$ is positive or negative constants, a free parameter to express external perturbations, adjusted optimally to produce numerical data.

The reaction and the recovery of a nonlinear interacting system from an external perturbation are shown in Figures 2(a)-(c). In Figure 2, an external perturbation starts at $t=700$ (Sp. 1), and the external strength of perturbation, $c_1$, with coefficient $f_1$ equals to −1260.0, and $c_2$ is set to zero. The black arrow is the starting point of perturbation, and the gray arrow exhibits the end of the external perturbation (Ep. 1). The nonlinear coefficients are listed in Table 1 (Condition 1). The solutions $\left(x_1\mathrm{,}{x}_2\right)$ are deformed by the perturbation but, the system does not disintegrate and finds a new stable phase-space and a conserved relation. The perturbation ends at $t=1200$ (Ep. 1), and the system recovers to, more or less, the original state $\left(x_1\mathrm{,}{x}_2\right)$. This is an example of the complete recovery (CR).

(a)(b)(c)

Figure 2. An external negative perturbation and a recovery. (a) The 2-variable CNIs with a negative perturbation on prey $x_1$. The perturbation is introduced from $t=700$ to $t=1200$, which is represented as gray background. (b) The $\left(x_1\mathrm{,}{x}_2\right)$ phase-space transition with the negative perturbation as in (a). Solid line (St. 1) is initial state, and St. 2 is the recovered state after Ep. 1. Dashed line is phase-space during Sp. 1 - Ep. 1. (c) The conservation law Ψ of 2-variable ND with one perturbation. The negative perturbation changed Ψ from ;60000 to ;30000. It recovers after Ep. 1.

It is necessary to distinguish maxima and minima driven by internal (endogenous) and external (exogenous) perturbations. The external perturbations are known in practice as pre-emergent or post-emergent pesticides, herbicides, fungicides and insecticides. It is important to understand the cause of the spray inefficiency or maximum yields of crops for better ecological control of species. The conservation of nature and biodiversity on Earth by protecting species and their habitats from excessive rates of extinction and human activities are one of the essential purposes for theoretical and numerical nonlinear analyses [48].

In Figure 3, the response of a strong negative perturbation to prey after the peak of endogenous maximum is shown. The values of coefficients are listed in Table 1 (Condition 1). The starting point of this perturbation is at $t=800$ and the end point of the perturbation is at $t=950$. The value of negative perturbation to $x_1$, is $f_1=-3175.3879$. The prey, $x_1$, rapidly declines with the negative perturbation, and $\left(x_1\mathrm{,}{x}_2\right)$ converges to zero for $t\gtrsim 1000$. The computer simulations are compatible with the known empirical results, for example, in pest control. A pest control is not so effective if it is performed in the time when harmful insects are in peak and active, because species are energetic enough to find food and mates to thrive. It is effective when a pest control is performed during the time when harmful insects are not so active after endogenous maximum. The acute and strong negative perturbation leads to an acute extinction (AEX) as shown in Figure 3.

Table 1. The list of nonlinear coefficients.

(a)(b)

Figure 3. A critical negative-perturbation leading to extinction. (a) The 2-variable CNIs with a critical negative perturbation on prey, $x_1$. Solutions converge to zero after the perturbation. (b) The conservation law Ψ with a critical negative-perturbation. Ψ converges to zero after the critical negative-perturbation.

In the nonlinear interacting system, positive perturbations which will increase $x_1$ or $x_2$ do not always mean a positive effect on stability of the system. There is a limit to the value of a positive perturbation, because an increase of $x_1$ may lead to a decrease of $x_2$ in the stable system $\left(x_1\mathrm{,}{x}_2\right)$, which indicates that a system has internally allowed maximum and minimum populations, and so, the existence of critical point is important. The behaviors of $\left(x_1\mathrm{,}{x}_2\right)$ at normal and critical values of positive perturbations, $c_1$ for $x_1$ are shown in Figure 4(a), resulting in the extinction of $x_1$. Figure 4(a) and Figure 4(b) show that the system cannot maintain a stable, interacting system when the positive perturbation surpasses a critical value. The unstable solutions branch out at $t\simeq 1100$ when the value of perturbation changes, for example, from $f_1=1160.0$ to $f_1=1599.9$ in the current numerical simulation.

The solutions of $\left(x_1\mathrm{,}{x}_2\right)$ are restricted to be positive integers, $x_i>0$ because a negative solution means the extinction. In Figure 3 and Figure 4(a), solutions $x_1$ and $x_2$ are positive and when they are less than zero, $x_1$ and $x_2$

(a)(b)

Figure 4. Critical perturbations and extinction of species. (a) The 2-variable conserving interaction with a critical perturbation. The species, $x_1$ and $x_2$, converge to zero after a critical perturbation. (b) The conservation law of the 2-variable conserving interaction with a critical perturbation. The Ψ converges to zero after a critical perturbation.

are considered extinct. It should be noted that if a dynamical prey-predator system is active, the rhythm of maxima and minima of interacting species is clearly maintained, which is well known behavior in real prey-predator systems. However, if an external perturbation (exogenous interaction) exceeds a certain critical value of a competitive system, the rhythm of maxima and minima will disappear and then after a time, the system will diverge (disintegrate) as shown in Figure 4. The rhythm of wild-life exhibits stable and dynamical interactions among species, which could be used as a measure for natural environment or wilderness. Wilderness would be defined, for example, by preserves, estates, farms, conservation preserves, ranches, national forests and national parks. When the rhythm of change of species disappears or does not come back, it may indicate that interacting species are in danger of extinction and affected or harmed by some factors usually caused by human activities.

It is possible to save species from extinction by changing the value of perturbations. Figure 5(a) is the result of a positive perturbation to save species $\left(x_1\mathrm{,}{x}_2\right)$ in danger of extinction in Figure 4(a). A positive perturbation is exerted after Sp. 1 - Ep. 1 in Figure 5, starting at $t=1000$ (Sp. 2) and ends at $t=1300$ (Ep. 2); the strengths of $c_1$ and $c_2$ changed to $f_1=200$, $f_2=-1000$. Figure 5 shows that species are in danger of extinction; however, if external perturbations are properly inserted, the system will come back to life again. The executions of appropriate perturbations can save endangered species from extinction, but it is difficult to determine what appropriate positive perturbations for real ecosystems are.

4. The 4-Variable and (n + m)-Variable CNIs

The model of 2n-variable CNIs has exhibited important relations among the conservation law, stability, extinction and restoration phenomena. Although nonlinear interactions have many parameters, values of parameters are not free to adjust because they are primarily constrained by the conservation law and stability. In reality, it is difficult to know even with accumulated big data what kind of stability and the conservation law are linked to initial values, external perturbations and human activities. The concepts of stability, conservation law are fundamental for natural phenomena in physical, biological systems, engineering technologies to realize clean and sustainable energies [49] - [56]. Nonlinearity exists in complicated interacting systems, such as the ecosystems of mammals [10] [11], microbes [12], systems of cells and organs [3]. These systems are highly evolved, interacting systems. We specifically discuss a practical technique for 4-variable and higher variable CNIs in the following discussions.

4.1. The 4-Variable Model for a Conserving Nonlinear Dynamical Interactions

The 4-variable CNI is derived from the 4-variable dynamical Lagrangian [7], and

(a)(b)

Figure 5. The critical behavior, restoration and recovering. (a) The 2-variable conserving interactions with perturbations to avoid converging to zero after a critical perturbation. The species, $x_1$ and $x_2$, come back to life after the second perturbation. (b) The conservation law Ψ with two perturbations in (a). The Ψ-function recovers from the perturbation after Sp. 2 - Ep. 2.

the numbering of nonlinear coefficients, ${\alpha }_i$ in the 4-variable Lagrangian is shown in order to avoid complications and confusions The nonlinear coefficients are numbered sequentially in the beginning as:

$\begin{array}{c}\mathcal{L}={\alpha }_1{\stackrel{\dot{}}{x}}_1{x}_2+{\alpha }_2{x}_1{\stackrel{\dot{}}{x}}_2+{\alpha }_3{\stackrel{\dot{}}{x}}_3{x}_4+{\alpha }_4{x}_3{\stackrel{\dot{}}{x}}_4+{\alpha }_5{x}_1^2+{\alpha }_6{x}_1{x}_2+{\alpha }_7{x}_1{x}_3+{\alpha }_8{x}_1{x}_4\\ +{\alpha }_9{x}_3{x}_1+{\alpha }_{10}{x}_3{x}_2+{\alpha }_{11}{x}_3^2+{\alpha }_{12}{x}_3{x}_4+{\alpha }_{13}{x}_2{x}_1+{\alpha }_{14}{x}_2^2+{\alpha }_{15}{x}_2{x}_3\\ +{\alpha }_{16}{x}_2{x}_4+{\alpha }_{17}{x}_4{x}_1+{\alpha }_{18}{x}_4{x}_2+{\alpha }_{19}{x}_4{x}_3+{\alpha }_{20}{x}_4^2+{\alpha }_{21}{x}_1^2{x}_2+{\alpha }_{22}{x}_1{x}_2^2\\ +{\alpha }_{23}{x}_1{x}_2{x}_3+{\alpha }_{24}{x}_1{x}_2{x}_4+{\alpha }_{25}{x}_3{x}_4{x}_1+{\alpha }_{26}{x}_3{x}_4{x}_2+{\alpha }_{27}{x}_3^2{x}_4+{\alpha }_{28}{x}_3{x}_4^2\\ +c_1{x}_2+c_2{x}_1+c_3{x}_4+c_4{x}_3.\end{array}$ (4.1)

Since the interactions such as $x_1{x}_2=x_2{x}_1$, $x_1{x}_2=x_2{x}_1,\cdots$ are the same in macroscopic systems, the above Lagrangian is expressed as:

$\begin{array}{c}\mathcal{L}={\alpha }_1{\stackrel{\dot{}}{x}}_1{x}_2+{\alpha }_2{x}_1{\stackrel{\dot{}}{x}}_2+{\alpha }_3{\stackrel{\dot{}}{x}}_3{x}_4+{\alpha }_4{x}_3{\stackrel{\dot{}}{x}}_4+\left({\alpha }_6+{\alpha }_{13}\right)x_1{x}_2+\left({\alpha }_7+{\alpha }_9\right)x_1{x}_3\\ +\left({\alpha }_8+{\alpha }_{17}\right)x_1{x}_4+\left({\alpha }_{10}+{\alpha }_{15}\right)x_2{x}_3+\left({\alpha }_{12}+{\alpha }_{19}\right)x_3{x}_4\\ +\left({\alpha }_{16}+{\alpha }_{18}\right)x_2{x}_4+{\alpha }_5{x}_1^2+{\alpha }_{11}{x}_3^2+{\alpha }_{14}{x}_2^2+{\alpha }_{20}{x}_4^2+{\alpha }_{21}{x}_1^2{x}_2+{\alpha }_{22}{x}_1{x}_2^2\\ +{\alpha }_{27}{x}_3^2{x}_4+{\alpha }_{28}{x}_3{x}_4^2+{\alpha }_{23}{x}_1{x}_2{x}_3+{\alpha }_{24}{x}_1{x}_2{x}_4+{\alpha }_{25}{x}_1{x}_3{x}_4+{\alpha }_{26}{x}_2{x}_3{x}_4\\ +c_1{x}_2+c_2{x}_1+c_3{x}_4+c_4{x}_3\mathrm{.}\end{array}$ (4.2)

The above Lagrangian (4.2) is used to derive equations of motion by using Euler-Lagrange method. The 4-variable CNI for complex systems is given by:

$\begin{array}{c}{\stackrel{\dot{}}{x}}_1=\frac{1}{d_{21}}\{\left({\alpha }_6+{\alpha }_{13}\right)x_1+2{\alpha }_{14}{x}_2+\left({\alpha }_{10}+{\alpha }_{15}\right)x_3+\left({\alpha }_{16}+{\alpha }_{18}\right)x_4\\ +{\alpha }_{21}{x}_1^2+2{\alpha }_{22}{x}_1{x}_2+{\alpha }_{23}{x}_1{x}_3+{\alpha }_{24}{x}_1{x}_4+{\alpha }_{26}{x}_3{x}_4\}+\frac{c_1}{d_{21}},\end{array}$ (4.3)

$\begin{array}{c}{\stackrel{\dot{}}{x}}_2=\frac{1}{d_{12}}\{2{\alpha }_5{x}_1+\left({\alpha }_6+{\alpha }_{13}\right)x_2+\left({\alpha }_7+{\alpha }_9\right)x_3+\left({\alpha }_8+{\alpha }_{17}\right)x_4\\ +2{\alpha }_{21}{x}_1{x}_2+{\alpha }_{22}{x}_2^2+{\alpha }_{23}{x}_2{x}_3+{\alpha }_{24}{x}_2{x}_4+{\alpha }_{25}{x}_3{x}_4\}+\frac{c_2}{d_{12}}\mathrm{,}\end{array}$ (4.4)

$\begin{array}{c}{\stackrel{\dot{}}{x}}_3=\frac{1}{d_{43}}\{\left({\alpha }_8+{\alpha }_{17}\right)x_1+\left({\alpha }_{16}+{\alpha }_{18}\right)x_2+\left({\alpha }_{12}+{\alpha }_{19}\right)x_3+2{\alpha }_{20}{x}_4\\ +{\alpha }_{24}{x}_1{x}_2+{\alpha }_{25}{x}_1{x}_3+{\alpha }_{26}{x}_2{x}_3+{\alpha }_{27}{x}_3^2+2{\alpha }_{28}{x}_3{x}_4\}+\frac{c_3}{d_{43}}\mathrm{,}\end{array}$ (4.5)

$\begin{array}{c}{\stackrel{\dot{}}{x}}_4=\frac{1}{d_{34}}\{\left({\alpha }_7+{\alpha }_9\right)x_1+\left({\alpha }_{10}+{\alpha }_{15}\right)x_2+2{\alpha }_{11}{x}_3+\left({\alpha }_{12}+{\alpha }_{19}\right)x_4\\ +{\alpha }_{23}{x}_1{x}_2+{\alpha }_{25}{x}_1{x}_4+{\alpha }_{26}{x}_2{x}_4+2{\alpha }_{27}{x}_3{x}_4+{\alpha }_{28}{x}_4^2\}+\frac{c_4}{d_{34}}\mathrm{,}\end{array}$ (4.6)

where coupling constants are, ${\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}\cdots \mathrm{,}{\alpha }_{28}$, and $d_{12}={\alpha }_1-{\alpha }_2$, $d_{21}={\alpha }_2-{\alpha }_1=-d_{12}$, and $d_{34}={\alpha }_3-{\alpha }_4$, $d_{43}={\alpha }_4-{\alpha }_3=-d_{34}$. The constants, $c_1\mathrm{,}{c}_2,\cdots$, represent strengths of external perturbations or effects, such as temperature changes, reductions of species by huntings or fishings, and so forth. The free parameters are 28 α’s and 4 ( $c_1$ to $c_4$ ) and resulting in 32 nonlinear coefficients to adjust. However, the number of adjustable coefficients can be decreased, since these parameters are restricted to positivity of variables, initial conditions and the conservation law.

Next, it is necessary for numerical analyses to roughly determine coupling constants of linear interactions, assuming coefficients of all nonlinear terms, $x^2$, $x_1{x}_2$, $x_2{x}_2$, $x_3{x}_1\mathrm{,}\cdots$ to be zero, and then, it is better to slowly turn on coefficients of nonlinear terms from 0 to $0\pm \epsilon$ ( $\epsilon$ is a small number) to simulate data and adjust coefficients of linear terms, and numerical simulations are shown in the following section.

The conserved quantity, Ψ, of a 4-variable CNI is constructed according to Noethers’ theorem in Section 2 and expressed as:

$\begin{array}{c}\Psi ={\alpha }_5{x}_1^2+\left({\alpha }_6+{\alpha }_{13}\right)x_1{x}_2+\left({\alpha }_7+{\alpha }_9\right)x_1{x}_3+\left({\alpha }_8+{\alpha }_{17}\right)x_1{x}_4+\left({\alpha }_{10}+{\alpha }_{15}\right)x_2{x}_3\\ +{\alpha }_{11}{x}_3^2+\left({\alpha }_{12}+{\alpha }_{19}\right)x_3{x}_4+{\alpha }_{14}{x}_2^2+\left({\alpha }_{16}+{\alpha }_{18}\right)x_2{x}_4+{\alpha }_{20}{x}_4^2+{\alpha }_{21}{x}_1^2\\ +{\alpha }_{22}{x}_1{x}_2^2+{\alpha }_{23}{x}_1{x}_2{x}_3+{\alpha }_{24}{x}_1{x}_2{x}_4+{\alpha }_{25}{x}_1{x}_3{x}_4+{\alpha }_{26}{x}_2{x}_3{x}_4\\ +{\alpha }_{27}{x}_3^2{x}_4+{\alpha }_{28}{x}_3{x}_4^2+c_1{x}_2+c_2{x}_1+c_3{x}_4+c_4{x}_3\mathrm{.}\end{array}$ (4.7)

The species ( $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4$ ) change in time, but the sum of the combination given by (4.7) is constant, or precisely speaking, piecewise continuous in time as shown in figures from 1 to 5. One should note that the Ψ-function does not have any velocity (time-derivative) terms and the coefficients, ${\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}{\alpha }_3\mathrm{,}{\alpha }_4$, do not exist in the Ψ-function (4.7), though the coefficients exist in the corresponding Lagrangian (4.2). The interpretation and physical meaning of Ψ-function (4.7) would depend on values of coupling constants of nonlinear dynamical interactions, the degree of nonlinearity, time-dependent thermodynamics and irreversiblity [36], which will be investigated in the near future.

The numerical solution to the system of nonlinear differential equations from (4.3) to (4.6) and the value of Ψ-function are shown in Figures 6(a)-(c). The 4 variables ( $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4$ ) are time-changing as in Figure 6(a), but Ψ-function is constant in time as shown in Figure 6(b), which confirms numerical solutions are accurate. The phase spaces of solutions for $\left(x_1\left(t\right)\mathrm{,}{x}_2\left(t\right)\right)$ and $\left(x_3\left(t\right)\mathrm{,}{x}_4\left(t\right)\right)$ are shown in Figure 6(c) and Figure 6(d), respectively. One can see that the phase space of solutions is confined in each volume of phase space in Figure 6(c) and Figure 6(d).

4.2. The Superposition of CNI Solutions

The superposition means the total sum of interacting species, and in practice, it is tacitly assumed that the big data of several interacting species is too big to analyze, and the sum of population data would not be of any use for physical analyses. The superposition of solutions is mathematically exact only in linear differential equations to obtain the general solution. However, we applied it to

(a)(b)(c)(d)

Figure 6. The solutions of 4-variable system and phase-space with respect to time. (a) represents solutions and (b) is the conservation law Ψ with respect to time. The (c) and (d) are phase-space solutions of $x_1$ and $x_2$, $x_1$ and $x_4$, respectively. (a) The 4-variable CNI solutions. (b) The Ψ-function of the 4-variable CNI. (c) The conservation law of the 4-variable CNI solutions. (d) The conservation law of the 4-variable CNI solutions.

the superposition of solutions to the 4-variable conserving nonlinear differential equation, whether it generates any meaningful results or not. Because each solution of CNI shows explicit maximum and minimum, it made us expect that the superposition of solutions:

${\displaystyle \underset{i=1}{\overset{4}{\sum }}}{x}_i=x_1+x_2+x_3+x_4,$ (4.8)

could produce a overall characteristic property of population change, and we found a physically reasonable result on population analyses that should be checked by way of the big data analysis. The superposition of CNI solutions is useful to analyze overall dynamical changes of population data.

The CNI solutions during the time, $0\le t\le 300$ is shown in Figure 7(a), which would impress an arbitrary population change unable to discern any rhythm of ecosystem, but as shown in Figure 7(b), the superposition of $x_1+x_2+x_3+x_4$ seems to show a relative extrema. In order to check if it is a coincidence, the time is extended to $0\le t\le 1000$ to examine changes of extrema in Figure 7(d). One can observe that the value and the time of extrema are different and varying arbitrary, but the extrema of the superposition is explicitly emerged as a long-time behavior in Figure 7(a). A formal rhythm of extrema can be observed by the superposition, $x_1+x_2+x_3+x_4$, and the time-scale change from Figures 7(b)-(d).

(a)(b)(c)(d)

Figure 7. The behavior of 4-variable CNI solutions ( $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4$ ) and the sum of solutions: $x_1+x_2+x_3+x_4$ in different time scales. The extrema of population dynamics could be seen by changing time scales. (a) A system of 4-variable CNI solutions, which would impress an arbitrary change in a time-range: $0\le t\le 300$. (b) The behavior of the sum of 4-variable solutions: $x_1+x_2+x_3+x_4$ (showing about 100 years cycle) in a time range: $0\le t\le 300$. (c) The Ψ-function for 4-variable CNI solutions in a time-range: $0\le t\le 1000$. (d) The behavior of the sum of 4-variable solutions: $x_1+x_2+x_3+x_4$ (showing about 100 years cycle) in a time range: $0\le t\le 1000$.

The emergence of a long-time pattern is a characteristic property in a CNI system. The superposition method may help analyze overall behavior and dynamics of big data, which is too complicated to comprehend by direct observations in short time.

The superposition of conserving nonlinear solutions seems to roughly display the net periodic time, $T_n$, which is about $T_n~100$ in Figure 7(d). The periodic time of component species, $T_i$ ( $i=1,2,3,4$ ), is respectively less than $T_n$, which is not constant like a periodic function, but it looks like fluctuating around the constant mean value of a net periodic time. The behavior persists when time is extended to a longer period. The amplitude and the net periodic time, $T_n$, could be important physical observables to understand ecosystems as conserving nonlinear dynamical systems. The important properties derived from the 4-variable CNI are summarized as follows:

1) The net periodic time $T_n$ defined by the sum of interacting species, ${\displaystyle {\sum }_i}{x}_i$, is observed and longer than the periodic time of component variables $T_i$ : $T_i\le T_n$.

The comparison of $T_n$ and $T_i$ and characteristic extrema could be found in computer analyses of big data of several interacting species in a long period of time. The property should be checked in big-data computer simulations if it exists or not. In computer simulations, one should be careful that the known external perturbations such as human activities, global climate changes, volcanic eruptions and impact of a large asteroid, $\cdots$, should be carefully included. It is very useful to understand population dynamics if long-term cycles of extrema are discovered in marine and land ecosystems.

2) The 4-variable system becomes more stable than the 2-variable system against negative external perturbations.

The property would be a consequence of many coupling terms of the 4-variable CNIs, compared to those of the 2-variable CNIs. In other words, if an interacting system has many interacting terms with other species, the system may not be affected so much by a breakdown of external negative effects, because other interacting terms could alleviate the breakdown and defects on the system, indicating fundamental importance of biodiversity and mutual interactions.

In general, nonlinear systems have different classes of solutions, such as chaos and bifurcation [40] [57] [58]. The property is known that a small change of initial conditions or coupling constants in one state of nonlinear system can result in large differences in a later state. In other words, the underlying patterns and the deterministic laws of chaos and bifurcation are sensitive to initial conditions or coupling constants. The CNI state is similarly sensitive to initial conditions or coupling constants, but it maintains stability of the net system. The results of CNI are different from those discussed in dissipative, nonconserving nonlinear interactions in many literatures which discuss limit cycles and attractors. The conservation law and the stability of CNI system could be a key to understand dynamics of complex systems.

4.3. A method to Construct the 4-Variable and Higher-Variable Solutions by Employing 2-Variable Solutions

The conserving nonlinear n-variable interactions become complicated when the CNI system with $n\ge 3$ is employed. For example, there are 32 nonlinear coefficients for the 4-variable CNI of the type (4.3)-(4.6). If coefficients of the type, ( ${\alpha }_6+{\alpha }_{13}$ ), $\cdots$, and $d_{21}\mathrm{,}{d}_{34}$ can be handled as one parameter, the number of independent parameters can be reduced to 22. Though there are still many free parameters to adjust, admissible solutions would confine values of parameters further when initial starting values and positivity of variables, $x_1\left(t\right)>0$, $x_2\left(t\right)>0$, $x_3\left(t\right)>0$, $x_4\left(t\right)>0$, $\cdots$, are considered in computer simulations.

In order to make the 4-variable CNI system useful for numerical simulations, a strategy to construct a 4-variable solution is explained by employing the 2-variable solutions, which makes finding CNI solutions a little easier than finding the 4-variable solutions directly. Let us denote the 4-variable CNI solution constructed by the sum of a pair of 2-variable CNI solutions as the (2 + 2)-CNI solution.

1) Let us suppose that a 4-variable interacting system is constructed by a pair of 2-variable interaction systems; one system for the interacting species ( $x_1\mathrm{,}{x}_2$ ) and another system for ( $x_3\mathrm{,}{x}_4$ ), and the 2-variable interacting systems begin to interact.

2) The 2-variable interacting species, ( $x_1\mathrm{,}{x}_2$ ), are regarded as solutions for (4.3) and (4.4) by setting coefficients related to $x_3$ and $x_4$ as zero. Specifically, coefficients, $\left({\alpha }_{10}+{\alpha }_{15}\right),\left({\alpha }_{16}+{\alpha }_{18}\right),{\alpha }_{23},{\alpha }_{24},{\alpha }_{25},{\alpha }_{26},\left({\alpha }_7+{\alpha }_9\right),\left({\alpha }_8+{\alpha }_{17}\right)$, are set zero in the beginning. Similarly, the other 2-variable interacting species, ( $x_3\mathrm{,}{x}_4$ ), are also regarded as solutions for (4.5) and (4.6). In other words, it is assumed that a 4-variable interaction ( $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4$ ) is composed of a non-interacting 2-variable ( $x_1\mathrm{,}{x}_2$ ) + 2-variable ( $x_3\mathrm{,}{x}_4$ ) CNIs, in the beginning.

3) The next step is to weakly couple the set of 2-variable CNI by observing coupling terms in the 4-variable solution (4.3)-(4.6). The weak coupling of ( $x_1\mathrm{,}{x}_2$ ) and ( $x_3\mathrm{,}{x}_4$ ) can be performed by adjusting slowly cross-coupling coefficients from 0 to $0\pm \epsilon$ (for instance, $\epsilon ~0.01$ ). This is a trial and error approach; the nonlinear coefficients are numerically restricted within certain small values, or solutions do not exist when values are not compatible with initial conditions and the conservation law, $\Psi \left(x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\right)$.

The reasonable values of nonlinear parameters would be readily searched in numerical simulations by the trial and error. Some numerical values of 2-variable conserving approximations are found in papers [5] [6] [7], and then, reasonable coupling constants can be found by the perturbation method explained above. The strategy can be extended step-by-step to higher variable solutions, for example, a 5-variable CNI solution as a (2 + 3)-variable CNI solution, and a 6-variable CNI solution as a (2 + 2 + 2)-variable CNI solution and so forth.

5. The 10-Year Population Cycles of Lynx and Hare in 200-Year Data and the Purpose of Nonlinear Analysis on Ecology

The competitive systems described by the 2-variable CNI model with external perturbations are applied to population cycles and recovering phenomena of systems from microbes to mammals. We find that nonlinear dynamical systems with a conservation law are stable and generate a characteristic rhythm (cycle) of population density, which we call the standard rhythm of a nonlinear dynamical system. The stability, balance and restoration phenomena of dynamical system are strongly related to the conservation law Ψ, and the standard rhythm of population density would be considered as a manifestation of the survival of strategy for living beings. The famous 10-year cycle of population density of Canadian lynx and snowshoe hare is well known for a prey-predator cycle and numerically analyzed by many researchers as one of nonlinear dynamical interactions. The recorded data prove that Canada lynxes respond to the cyclic rise-and-fall in snowshoe hare populations over the years in Alaska and central Canada. The conservation law would help investigate fundamental relations among environmental problems and managing forests occupied with threatened and endangered species [59].

5.1. The 2-Variable Conserving Approximation Applied to Nonlinear Interactions between Canadian Lynx and Snowshoe Hare

Lynx numbers fluctuate in synchrony with population changes of hares over vast areas in Alaska and central Canada, and the lynx numbers generally delay 1 - 2 years behind the decline in hare numbers. Lynx densities in most central and northern populations are reported to change from 3 to 17-fold during an 8-11 year cycle, tracking the abundance of Snowshoe Hares with 1 - 2 year lag. Field studies have documented from 2 - 45 lynx/100 km² at various times in the cycle and in various habitats. Conservation planning for threatened and endangered species requires reliable information on distributions of species’ habitat [60].

It is reported that there is no evidence to suggest that overall lynx numbers or distribution across Canada have declined significantly over the past two decades, although loss of habitat through increased urbanization and development and forestry is likely affecting lynx populations. Despite the dramatic decrease in harvest through the 1990 cyclic peak, there is evidence that lynx populations in most of the northern range were cycling normally. It is noticeable that field researches report that there is little evidence to conclude that the harvest during the 1980s had a long-term impact on contiguous northern lynx populations [8] [30]. It is necessary to have such a field research to understand ecological environment as well as impact of activity from human societies.

The snowshoe hares were the dominant herbivore, and the changes in their numbers were correlated with those in numbers of arctic ground squirrel, spruce grouse, ptarmigan, lynx, coyote, great horned owl, goshawk, raven and hawk owl. The hare numbers were not correlated with numbers of red-backed vole [59] [60] [61] [62]. The lynx-hare population change is well examined by computer simulations of the 2-variable conserving CNI model with external perturbations as shown in Figure 8, showing that the CNI system keeps stable phase space, Ψ-function in 7.

5.2. The Population Regulation in Canadian Lynx and Snowshoe Hare in 200 Years

It is difficult to identify population regulation mechanism concerning the prey-predator patterns of large mammals because large mammals’ life span is relatively long compared with that of microbes. The prey-predator cycle for wolves

(a)(b)(c)

Figure 8. Several external perturbations and recoveries. (a) The 2-variable ND solutions with three external perturbations. The rhythm of $x_1$ and $x_2$ recovers from several perturbations. The gray background represents periods of external perturbations. (b) The phase-space transitions of $\left(x_1\mathrm{,}{x}_2\right)$. The dashed lines represent solutions, $\left(x_1\mathrm{,}{x}_2\right)$, during perturbations. The solid line represents solutions without perturbations. (c) The conservation law Ψ with three perturbations. It recovers from three perturbations (the gray background in (a)). $\Psi \simeq 60000$ in the St. 1 and St. 2.

and caribou takes some decades to observe; their interacting relation and behaviors were observed with modern technology (GPS-colored animals) [10]. However, the food web configuration between snowshoe hare and Canadian lynx is a well-known prey-predator-type phenomenon, and the ten-year cycle of Canadian lynx was examined from the data of Canadian lynx fur trade returns of the Northern Department of the Hudson’s Bay Company [8]. The Canadian lynx and snowshoe hare have a synchronous ten-year cycle in population numbers [9] [10]. The fundamental mechanism for these cycles is affected by factors, such as nutrient, predation, and environmental factors and human social activities [11]. The prey-predator and interactions among species are classified as the direct interaction, whereas temperature, climate and environmental changes are as the external perturbation.

The lynx and hare population data were provided by the courtesy of the Canadian lynx data of the Environment and Natural Resources, Government of the NWT, Yellowknife, NT, in 2017. We have applied the 2-variable conserving nonlinear equation and obtained compatible results with researchers’ conclusions and remarks on population cycles, but in addition, we obtained different perceptions and implications for population dynamics of lynx and hare. The nonlinear conserving model suggests that species of a system consequently find a strategy or a mechanism to survive for long time periods. In other words, the cycle of population density is a manifestation of the strategy or mechanism to survive, suggested by the stability of phase space solutions of the system.

The data of Canadian lynx and snowshoe hare (1800-2018) is shown in Figure 9(a) and Figure 9(b). The lynx population data is clearly characterized in two periods between (1800-1920) and (1920-2018), which would be related to the modern world history of human societies, while the hare population is in a stable rhythm of the 10-year population cycles. The 2-variable nonlinear conserving model is employed to the lynx-hare data with external perturbations, shown in the Ψ-function in Figure 10. The Lynx data is reproduced in the 2-variable model with external perturbations as in Figure 9(c). The simulation of lynx data around 1800-1920 shows large and abrupt change in the Ψ-function, resulted from large external perturbations. Although the 10-year cycle of lynx population

(a)(b)(c)(d)

Figure 9. The simulation of Canadian lynx and snowshoe hare. (a) The Canadian lynx data of the Environment and Natural Resources 2017. Tundra Ecological Research Station Small Mammal Database. Government of the NWT, Yellowknife, NT. (b) The snowshoe hare data of the Canadian lynx data of the Environment and Natural Resources 2017. Tundra Ecological Research Station Small Mammal Database. Government of the NWT, Yellowknife, NT. (c) The dashed line represents Canadian lynx population simulated by the 2-variable CNI model with external perturbations. (d) The dashed line represents a predicted population of snowshoe hare to that of lynx predicted data of (c).

Figure 10. The change of Ψ-function with external perturbations necessary for simulations of Canadian lynx and snowshoe hare (1800-2018).

has not been disappeared, the value of the Ψ-function has gradually decreased, readily observed by comparing 1821-1921 and 1960-2018 data.

It indicates that human activities affected the lynx population changes, which should be clearly observed, for instance, from the data of Canada lynx fur-trades return of the Hudson’s Bay Company [31]. However, the value of Ψ-function becomes small but changes slowly during 1920-2018, which may indicate an outset of environmental conservation of natural species and resources. The number of lynx is very slowly decreasing, but it is clearly responding to the hare population cycle. In Figure 9(d), it is shown how the hare population can be affected by the lynx population cycle and suggests the hare population is not quantitatively affected by the lynx population cycle, which may express the hares’ vitality and multiple dynamical links to other species.

The feeding and nutrient experiments reported in [11] are considered as external perturbations to the system. As shown in Figure 8, the perturbations cause certain effects on the system, but the system finds a rhythm to maintain the dynamics of species, which is related to the stability of nonlinear dynamical systems. The shapes of phase space $\left(x_1\mathrm{,}{x}_2\right)$ are not so different from the original standard rhythm after several perturbations. Our numerical results agree with conclusions derived from feeding experiments and nutrient-addition experiments. Therefore, the properties of a system with a conservation law have a key to understand the unanswered question why these cycles exist.

The timing of perturbation leads to an essential difference to the feeding experiment of snowshoe hare: ... during the peak of the cycle in 1989 and 1990 (the feeding experiment) had no impact on reproductive output ... however, during the decline phase in 1991 and 1992, the predator exposure plus food treatment caused a dramatic increase in reproductive output ... [11]. The experiment can be studied in our model calculations, showing that perturbations in the increasing and decreasing phase do not cause large effects on standard rhythm. The nonlinear simulations show that negative external perturbations just at the time of a decreasing phase or positive perturbations at an increasing phase induce dramatic effects. For example, in order to increase populations of some species, it is better to let them find sufficient food (positive perturbation) during species’ reproductive, offspring season.

The cycle of standard rhythm for Canadian lynx and snowshoe hare indicates that the stable dynamical system of lynx and hare functions persistently even if environmental nature is ruined by some degree. This is also compatible with the empirical fact that the ten-year cycle in snowshoe hare is resilient to a variety of natural disturbances from forest fires to short-term climatic fluctuations. But it has a limit as we have shown explicitly in Section 3, as recovery, disintegration and extinction phenomena. If a strong negative perturbation is applied persistently for a long period, the system would fall into danger of extinction. The important results of our simulation tell that before a system gets in danger of extinction, the standard rhythm of the system will disappear or tend to become ambiguous. Even if the population of lynx seems to be decreasing in Figure 9(d), the interacting hare species are stable as theoretically predicted in Figure 9(b), which indicates a possible recovery of the lynx species in terms of food and predation. Hence, the lynx population decrease in 1920-2018 should be concluded to be induced by other causes beside food and predation.

We cannot be sure at present what kinds of external perturbations caused the lynx population decrease, which is important to be investigated by local field researchers. A long-term (more than ten years) negative perturbation and a vast environmental change that humans could cause would definitely endanger the standard rhythm of snowshoe hare, lynx and related species. If we carefully observe the standard rhythm of specific ecological systems of species, we could help sustain the dynamics of ecological system and study specific solutions to preserve natural environment and sustainability.

6. Notes on Nonlinear Interactions, Lyapunov Function, Kolmogorov’s Predator-Prey Model, Atto-Fox Problem

The CNI restricts the form of nonlinear interactions and reveals some characteristic properties related to nonlinear dynamics, which helps us understand biological mechanism by way of possible physical analyses. The applicability and physical meanings of nonlinear differential equations with time-dependent coefficients should be studied further for science from microscopic to macroscopic thermomechanical and biological systems.

As physical reasons of nonlinear differential equation have been studied in terms of conserving and nonconserving interactions, nonlinear differential equations with time-dependent coefficients may reveal new mechanism and improve understanding complex phenomena scientifically [36]. While nonlinear differential equations have been used for mathematical analyses and computer simulations to get coupling constants to produce big-data, one should know that nonlinear equations are just a mathematical and experimental tool for classification and systematics (taxonomy) of natural phenomena.

There are unphysical and scientifically invalid problems, misunderstandings on nonlinear problems, such as the atto-fox problem [13] [14] [15] [16] [17]. Specifically, some researchers still insist that non-existing prey (negative number of prey) can be eaten by predators, and predators can increase or decrease by eating or interacting with negative prey, which comes from typical nonconserving nonlinear interactions, even with Lotka-Volterra and other general nonlinear equations. It is important to answer properties of conserving and nonconserving nonlinear dissipative problems in order to study nonlinear dynamics and some typical questions in this section.

6.1. The Conserving and Nonconserving Nonlinear Equations

It is often explained that a model of interacting populations is written in general as,

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=xf\left(x,y\right)\\ \frac{\text{d}y}{\text{d}t}=yg\left(x,y\right)\end{array}$ (6.1)

where $f\left(x,y\right)$ and $g\left(x,y\right)$ are any functions appropriately chosen to study an interacting system, called the autonomous (self-ruling or self-governing) differential equation. The time variable t does not appear explicitly in the functions $f\left(x,y\right)$ and $g\left(x,y\right)$ with the required condition $f\left(0,0\right)=0$ and $g\left(0,0\right)=0$. The meaning of functions depends on problems on which the nonlinear equations are applied, and this model is known as a general differential equation often called Kolmogorov’s predator-prey model [63]. It is not possible to determine whether Equation (6.1) is conserving or nonconserving nonlinear equation, because of the unknown functional form of $f\left(x,y\right)$ and $g\left(x,y\right)$.

Well-known examples in the field of ecology are those of Malthus’ research for a population analysis, A. Lotka (1925) [37] and V. Volterra (1926) [38] for predator-prey differential equations known as Lotka-Volterra (LV) equation. The Kolmogorov’s equation can encompass wide ranges of nonlinear equations with some degrees of freedom to choose functions $f\left(x,y\right)$ and $g\left(x,y\right)$, however, it will not distinguish conserving and nonconserving nonlinear equations at the outset. Hence, it cannot avoid an unphysical problem like the atto-fox problem [13] [14] [15] [16] [17], depending on the form of functions, $f\left(x,y\right)$ and $g\left(x,y\right)$ and initial conditions, $\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)$. If initial conditions and appropriate form of autonomous functions are chosen in terms of Ψ-function in conserving approximations, unphysical problems would be avoided and contribute to useful physical analyses.

6.2. Lyapunov Function Produced by Noether’s Theorem

We showed that the conserved quantity, Ψ-function, can reproduce Lyapunov function of the classical Lotka-Volterra equations in the previous work [5] [6] [7]. The mathematical expressions and physical meanings of $f\left(x,y\right)$, $g\left(x,y\right)$ and stability in differential Equations (6.1) have been investigated by many researchers in terms of limit cycles and atractors of Lyapunov functions [47]. Lyapunov functions are scalar functions denoted as $V_1\left(t\mathrm{,}x\mathrm{,}y\right)$ and $V_2\left(t\mathrm{,}x\mathrm{,}y\right)$, instead of $f\left(x,y\right)$ and $g\left(x,y\right)$ in (6.1), and nonlinear functions are generalized to explicitly include t with the condition $V_1\left(t,0,0\right)=0$ and $V_2\left(t,0,0\right)=0$. There are two types of Lyapunov functions, which are strict Lyapunov (having negative definite time derivative along trajectories) and non-strict Lyapunov functions (negative semi-definite time derivatives along trajectories) [47]. The notion indicates how strongly a nonlinear system approaches its equilibrium state when $t\to \infty$. In other words, a nonlinear interacting system in the Lyapunov type corresponds to a dissipative system or a nonconservative system leading to an equilibrium after a long time.

The systems with Lyapunov function have limit cycles and attractors, which designates energy dissipations of the system. The systems with Ψ-functions are strictly conserving systems corresponding to limit cycles at given time. It is essential to understand that Kolmogorov or Lyapunov functions for systems of differential equations can be derived from Noether’s theorem when a system has an equilibrium state corresponding to certain dynamical conservation laws. The strict and non-strict Lyapunov function could be expressed by nonlinear conserving differential equation corresponding to the global property of Lagrangian approach based on conservation laws. Though Noether’s theorem relates stability and conserved quantity analogous to symmetries and conservation laws of energy and momentum in physics, it is not yet clear what physical meanings the conservation laws have in nonlinear biological and ecosystems. Conservation laws are helpful to understand a system of nonlinear differential equations in view of scientific or physical terms.

6.3. A Note on the Atto-Fox (10⁻¹⁸-Fox) Problem

It should be emphasized that an irrational problem known as atto-Fox (10⁻¹⁸-Fox) problem [16] [17] in a system of nonlinear differential equations should not occur in a realistic ecosystem, or a conserved system of differential equations, because the conservation law, initial conditions and the Ψ-function restrict admissible solutions for nonlinear interacting systems. The nonlinear ordinary differential equations of the 2n-CNI system with a realistic initial condition can have stable solutions, and the solutions consist a stable closed hyper-surface of ( $x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_{2n}$ ), as in Figure 6(c) and Figure 6(d), for example. The problem is restricted by properties of conservation laws, and dissipative non-conserved systems progress to a realistic equilibrium. One can directly observe in the CNI model that the Ψ-function cannot be constant or diverge when a physical solution does not exist with given coupling constants and initial conditions [5] [6] [7].

It is known that typical Lotka-Volterra (LV) equations used in applications produce the atto-Fox (10⁻¹⁸-Fox) problem. Let us prove it by the following LV type nonlinear equation:

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(\alpha -\beta y\right)\\ \frac{\text{d}y}{\text{d}t}=y\left(-\gamma +\delta x\right)\end{array}$ (6.2)

where x and y are the population number of prey and predator respectively; $\text{d}x/\text{d}t$ and $\text{d}y/\text{d}t$ represent variations of populations with respect to time. The constants, $\alpha \mathrm{,}\beta \mathrm{,}\delta \mathrm{,}\gamma$, are arbitrary nonlinear coefficients whose interpretations depend on problems at hand, and different nonlinear terms are readily included in the equation. It is not generally possible to obtain an analytical solution to the nonlinear differential equation, and a solution to (6.2) must be obtained by numerical computations with given coefficients as shown in Figure 11 and Figure 12.

Living animals interact with other animals through reasonably large integer numbers as we know in a common sense. The lower left of Figure 12 shows that the small numbers such as $\left(x\mathrm{,}y\right)~\left(\mathrm{1,90}\right)$ at certain time t can interact with each other and continue to live, which is absurd and meaningless. More specifically, the LV Equation (6.2) is invariant with the scale change:

Figure 11. The values of nonlinear coefficients are, $\alpha =0.35$, $\beta =0.0035$, $\gamma =0.25$, $\delta =0.0020$ with initial values, $x\left(0\right)=200$ and $y\left(0\right)=500$.

Figure 12. The x-y phase space solution. The LV equation develops from (0, 0) to finite numbers.

$X\left(t\right)=ax\left(t\right),Y\left(t\right)=by\left(t\right),$ (6.3)

where a and b are arbitrary numbers. The equation for $X\left(t\right)$ and $Y\left(t\right)$ becomes:

$\begin{array}{l}\frac{\text{d}X}{\text{d}t}=X\left(\alpha -{\beta }^{\prime }Y\right)\\ \frac{\text{d}Y}{\text{d}t}=Y\left(-\gamma +{\delta }^{\prime }X\right)\end{array}$ (6.4)

with ${\beta }^{\prime }=\beta /b$ and ${\delta }^{\prime }=\delta /a$. Hence, the solution $\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)$ in Figure 11 and Figure 12 can be equivalently constructed by the solution $\left(X\left(t\right)\mathrm{,}Y\left(t\right)\right)$.

If we choose, for example, $a={10}^{18}$ and $b=1$, the solution can be expressed as $x\left(t\right)={10}^{-18}X\left(t\right)$ and $y\left(t\right)=Y\left(t\right)$. The solution $\left(x\left(t\right),y\left(t\right)\right)=\left({10}^{-18}X\left(t\right),Y\left(t\right)\right)$ is obtained, and the system develops from (0, 0) to finite numbers of species (X, Y) with any initial conditions, which proves the atto-fox (10⁻¹⁸-fox) problem. The LV equation of the type (1.2) has a serious problem when it is applied to real biological and ecological problems. It should be used as an qualitative analysis at most. However, it is a misunderstanding that the nonlinear differential equation of interacting species is always constrained by the atto-fox (0⁻¹⁸-fox) problem.

The LV equations may be summarized as: 1) the LV Equation (6.2) is too simple to apply to ecological and interactive systems. 2) density-independent external perturbations should be included in numerical simulations in order to include changes of climate, temperature, seasons and landscape and so forth. The external perturbations are independent of $x\left(t\right)$ and $y\left(t\right)$. 3) the nonlinearity, initial conditions and nonlinear coefficients restrict values of coupling constants and avoid the atto-fox problem. In addition, nonlinear differential equations with or without time-dependent coefficients should be categorized in different classes of differential equations, which could provide with new and essential information of interacting systems [36].

6.4. Answers to Some Questions on Nonlinear Equations

Some typical questions concerning difficulties on nonlinear differential equations are answered in order to elucidate applicability of the 2n-variable CNI equations.

● (Question 1) There may be many periodic solutions in the n-variable nonlinear model that depend on the values of nonlinear coefficients, ${\alpha }_1\mathrm{,}\cdots$, and initial starting values, $x_1\left(0\right)\mathrm{,}\cdots \mathrm{,}{x}_n\left(0\right)$.

● (Answer) It is indeterminable from the beginning to know how many independent solutions a nonlinear differential equation can produce, or it could have none. This is restated such that solutions ( $x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_n$ ) would be transformed in the phase space like Figure 12, resulting in the atto-fox problem or a physically meaningless problem, which may be correct for simple LV type nonlinear equations. However, one should be careful with the n-variable CNI model. The nonlinear coefficients and initial values in the n-variable model are restricted with the conservation law Ψ and positivity $x_i\ge 0$. The adjustable nonlinear parameters are not entirely free, restricting physically meaningless atto-fox problems. The consistent values of initial starting values and nonlinear coefficients should be determined from numerical simulations by including external perturbations.

● (Question 2) Selecting initial conditions of scarce prey and plentiful predators, unphysical solutions to nonlinear equations could occur. For example, $x_1$ may decrease even if it is zero when $x_2$ is positive and increasing, and $x_2$ may decrease or increase even if $x_1$ is negative. This means that non-existing prey can be eaten by predators, and predators can increase or decrease by eating or interacting with negative prey.

● (Answer) This is a typical conjecture derived from LV type nonlinear equations and the comment is not correct at all regarding the conserving nonlinear model. Because nonlinear coefficients are not entirely free parameters, it is not possible to find solutions by selecting any nonlinear coefficients and initial conditions of scarce prey and plentiful predators, and atto-fox type problems are prohibited. A specific numerical example is shown in the paper [6] [7]. The solutions to n-variable CNI equations are different from those of conventional, nonconserving nonlinear equations.

It is important to stress that the atto-fox problem found in a simple Lotka-Volterra equation is not intrinsic to conserving nonlinear differential equations. The n-variable CNI equations with external perturbations are useful to simulate real data numerically.

7. Conclusions on the Survival of the Fittest and Symbiosis from Conserving Nonlinear Interactions

We examined characteristic properties of several ecological systems based on CNIs which include generalized Lotka-Volterra type prey-predator, competitive interactions. The practical construction of solutions from 2-variable to n-variable CNIs is discussed in order to apply the model to more realistic biological phenomena and responses to external environmental changes. It is difficult to adjust values of coefficients for nonlinear differential equations to satisfy given initial conditions and real population ecological data. It is emphasized in the current CNI model that nonlinear differential equations are more than a convenient computer analysis of big data. When an appropriate set of nonlinear parameters in the conserving nonlinear model is found, it can be used to examine characteristic properties, such as stability and standard rhythm internally related to dynamics of a corresponding nonlinear system.

The important factors (nutrient, predation and social interactions) are needed for all species to survive in nature, but they are easily altered with environmental and natural conditions. An abnormal increase in population numbers of a species would endanger the survival of a species itself as well as other species, which could be observed by examining the population rhythm. Even if a mechanism of population increase is unknown, the CNI method would help study a cause in terms of the conservation law of the ecosystem, stability and recovering strength to external perturbations. It seems that as a predator needs a prey for its food, a prey needs a predator for the conservation of its own species. Nature simultaneously nurses and demolishes life because of the finite resources of food, energy and environment of the Earth. The conservation law and rhythm of species could be considered to have been constructed by species to survive in natural conditions for a long time with evolutionary strategies. Hence, the cycle (rhythm) of species would be interpreted as manifestation of the survival of the fittest to the balance of a biological system.

Conservation of forest-associated animals requires an understanding of habitat quality, along with the development of environmental spatio-temporal management strategies that are consistent with high-quality habitat [59]. The natural environment and habitat quality should be conserved for species, and so such conditions would manifest themselves as a dynamical rhythm of interaction and be related to conservation laws, whereas temporal and spatio-temporal conditions should be considered as external perturbations. The current conserving approximations with external perturbations help us understand cause and effect for complicated natural phenomena of ecosystems.

Nonlinear differential equations have bifurcation solutions depending on coupling strengths which cause a sudden qualitative and quantitative change into solutions, even when a small smooth change is made to parameter values (coefficients of given nonlinear differential equations) of nonlinear differential equations. In addition, nonlinear differential equations with time-dependent coefficients generate a different class of independent solutions compared to nonlinear differential equations with constant coefficients [36]. Nonlinear differential equations seem to have yet unexplored properties, which should be investigated further. The nonlinear differential equations with time-dependent coefficients would be related to the reaction-diffusion interactions of density-dependent growth and dispersal of viruses, density-dependent diffusion coefficient of cooperative phenomena and the evolutionary population dynamics [64]. The applications of CNI model to other dynamical fields will be investigated in the future.

We conclude that stability and conservation law are constructed by species in mutual dependence or cooperation to survive for long-time periods in severe nature. The standard rhythm should be regarded as the result of strategy for species to survive in nature. Therefore, the conserving nonlinear differential equations are suitable to reveal properties of complicated interacting ecological systems.

Whatever roles plants and animals have to play, species that can fit and balance other creatures can survive in nature. A strong predator cannot even survive if it ignores the law of the standard rhythm and conservation law of ecosystems constructed by other members and the environment. We hope that this study will help understand dynamical balance and activities for living animals and importance of natural environment for life.

Acknowledgements

The authorized use of data is acknowledged by: Environment and Natural Resources, 2017; Tundra Ecological Research Station Small Mammal Database; NWT Wildlife Management Information System; Government of the NWT, Yellowknife, NT (Data-use-agreement No: 8515-05 DR367).

The author, H. Uechi, acknowledges that the conserving nonlinear analyses and thermomechanical dynamics employed here are partly supported by Kansai Research Foundation for Technology Promotion (KRF), 2021, Osaka, Japan.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	CITES (Accessed in 2021) The Convention on International Trade in Endangered Species of Wild Fauna and Flora. https://cites.org/eng
[2]	Schindler, D.E., Carter, J.L., Francis, T.B., Lisi, P.J., Askey, P.J. and Sebastian, D.C. (2012) Mysis in the Okanagan Lake Food Web: A Time-Series Analysis of Interaction Strengths in an Invaded Plankton Community. Aquatic Ecology, 46, 215-227. https://doi.org/10.1007/s10452-012-9393-0
[3]	Chyu, J., Fonarow, G.C., Tseng, C.H. and Horwich, T.B. (2014) Four-Variable Risk Model in Men and Women with Heart Failure. Circulation: Heart Failure, 7, 88-95. https://doi.org/10.1161/CIRCHEARTFAILURE.113.000404
[4]	Nair, S., Mishra, V., Hayden, K., Lisboa, P.J.G., Pandya, B., Vinjamuri, S., Hardy, K.J. and Wilding, J.P.H. (2011) The Four-Variable Modification of Diet in Renal Disease Formula Underestimates Glomerular Filtration Rate in Obese Type 2 Diabetic Individuals with Chronic Kidney Disease. Diabetologia, 54, 1304-1307. https://doi.org/10.1007/s00125-011-2085-9
[5]	Uechi, L. and Akutsu, T. (2012) Conservation Laws and Symmetries in Competitive Systems. Progress of Theoretical Physics Supplement, 194, 210-222. https://doi.org/10.1143/PTPS.194.210
[6]	Uechi, L. and Akutsu, T. (2013) Stability and Restoration Phenomena in Competitive Systems. Progress of Theoretical and Experimental Physics, 2013, Article No. 103J01. https://doi.org/10.1093/ptep/ptt077
[7]	Uechi, L. and Uechi, H. (2016) Noether’s Conservation Laws and Stability in Nonlinear Conservative Interactions. Open Access Library Journal, 3, Article ID: e2592. https://doi.org/10.4236/oalib.1102592
[8]	Elton, C. and Nicholson, M. (1942) The Ten-Year Cycle in Numbers of the Lynx in Canada. The Journal of Animal Ecology, 11, 215-244. https://doi.org/10.2307/1358
[9]	Stenseth, N.C., Falck, W., Bjørnstad, O.N. and Krebs, C.J. (1997) Population Regulation in Snowshoe Hare and Canadian Lynx: Asymmetric Food Web Configurations between Hare and Lynx. Proceedings of the National Academy of Sciences of the United States of America, 94, 5147-5152. https://doi.org/10.1073/pnas.94.10.5147
[10]	Stenseth, N.C., Falck, W., et al. (1998) From Patterns to Processes: Phase and Density Dependencies in the Canadian Lynx Cycle. Proceedings of the National Academy of Sciences of the United States of America, 95, 15430-15435. https://doi.org/10.1073/pnas.95.26.15430
[11]	Krebs, C.J., Boonstra, R., Boutin, S. and Sinclair, A.R.E. (2001) What Drives the 10-Year Cycle of Snowshoe Hares? BioScience, 51, 25-35. https://doi.org/10.1641/0006-3568(2001)051[0025:WDTYCO]2.0.CO;2
[12]	Gazi, N.H. (2012) Dynamics of a Marine Plankton System: Diffusive Instability and Pattern Formation. Applied Mathematics and Computation, 218, 8895-8905. https://doi.org/10.1016/j.amc.2012.02.048
[13]	Lobry, C. and Sari, T. (2015) Migrations in the Rosenzweig-MacArthur Model and the Atto-Fox Problem. ARIMA Journal, 20, 95-125. https://doi.org/10.46298/arima.1990
[14]	Mainhardt, H. (1982) Models of Biological Pattern Formation. Academic Press, Cambridge.
[15]	Rubin, A. and Riznichenko, G. (2014) Mathematical Biophysics. Springer, Boston. https://doi.org/10.1007/978-1-4614-8702-9
[16]	Mollison, D. (1991) Dependence of Epidemic and Population Velocities on Basic Parameters. Mathematical Biosciences, 107, 255-287. https://doi.org/10.1016/0025-5564(91)90009-8
[17]	Kloeden, P.E. and Pötzsche, C. (2010) Dynamics of Modified Predator-Prey Models. International Journal of Bifurcation and Chaos, 20, 2657-2669. https://doi.org/10.1142/S0218127410027271
[18]	Margulis, L. (1998) Symbiotic Planet. Basic Books, New York.
[19]	Margulis, L. and Sagan, D. (1995) What Is Life. California University Press, Berkeley.
[20]	Wade, M.J. (2007) The Co-Evolutionary Genetics of Ecological Communities. Nature Reviews Genetics, 8, 185-195. https://doi.org/10.1038/nrg2031
[21]	Allesina, S., Alonso, D. and Pascual, M. (2008) A General Model for Food Web Structure. Science, 320, 6858-661. https://doi.org/10.1126/science.1156269
[22]	Lafferty, K.D., Dobson, A.P. and Kuris, A.M. (2006) Parasites Dominate Food Web Links. Proceedings of the National Academy of Sciences of the United States of America, 103, 11211-11216. https://doi.org/10.1073/pnas.0604755103
[23]	Zhang, L., Hou, D., Chen, X., Li, D., Zhu, L., Zhang, Y., et al. (2012) Exogenous Plant MIR168a Specifically Targets Mammalian LDLRAP1: Evidence of Cross-Kingdom Regulation by MicroRNA. Cell Research, 22, 107-126. https://doi.org/10.1038/cr.2011.158
[24]	Takefumi, M., Yamaguchi, S., Yamamoto, I., Yamaoka, R. and Akino, T. (2014) “Double-Trick” Visual and Chemical Mimicry by the Juvenile Orchid Mantis Hymenopus coronatus Used in Predation of the Oriental Honeybee Apis cerana. Zoological Science, 31, 795-801. https://doi.org/10.2108/zs140126
[25]	Rees, A. (2006) Generically Modified Food. Pluto Press, London.
[26]	Smith, J.M. (2018) Healing from GMOs and Roundup. Smith & Ostara.
[27]	Coppex, F., Droz, M. and Lipowski, A. (2004) Extinction Dynamics of Lotka-Volterra Ecosystems on Evolving Networks. Physical Review E, 69, Article ID: 061901. https://doi.org/10.1103/PhysRevE.69.061901
[28]	Parker, M. and Kamenev, A. (2009) Extinction in the Lotka-Volterra Model. Physical Review E, 80, Article ID: 021129. https://doi.org/10.1103/PhysRevE.80.021129
[29]	Hunter, L. (2015) Wild Cats of the World. Bloomsbury Publishing, London.
[30]	Brand, C.J. and Keith, L.B. (1979) Lynx Demography during a Snowshoe Hare Decline in Alberta. The Journal of Wildlife Management, 43, 827-849. https://doi.org/10.2307/3808267
[31]	Hudson’s Bay Company (HBC) (2019). http://www3.hbc.com/
[32]	MacLulich, D.A. (1937) Fluctuations in the Numbers of Varying Hare (Lepus Americanus). No. 43. The University of Toronto Press, Toronto. https://doi.org/10.3138/9781487583064
[33]	Government of Northwest Territories (2018) Database Release Agreement.
[34]	Ikeda, M. and Siljak, D.D. (1980) Lotka-Volterra Equations: Decomposition, Stability, and Structure. Journal of Mathematical Biology, 9, 65-83. https://doi.org/10.1007/BF00276036
[35]	Matsuda, H., Ogita, N., Sasaki, A. and Satō, K. (1992) Statistical Mechanics of Population. Progress of Theoretical Physics, 88, 1035-1048. https://doi.org/10.1143/ptp/88.6.1035
[36]	Uechi, H, Uechi, L. and Uechi, S.T. (2021) Thermodynamic Consistency and Thermomechanical Dynamics (TMD) for Nonequilibrium Irreversible Mechanism of Heat Engines. Journal of Applied Mathematics and Physics, 9, 1364-1390. https://doi.org/10.4236/jamp.2021.96093
[37]	Lotka, A. (1925) Elements of Physical Biology. Williams and Wilkins, Baltimore.
[38]	Volterra, V. (1926) Variazioni e fluttuazioni del numero di individui in specie animali conviventi. Memoria della Reale Accademia Nazionale dei Lincei, 2, 31-113.
[39]	Liu, B., Zhang, Y. and Chen, L. (2005) The Dynamical Behaviors of a Lotka-Volterra Predator-Prey Model Concerning Integrated Pest Management. Nonlinear Analysis Real world applications, 6, 227-243. https://doi.org/10.1016/j.nonrwa.2004.08.001
[40]	Liu, B., Teng, Z. and Liu, W. (2007) Dynamic Behaviors of the Periodic Lotka-Volterra Competing System with Impulsive Perturbations. Chaos, Solitons and Fractals, 31, 356-370. https://doi.org/10.1016/j.chaos.2005.09.059
[41]	Brauer, F. and Castillo-Chavez, C. (2012) Mathematical Models in Population Biology and Epidemiology. Springer, New York. https://doi.org/10.1007/978-1-4614-1686-9
[42]	Logan, J.D. (1977) Invariant Variational Principles. Academic Press, New York.
[43]	Goldstein, H., Poole, C.P. and Safko, J.L. (2002) Classical Mechanics. 3rd Edition, Addison-Wesley, New York. https://doi.org/10.1119/1.1484149
[44]	Sieniutycz, S. (1994) Conservation Laws in Variational Thermo-Hydrodynamics. Springer, Kluwer Academic Pub., Norwell. https://doi.org/10.1007/978-94-011-1084-6
[45]	Sieniutycz, S. and Farkas, H. (2005) Variational and Extremum Principles in Macroscopic Systems. Elsevier, Amsterdam.
[46]	Vujanovic, B.D. and Jones, S.E. (1989) Variational Methods in Nonconservative Phenomena. Academic Press, Cambridge.
[47]	Mazenc, F. and Malisoff, M. (2010) Strict Lyapunov Function Constructions under LaSalle Conditions with an Application to Lotka-Volterra Systems. IEEE Transactions on Automatic Control, 55, 841-854. https://doi.org/10.1109/TAC.2010.2041995
[48]	Sahney, S. and Benton, M.J. (2008) Recovery from the Most Profound Mass Extinction of All Time. Proceedings of the Royal Society B Biological Sciences, 275, 759-765. https://doi.org/10.1098/rspb.2007.1370
[49]	Ramos, A.F., Innocentini, G.C.P., Forger, F.M. and Hornos, J.E.M. (2010) Symmetry in Biology: From Genetic Code to Stochastic Gene Regulation. IET systems biology, 4, 311-329. https://doi.org/10.1049/iet-syb.2010.0058
[50]	Uechi, S.T. and Uechi, H. (2018) The Profiling of International Roughness Index (IRI) Based on Lagrangian Method. World Journal of Engineering and Technology, 6, 885-902. https://doi.org/10.4236/wjet.2018.64059
[51]	Uechi, S.T. and Uechi, H. (2019) A Mechanical Vibration-Induced Electric Energy Generation (MVEG) and Applications to Ride Quality of Vehicles and International Roughness Index (IRI). Studies in Engineering and Technology, 6, 59-69. https://doi.org/10.11114/set.v6i1.4301
[52]	Førland, K.S., Førland, T. and Ratkje, S.K. (1988) Irreversible Thermodynamics; Theory and Applications. John Wiley & Sons.
[53]	Ives, A.R. and Carpenter, S.R. (2007) Stability and Diversity of Ecosystems. Science, 317, 58-62. https://doi.org/10.1126/science.1133258
[54]	Allesina, S. and Tang, S. (2012) Stability Criteria for Complex Ecosystems. Nature, 483, 205-208. https://doi.org/10.1038/nature10832
[55]	Thébault, E. and Fontaine, C. (2010) Stability of Ecological Communities and the Architecture of Mutualistic and Trophic Networks. Science, 329, 853-856. https://doi.org/10.1126/science.1188321
[56]	Daafouz, J., Riedinger, P. and Iung, C. (2002) Stability Analysis and Control Synthesis for Switched Systems: A Switched Lyapunov Function Approach. IEEE Transactions on Automatic Control, 47, 1883-1887. https://doi.org/10.1109/TAC.2002.804474
[57]	Strogatz, S.H. (1994) Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Westview Press, Boulder.
[58]	Smith, H.L. and Waltman, P. (1995) The Theory of the Chemostat. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511530043
[59]	Holbrook, J.D., Squires, J.R., Bollenbacher, B., Graham, R., Olson, L.E., Hanvey, G., et al. (2019) Management of Forests and Forest Carnivores: Relating Landscape Mosaics to Habitat Quality of Canada Lynx at Their Range Periphery. Forest Ecology and Management, 437, 411-425. https://doi.org/10.1016/j.foreco.2019.01.011
[60]	Simons-Legaard, E.M., Harrison, D.J., Krohn, W.B. and Vashon, J.H. (2013) Canada Lynx Occurrence and Forest Management in the Acadian Forest. The Journal of Wildlife Management, 77, 567-578. https://doi.org/10.1002/jwmg.508
[61]	Boutin, S., Krebs, C.J., Boonstra, R., Dale, M.R.T., Hannon, S.J., Martin, K., et al. (1995) Population Changes of the Vertebrate Community during a Snowshoe Hare Cycle in Canada’s Boreal Forest, Oikos, 74, 69-80. https://doi.org/10.2307/3545676
[62]	Poole, K.G. (2003) A Review of the Canada Lynx, Lynx Canadensis, in Canada. Canadian Field-Naturalist, 117, 360-376. https://doi.org/10.22621/cfn.v117i3.738
[63]	Freedman, H.I. (1980) Deterministic Mathematical Models in Population Ecology. Marcel Dekker Inc., New York.
[64]	Hunter, M., Krishnan, N., Liu, T., Möbius, W. and Fusco, D. (2021) Virus-Host Interactions Shape Viral Dispersal Giving Rise to Distinct Classes of Traveling Waves in Spatial Expansions. Physical Review X, 11, Article ID: 021066. https://doi.org/10.1103/PhysRevX.11.021066