Environmental sustainability is an important concept that aims to preserve the Earth’s natural and environmental resources to ensure the quality of life today and in the future, as it depends on achieving a balance between human needs, preserving the environment, and eliminating negative impacts on the ecosystem. Diverse activities and tasks across disciplines can be represented in a mathematical model that encompasses the diversity of studies across disciplines in the environmental sustainability system. This principle (sustainability) has become widespread and can be applied to all aspects of life on Earth, at different times, and in different places, such as forests, green areas, or wetlands. Healthy forests are among the most important examples of healthy and vital systems for sustainability. ^{1– 3} To provide healthy and vital forests to ensure environmental sustainability, the balance in the ecosystem must be studied through a mathematical model that examines organisms and their relationship with each other, especially the relationship of coexistence between living organisms, as well as reducing the severity of predation by predators on organisms that are considered easy prey. This was the basis of our study. Many researchers in the field of ecosystems and the study of cases of environmental stability and equilibrium have expressed the system using a mathematical model that studies the relationship between living organisms. These systems have become widespread in all sciences and disciplines, and include diverse sciences in the biological sciences, infectious diseases, physiology, environmental sciences, immunology, economics, chemical reactions, genetics, and pharmacodynamics. ^{4, 5} In addition to these systems, there are universally agreed upon principles that the relationship between living organisms (humans, animals, and plants) constitutes a key link in the chain of universal life. This relationship is represented by an interaction represented by a mathematical model, which is an important application of applied mathematics in the animal world; weak animals (prey) can become a major food source for other animals (predators). In this case, the role of mathematicians is to choose the model that achieves the principles of safety and environmental sustainability. ^{6, 7} This interaction is called a predator-prey relationship. ^{8– 11} Two different types of animals and organisms can also interact in a mutually beneficial relationship, known as symbiosis. ¹² An example is the nectar of flowering plants and bees that feed on it. Nectar is the food of bees, whereas bees pollinate flowers. ^{13, 14} The other interaction, which is one of the most important, can be represented by a single food source shared between two or more species of organisms competing for it or organisms competing for survival. ¹⁵ This type can be described by a set of mostly nonlinear equations (perhaps algebraic or differential equations). Therefore, the science of dynamic behavior in biological systems has become dominant in modeling science over the last few years. One of the important issues that must be highlighted when studying the predator-prey model is that logistic growth reaches the highest carrying capacity in the prey community if the predator is absent, and gives complete equilibrium and an ascending growth curve for the prey in the absence of predators. For this reason, most modelers and biologists have devoted their attention and time to this species, studying it closely and among them. ^{16– 21} Hence, biologists have begun to use these models to help them understand and study the biological aspects of organisms under their natural conditions. ^{22– 25} Recently, many researchers have been interested in a specific issue, which is the fear of the prey and the effect of its fear of the predator, and discovered through mathematical results that this effect is very negative. ^{26, 27} In addition, a study was conducted in ^{28, 29} on the impact of direct killing and hunting on birds and invertebrates during their breeding seasons to clarify how predators are affected. A predator-prey model with non-local competition can produce complex dynamics, such as spatiotemporal patterns and spatially heterogeneous periodic solutions. See ³⁰ Ghanem, H. and Majeed studied the dynamics of a predator-prey model including two types of epidemics with harvest. The basic structure of a mathematical model with a biological meaning and character that can be studied accurately and carefully, based on the most important members of the population, is the prey and predator, which is described as the density of the prey consumed per unit of time by the predators. ^{31– 33} Latest I have come up with Ghanem, H. and Khaleefah, K in, ³⁴ a special study on the prey and predator patterns of Lake Victoria tilapia, snook, and coptodon, and recommendations for their conservation from overfishing. The study in this research is based on proposing a mathematical environmental (biological) model for living organisms, achieving a large part of the principles of environmental sustainability by creating the necessary conditions for coexistence and exchange of benefits between organisms that were considered prey for the purpose of controlling their extinction, which causes the death of a large part of the environment and an effective resource in the environmental balance, which consists of a system of differential equations. The local and global stability were analyzed and discussed, and the results were analyzed using numerical simulations related to the conditions of equilibrium and environmental sustainability. The results were positive for the purpose of providing a sustainable ecosystem.

Consider the proposed cooperative ecological mathematical model of three organisms, which involves a predator feeding on two prey species in a symbiotic relationship and exchanging food and resources for survival. Lotka-Volterra considered the response function (Holling function of the first kind) as a response dynamic process. Therefore, the population density hypothesis represented by the symbols X(T), Y(T), and Z(T) is that the first, second, and predator feed on them at any time (T) respectively. The growth density of the two prey groups represented a positive natural growth rate for both species (logistic growth). The relationship between organisms was described according to their effect on each other, where coexistence was considered a close relationship between the first and second prey, and the third organism was considered a predator of the first two species, in addition to the intense competition within each community, which may have been for food or for females during the mating season. Therefore, we describe how population density changes over time for the proposed species by partitioning the system into a set of three independent ordinary differential equations. Based on this model, see the following system: dX dT = a 1 ( 1 − X ) X + bXY − λ 1 XZ − d 1 X dY dT = a 2 ( 1 − Y ) Y + cXY − λ 2 YZ − d 2 Y dZ dT = e 1 λ 1 XZ + e 2 λ 2 YZ − ( d 3 + d 4 Z ) Z (1)

With initial values of variables (X(0) ⩾ 0, Y(0) ⩾ 0, Z(0) ⩾ 0) .

The Table 1 describes the basic parameters of the model (1).

We study the initial conditions of model (1) in space ( R + 3 ) and show that the functions of model (1) are continuous and differentiable on ( R + 3 ) , thus, the Lipschitz condition is satisfied, which is proof of the existence and uniqueness of model (1). In addition, all solutions of model (1) with its nonnegative initial condition are completely determined, as shown in the following theorem: Theorem 1:

The solutions to Model (1) are uniformly bounded.

The study in this part is concerned with searching for the existence of possible and expected (fixed) equilibrium points for Model (1), where it was noted that no more than seven points are the possible equilibrium points. 1-

The point that always exists is the simple equilibrium point, E 0 = (0, 0, 0).

This point is present if we adopt conditions ( 2 ) and ( 3 ) in addition to the following condition: ( x ~ > 1 ) and ( c b > a 1 a 2 ) (4) 5-

The fifth point of balance E ̇ 4 = ( 0 , y ̇ , z ̇ ) where y ̇ = a 2 d 4 − d 2 d 4 + λ 2 d 3 a 2 d 4 + e 2 ( λ 2 ) 2 and z ̇ = e 2 λ 2 y ̇ − d 3 d 4

E ̇ 4 is found based on condition ( 3) add the following condition: ( e 2 λ 2 y ̇ > d 3 ) (5) 6-

Sixth balance point E ̿ 5 = ( x ̿ , 0 , z ̿ ) where:

E ̿ 5 is found based on condition ( 2) add the following condition: ( e 1 λ 1 x ̿ > d 3 ) (6) 7-

The seventh point is called the positive balance point or the coexistence point E 6 ∗ = ( x ∗ , y ∗ , z ∗ )

Where: x ∗ = ( a 1 − d 1 − λ 2 z ∗ ) ( c b − a 1 a 2 ) − a 1 b ( λ 2 z ∗ + d 2 ) + a 1 b ( a 2 + c ) − cb ( λ 1 z ∗ + d 1 ) ( c b − a 1 a 2 ) y ∗ = a 1 a 2 + a 1 c − c λ 1 z ∗ − a 1 λ 2 z ∗ − c d 1 − a 1 d 2 ( c b − a 1 a 2 ) z ∗ = S 1 ∗ S 2 ∗ + S 3 ∗ S 1 ∗ = − ( c e 2 λ 1 λ 2 + a 1 e 2 ( λ 2 ) 2 − ( c b − a 1 a 2 ) ( e 1 λ 1 λ 2 − d 4 ) + c b e 2 ( λ 1 ) 2 + b a 1 e 1 λ 1 λ 2 ) S 2 ∗ = a 1 a 2 e 2 λ 2 + c a 1 e 2 λ 2 − c e 2 λ 2 d 1 − a 1 e 2 λ 2 d 2 + ( cb − a 1 a 2 ) ( e 1 λ 1 d 1 − d 3 ) S 3 ∗ = a 1 a 2 e 1 λ 1 b + c a 1 e 1 λ 1 b − a 1 e 1 λ 1 b d 2

E 6 ∗ is found based on conditions ( 2), ( 3) and ( 4) add the following condition:

Where d 4 λ 2 < e 1 λ 1 < d 3 ( a 1 − d 1 ) (7) ( c b − a 1 a 2 ) + ( a 1 e 1 λ 1 − e 1 λ 1 d 1 − d 3 ) − a 1 e 1 λ 1 b d 2 > S 4 ∗ S 4 ∗ = ( a 1 a 2 − a 1 d 2 ) e 2 λ 2 + ( a 1 c − d 1 c ) e 2 λ 2 + ( a 1 a 2 + a 1 c ) e 1 λ 1 . (8)

The stability analysis of model (1) is dealt with in this section for all stability points, and by calculating the Jacobian matrix J ( X , Y , Z ) and studying all cases of convergence and stability depending on the resulting eigenvalues to compensate for the points of model (1) that were obtained, where the matrix can be written as follows: J ( X , Y , Z ) = [ a 1 − 2 a 1 x + by − λ 1 z − d 1 bx − λ 1 x c y a 2 − 2 a 2 y + c x − λ 2 z − d 2 − λ 2 y e 1 λ 1 z e 2 λ 2 z e 1 λ 1 x + e 2 λ 2 y − d 3 − 2 d 4 z ]

System (1) contains a distinct equation for matrix close to E 0 in the following form: J ( E 0 ) = [ ( a 1 − d 1 ) 0 0 0 ( a 2 − d 2 ) 0 0 0 − d 3 ]

Given by: ( Г 1 − ( a 1 − d 1 ) ) ( Г 2 − ( a 2 − d 2 ) ) ( Г 3 + d 3 ) = 0 (9)

Consider Equation (9) based on condition (2) we get ( Г 1 > 0 ) and by ( 3) we have that ( Г 2 > 0 ) and ( Г 3 = − d 3 < 0 ) so , E 0 unstable point.

The Jacobian matrix for equilibrium point E ̂ 1 in model (1) is written as follows: J ( E ̂ 1 ) = [ − a 1 x ̂ b x ̂ − λ 1 x ̂ 0 a 2 + c x ̂ − d 2 0 0 0 e 1 x ̂ λ 1 − d 3 ]

The following equation is the characteristic equation of J ( E ̂ 1 ) and is defined as follows: ( − a 1 x ̂ − Г ̂ 1 ) ( a 2 + c x ̂ − d 2 − Г ̂ 2 ) ( e 1 x ̂ λ 1 − d 3 − Г ̂ 3 ) = 0 (10)

Where Г ̂ 1 = − a 1 x ̂ , Г ̂ 2 = a 2 + c x ̂ − d 2 , Г ̂ 3 = e 1 x ̂ λ 1 − d 3

Consider Equation (10) based on condition (3) we get ( Г ̂ 2 > 0 ) and ( Г ̂ 1 = − a 1 x ̂ < 0 )

Without the need for an eigenvalue indication Г ̂ 3 , so , E ̂ 1 (Saddle) unstable point.

The Jacobian matrix for equilibrium point E ¯ 2 in model (1) is written as follows: J ( E ¯ 2 ) = [ a 1 − d 1 + b y ¯ 0 0 c y ¯ − a 2 y ¯ − λ 2 y ¯ 0 0 e 2 λ 2 y ¯ − d 3 ]

The following equation is the characteristic equation of J ( E ¯ 2 ) and is defined as follows: ( a 1 − d 1 + b y ¯ − Г ¯ 1 ) ( − a 2 y ¯ − Г ¯ 2 ) ( e 2 y ¯ λ 2 − d 3 − Г ¯ 3 ) = 0 (11)

Consider Equation (11) based on condition (2) we get ( Г ¯ 1 > 0 ) and ( Г ¯ 2 = − a 2 y ¯ < 0 ) so , E ¯ 2 (Saddle) unstable point.

The Jacobian matrix for equilibrium point E ~ 3 in model (1) is written as follows: J ( E ~ 3 ) = [ a 1 − 2 a 1 x ~ + b y ~ − d 1 b x ~ − λ 1 x ~ c y ~ a 2 − 2 a 2 y ~ + c x ~ − d 2 − λ 2 y ~ 0 0 e 1 λ 1 x ~ + e 2 λ 2 y ~ − d 3 ]

The following equation is the characteristic equation of J ( E ~ 3 ) and is defined as follows: ( Г ~ ) 3 + A ~ 1 ( Г ~ ) 2 + A ~ 2 ( Г ~ ) + A ~ 3 = 0 (12) where A ~ 1 = d 3 − e 1 λ 1 x ~ − e 2 λ 2 y ~ + a 1 x ~ + a 2 y ~ . A ~ 2 = a 2 d 3 y ~ + a 2 e 1 λ 1 x ~ y ~ + a 2 e 2 λ 2 ( y ~ ) 2 + a 1 x ~ d 3 + a 1 e 1 λ 1 ( x ~ ) 2 − a 1 e 2 λ 2 x ~ y ~ + a 1 a 2 x ~ y ~ + cb x ~ y ~ . A ~ 3 = a 1 a 2 x ~ y ~ d 3 − e 1 a 1 a 2 x ~ y ~ − a 1 a 2 e 2 λ 2 x ~ ( y ~ ) 2 + bc e 1 λ 1 y ~ ( x ~ ) 2 + bc e 2 λ 2 x ~ ( y ~ ) 2 + bc x ~ y ~ d 3 .

According to the Routh-Hurwitz criterion, there are roots with a negative real part of Equation (12) if and only if A ~ 1 > 0 , A ~ 3 > 0 and ∆ ~ = A ~ 1 A ~ 2 − A ~ 3 > 0 .

( A ~ 1 > 0 and A ~ 3 > 0 ) based on the condition (4) add the following conditions: a 1 > e 1 λ 1 a 2 > e 2 λ 2 } (13) e 1 λ 1 x ~ > d 3 > e 1 λ 1 (14)

Now , ∆ ~ = A ~ 1 A ~ 2 − A ~ 3 .

We have that ( A ~ 1 > 0 and A ~ 3 > 0 ) with the above conditions.

And accordingly, ∆ ~ = A ~ 1 A ~ 2 − A ~ 3 > 0 . If the condition is met A ~ 1 A ~ 2 > A ~ 3 .

So, E ~ 3 it is locally asymptotical stable, otherwise unstable.

The Jacobian matrix for equilibrium point E ̇ 4 in model (1) is written as follows: J ( E ̇ 4 ) = [ M ̇ 0 0 c y ̇ − a 2 y ̇ − λ 2 y ̇ e 1 λ 1 z ̇ e 2 λ 2 z ̇ − d 4 z ̇ ] , where M ̇ = a 1 + b y ̇ − λ 1 z ̇ − d 1 .

The following equation is the characteristic equation of J ( E ̇ 4 ) and is defined as follows: ( Г ̇ ) 3 + R ̇ 1 ( Г ̇ ) 2 + R ̇ 2 ( Г ̇ ) + R ̇ 3 = 0 (15)

Where R ̇ 1 = a 2 y ̇ + d 4 z ̇ − M ̇ , R ̇ 2 = e 2 ( λ 2 ) 2 z ̇ y ̇ + a 2 y ̇ d 4 z − ̇ ( a 2 y ̇ + d 4 z ̇ ) M ̇ , R ̇ 3 = e 2 ( λ 2 ) 2 z ̇ y ̇ M ̇ − a 2 y ̇ d 4 z ̇ M ̇ .

According to the Routh-Hurwitz criterion, there are roots with a negative real part of Equation (15) if and only if R ̇ 1 > 0 , R ̇ 3 > 0 and ∆ ̇ = R ̇ 1 R ̇ 2 − R ̇ 3 > 0 .

( R ̇ 1 > 0 and R ̇ 3 > 0 ) if the following conditions hold: a 1 + b y ̇ < λ 1 z ̇ + d 1 (16) a 2 d 4 > e 2 ( λ 2 ) 2 (17)

Now , ∆ ̇ = R ̇ 1 R ̇ 2 − R ̇ 3

We have that R ̇ 1 > 0 and R ̇ 3 > 0 with the above conditions.

And accordingly, ∆ ̇ = R ̇ 1 R ̇ 2 − R ̇ 3 > 0 . If the condition is met ( R ̇ 1 R ̇ 2 > R ̇ 3 ) .

So, E ̇ 4 it is locally asymptotical stable, otherwise unstable.

The Jacobian matrix for equilibrium point E ̿ 5 in model (1) is written as follows: J ( E ̿ 5 ) = [ − a 1 x ̿ b x ̿ − λ 1 x ̿ 0 W ̿ 0 e 1 λ 1 z ̿ e 2 λ 2 z ̿ − d 4 z ̿ ] , where W ̿ = a 2 + c x ̿ − λ 2 z ̿ − d 2 .

The following equation is the characteristic equation of J ( E ̿ 5 ) and is defined as follows: ( Г ̿ ) 3 + K ̿ 1 ( Г ̿ ) 2 + K ̿ 2 ( Г ̿ ) + K ̿ 3 = 0 (18)

Where K ̿ 1 = d 4 z ̿ + a 1 x ̿ − W ̿ , K ̿ 2 = a 1 d 4 z ̿ x ̿ − W ̿ d 4 z ̿ − a 1 x ̿ W ̿ − e 1 z ̿ x ̿ ( λ 1 ) 2 , K ̿ 3 = e 1 z ̿ x ̿ ( λ 1 ) 2 W ̿ − a 1 d 4 z ̿ x ̿ W ̿ .

According to the Routh-Hurwitz criterion, there are roots with a negative real part of Equation (18) if and only if K ̿ 1 > 0 , K ̿ 3 > 0 and ∆ ̿ = K ̿ 1 K ̿ 2 − K ̿ 3 > 0 .

( K ̿ 1 > 0 and K ̿ 3 > 0 ) if the following conditions hold: a 2 + c x ̿ < d 2 + λ 2 z ̿ (19) e 1 ( λ 1 ) 2 < a 1 d 4 (20)

Now, ∆ ̿ = K ̿ 1 K ̿ 2 − K ̿ 3

We have that ( K ̿ 1 > 0 and K ̿ 3 > 0 ) with the above conditions.

And accordingly, ∆ ̿ = K ̿ 1 K ̿ 2 − K ̿ 3 > 0 If the condition is met ( K ̿ 1 K ̿ 2 > K ̿ 3 ) .

So, E ̿ 5 it is locally asymptotical stable, otherwise unstable.

The Jacobian matrix for equilibrium point E 6 ∗ in Model (1) is expressed as follows: J ( E 6 ∗ ) = [ − a 1 x ∗ b x ∗ − λ 1 x ∗ c y ∗ − a 2 y ∗ − λ 2 y ∗ e 1 λ 1 z ∗ e 2 λ 2 z ∗ − d 4 z ∗ ] .

The following equation is the characteristic equation of J ( E 6 ∗ ) and is defined as follows: ( Г ∗ ) 3 + V ∗ 1 ( Г ∗ ) 2 + V ∗ 2 ( Г ∗ ) + V ∗ 3 = 0 (21)

Where V ∗ 1 = a 1 x ∗ + a 2 y ∗ + d 4 z ∗ , V ∗ 2 = a 1 d 4 z ∗ x ∗ + a 2 d 4 z ∗ y ∗ − e 2 z ∗ y ∗ ( λ 2 ) 2 + a 1 a 2 y ∗ x ∗ + b c y ∗ x ∗ − e 1 z ∗ x ∗ ( λ 1 ) 2 , V ∗ 3 = ( b e 1 − c e 2 ) λ 1 λ 2 z ∗ x ∗ y ∗ + ( a 2 d 4 – e 2 λ 2 ) a 1 z ∗ x ∗ y ∗ + ( b c d 4 − a 2 e 1 ( λ 1 ) 2 ) z ∗ x ∗ y ∗ .

According to the Routh-Hurwitz criterion, there are roots with a negative real part of Equation (21) if and only if V ∗ 1 > 0 , V ∗ 3 > 0 and ∆ ∗ = V ∗ 1 V ∗ 2 − V ∗ 3 > 0 .

( V ∗ 1 > 0 and V ∗ 3 > 0 ) based on the condition (17) add the following conditions: b e 1 > c e 2 (22) b c d 4 > a 2 e 1 ( λ 1 ) 2 (23)

Now, ∆ ∗ = V ∗ 1 V ∗ 2 − V ∗ 3 .

We have that ( V ∗ 1 > 0 and V ∗ 3 > 0 ) with the above conditions.

And accordingly, ∆ ∗ = V ∗ 1 V ∗ 2 − V ∗ 3 > 0 . If the condition is met ( V ∗ 1 V ∗ 2 > V ∗ 3 ) .

So, ( E 6 ∗ ) it is locally asymptotical stable, otherwise unstable.

The comprehensive study and analysis in particular is studied in this section for the balanced points found in the previous section, called the local stability of model (1), and that is by using the Lyapunov function that satisfies the conditions of global equilibrium, as shown in the following theorems: Lemma:

(Lyapunov stability theorem):

Suppose that x ∗ be an equilibrium point of system (1.12) and V : N → R be an C ¹ real-valued function defined on some neighborhood N of x ∗ with ( N ⊆ R n ) such that ⋄

V ( x ∗ ) = 0 and V ( x ) > 0 for all x ≠ x ∗ and x ∈ N (i.e., mean V is a positive definite function).

Note that the function V given above is known as the Lyapunov function (weak Lyapunov function) if and only if the first and second conditions hold, whereas it is known as a strict Lyapunov function (strong Lyapunov function) if and only if the first and third conditions hold. In addition, if N can be chosen to be all of R n , then x ∗ is said to be globally asymptotically stable if the first and third conditions hold.