www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Qualitative analysis of a mathematical model about population of green turtles on the Galapagos island Candy Herrera∗, Cosme Duque† and Hugo Leiva‡ ∗ Yachay Tech University, School of Biological Science and Engineering San Miguel de Urcuqui-100119. Imbabura-Ecuador candy.herrera@yachaytech.edu.ec † Universidad de Los Andes, Facultad de Ciencias, Departamento de Matemáticas Mérida 5101-Venezuela cosduq@gmail.com, duquec@ula.ve ‡Yachay Tech University, School of Mathematical and Computational Sciences San Miguel de Urcuqui-100119. Imbabura, Ecuador hleiva@yachaytech.edu.ec Received: 8 March 2021, accepted: 29 July 2021, published: 1 October 2021 Abstract— According to the IUCN, most sea turtles fall into one of the endangered categories. Since, sea turtles, like many other reptiles, present an unusual developmental process, marked by the determination of the sex of the offspring by envi- ronmental factors, more specifically by temperature. In the temperature sex determination (TSD) system the temperature of an embryo’s environment during incubation period will dictate the embryo’s sex development. This developmental process, together with the complex mating and nesting behavior and the vulnerability of sea turtles to threats of a natural or anthropogenic nature, naturally lead to the study of the population dynamics of the species. For this reason, in this paper, we have developed a contin- uous model given by a system of three ordinary differential equations to study the dynamics of the green sea turtle population long-term, focusing the mathematical simulations on the data obtained for the nesting species of Galapagos Islands. Through the qualitative analysis of the model, the following is demonstrated: 1) The flow induced by the system is positively invariant on the region of biological inter- est (Ω); and 2) The given condition on f̂ is necessary and sufficient for the unique nontrivial equilibrium point (I∗) to be globally asymptotically stable in that region. When implementing the estimated values for our parameters in the numerical simulations, it was observed that indeed the population of Galapagos green sea turtles complies with the condition for which the nontrivial critical point (I∗) is globally asymptotically stable; that is, the asymptotic stabil- Copyright: © 2021 Herrera et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Candy Herrera, Cosme Duque, Hugo Leiva, Qualitative analysis of a mathematical model about population of green turtles on the Galapagos island, Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 1 of 10 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... ity is maintained for any initial value within the set Ω. In contrast, when altering the estimated values of the parameters so that the established condition is not met, the trivial critical point (I0) becomes globally stable, and the population falls towards extinction regardless of the values taken within the positively invariant Ω set. Keywords-population dynamics; sex-structured continuous time model; Chelonia mydas; equilibrium point; local stability; global stability. I. INTRODUCTION Sea turtles present a temperature sex determina- tion (TSD) system, this is an evolutionary condi- tion that turtles, like other reptiles, have adopted throughout time. In the TSD system, the tempera- ture of an embryo’s environment during incubation period will dictate the embryo’s sex development. Many species of turtles [1], tortoises [2], lizards [3], and crocodiles [4] that exhibit TSD have a thermosensitive period (TSP) during which the embryo sex is developed. For turtles, this period has been observed to take place during the mid- trimester of the embryo incubation period [5]. The temperature that defines the 1:1 sex ratio balance, known as pivotal temperature, is ∼ 29.4◦C [6], [7]. When the mean temperature of the nest during TSD is around the pivotal temperature an even distribution of male and female hatchlings occur [6], [7]. Below the pivotal temperature, hatchling sex population will be mostly male; and above, it will be mostly female[6], [7], [8]. With a trend towards increasing mean global temperature, species with TSD are particularly af- fected. The population of sea turtles is facing high egg mortality and feminization of the offspring[9], [10]. In recent years a disproportionate ratio of female to male turtle eggs has been observed in a number of different studies [11]. Unfortunately, monitoring sex ratios involves a series of method- ological and ethical complications. Sex ratios data obtained from the study of adult populations have predicted a complete feminization and a possibly extinction of marine turtles in the future if tem- perature continues rising [10], [7], [12], [13]. Due to the complex nesting behavior of sea turtles and the mechanism by which the sex of the offspring is determined, it is reasonable to look into the factors that could be directly affecting the sex ratios of eggs, such as sand temperatures. In order to understand why there is a bias towards the fe- male population in sea turtles. The most prevalent hypothesis points at climate change as the main factor leading the rising at sand temperatures in nesting sites [11], [10]. The lack of males within the sea turtle population will eventually affect the population dynamics. It is unknown the minimum proportion of males sufficient to support the sea turtle population in such a way to avoid population collapse or even extinction. This problem naturally lends itself to investigation via a sex-structured model to analyse the dynamic of the population. II. MODEL FORMULATION The green turtle, like other sea turtles, has a remarkable life cycle. Individuals inhabit widely separated localities during the course of their lives. These habitats include foraging, migration, breed- ing, and nesting areas (Figure 1). After hatchlings emerge from their nests, they immediately travel to the sea. Once in the ocean, the hatchlings are washed away by ocean currents, live a pelagic phase in the open ocean, and are not seen again until they appear as juveniles in foraging areas, probably a decade later [14]. This period of time between the hatchling and juvenile stages is known as the ”lost years” because the migratory route taken by the hatchlings and their behavior remains a mystery [14], [15]. In the foraging areas, they continue to mature until they become subadults. Subadults are occasionally seen foraging in the open ocean [14]. Once the turtles reach sexual maturity, males and females return to their home beach to mate and nest. This behavior is known as philopatry and has been documented for several species [19]. Because during the early stages of the turtle life cycle, hatchlings, juveniles, and subadults are not reproductively active, we can simplify the life cycle of the green turtle and consider only two main stages: adults and eggs. In turn, the adult stage can be divided into two populations: female Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 2 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... Fig. 1. Schematic diagram of the Chelonia mydas life cycle. adults and male adults, which are interesting to analyze by themselves. An schematic diagram of the simplified life cycle of green sea turtle is shown in Figure 2. This scheme of the life cycle (Figure 2) makes a clear distinction between adult males and females, focusing on their interaction during mating pro- cess that eventually will result in the production of eggs. After an incubation period of 7 to 9 weeks [17] the eggs develop into either female or male sea turtles. It also shows the outflows of each stage: the death of adult females, males, and eggs respectively. Successful mating of males and females will produce and increase the egg population which in turn will develop into males and females. What determines the proportion of eggs that develop as males or females is the incubation temperature of the eggs. The reproductive biology and nesting behavior of sea turtles comprises a number of aspects that make modeling this interaction truly challenging. The interaction itself between males and females involves a variety of factors that can strongly affect the birth rate of the population and therefore the long-term dynamics. Density of adult female (AF ) and male populations (AM ), the behavioral responses during mating process, and searching efficiency [18], are some of the variables to take into account for a successful mating. This study is particularly concerned on the proportion of eggs allocated to males and females in the population. At time t, we denote the egg population as E(t) and adult population for males and females as AM (t), AF (t) respectively. The dynamics of the egg population is governed by the following first order ordinary differential equation: dE dt = g(AF ,AM ) − (α + µe)E, (1) where α is the maturity rate of eggs that become adult males or females. µe is the mortality rate of the egg stage. The function g(AF ,AM ) is the recruitmant rate of eggs and is represented by: g(AF ,AM ) = f̂(AM + AF ) ( 1 − AM + AF K ) . For this particular case, function g is dependent on the size of the population, so it is defined by a logistic type function, where f̂ is the fecundity rate, K is the environmental carrying capacity, the maximum population size that the environment can sustain indefinitely. The saturation constraint shows that when the adult population approaches the carrying capacity, then the egg per-capita ap- proaches zero. Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 3 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... Fig. 2. Schematic diagram of the mathematical model. For changes in adult male and female population over time linear ordinary differential equations were developed: dAM dt = qαEM −µAM, (2) dAF dt = (1 −q)αEF −µAF , (3) where the proportion of male and female eggs are given by q and 1−q receptively. µ tell us the death rate for the adult stage. As there is not evidence of different maturity and death rates for males and females, α and µ remains the same for both sexes. From (1) to (3) the next system of differential equations is proposed:   dAM dt = qαE −µAM ≡ F1(AM,AF ,E) dAF dt = (1 −q)αE −µAF ≡ F2(AM,AF ,E) dE dt = f̂ (AF + AM ) ( 1 − AM + AF K ) − (α + µe) E ≡ F3(AM,AF ,E) (4) III. QUALITATIVE ANALYSIS OF THE MODEL In this section we will give a complete quali- tative description of the dynamics of system (4), concretely, we characterize the region where the system is positively invariant and we shall describe completely the global stability of the nontrivial equilibrium point. Since the vector field of system (4) is con- tinuously differentiable, our first result, follows from the fundamental existence-uniqueness theo- rem (see for instance L. Perko [24], Hale [25]). Theorem 3.1: For any initial condition AM (0) ≥ 0, AF (0) ≥ 0, and E(0) ≥ 0, there exists β > 0 such that system (4) has a unique solution defined on [0,β). The next theorem guarantees that the system (4) is biologically well behaved and that the dynamic of the system is concentrated on a bounded region of R3+. Concretely, the following results holds: Theorem 3.2: Suppose that f̂ < 2µ(α+µe) α , then the region Ω defined by Ω = { (AM,AF ,E) ∈ R3 : 0 < AF < K, 0 < AM < K −AF , 0 < E < µK 2α } (5) is positively invariant under the flow induced by (4) (see Figure 3). Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 4 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... Fig. 3. Region where (4) is positively invariant Proof: In order to proof this theorem, we will analyze the vector field defined by (4) on ∂Ω, this analyze is made in the following table. Note that the prism given by Ω has five faces and six vertices. i) AM ∈ (0,K), AF = 0, E = 0 A′M = −µAM < 0, A′F = 0, E′ = f̂ ( 1 − AM K ) AM > 0 ii) AM = 0, AF ∈ (0,K), E = 0 A′M = 0, A′F = −µAF < 0 E′ = f̂ ( 1 − AF K ) AF > 0 iii) AM = 0, AF = 0, E ∈ (0, µK2α ) A′M = qαE > 0 A′F = (1 −q)αE > 0 E′ = −(α + µe)E < 0 iv) AM = K, AF = 0, E ∈ (0, µK2α ) A′M = −µK + qαE ≤− µK 2 < 0, A′F = (1 −q)αE > 0, E′ = −(α + µe)E < 0, v) AM = 0, AF = K, E ∈ (0, µK2α ) A′M = −qαE > 0, A′F = −µK + (1 −q)αE ≤− µK 2 < 0 E′ = −(α + µe)E < 0, vi) AM = 0, AF ∈ (0,K), E = µK2α A′M = qµK 2 > 0, A′F = no matter sign E′ = −(α + µe) µK2α + f̂ ( 1 − AF K ) AF ≤−(α + µe) µK2α + f̂ K 4 < 0 vii) AM ∈ (0,K), AF = 0, E = µK2α A′M = no matter sign A′F = (1 −q) µK 2 > 0 E′ = −(α + µe) µK2α + f̂ ( 1 − AM K ) AM ≤−(α + µe) µK2α + f̂ K 4 < 0 viii) AM ∈ (0,K), AF = K −AM , E = 0 A′M = −µAM < 0, AF = −µAF < 0, E′ = 0 ix) AM ∈ (0,K), AF = K −AM , E = µK2α Note that 〈1, 1, 0〉 is a normal vector of the plane AM + AF = K. Therefore, 〈F1,F2,F3〉 · 〈1, 1, 0〉 = −µK2 < 0, and E′ = −(α + µe) µK2α < 0. So, 〈F1,F2,F3〉 is directed to the interior of Ω x) AM ∈ (0,K), AF = 0, E ∈ (0, µK2α ) A′M = no matter sign A′F = (1 −q)αE > 0 E′ = no matter sign xi) AM = 0, AF ∈ (0,K), E ∈ (0, µK2α ) A′M = qαE > 0, A′F = no matter sign, E′ = no matter sign. xii) AM ∈ (0,K), 0 < AF < K −AM , E = 0 A′M = no matter sign A′F = no matter sign E′ = f̂ ( 1 − AM +AF K ) (AM + AF ) > 0 xiii) AM ∈ (0,K), 0 < AF < K −AM , E = µK2α A′M = no matter sign A′F = no matter sign E′ = −(α + µe) µK2α + f̂ ( 1 − AM +AF K ) (AM + AF ) ≤−(α + µe) µK2α + f̂ K 4 < 0 xiv) AM ∈ (0,K), AF = K −AM , E ∈ (0, µK2α ) This case is analogous to (ix) Therefore, we have that the vector field given by the right side of system (4) on the boundary of Ω is directed to the interior of the set Ω, in consequence solutions with initial data in Ω remain there for any t ≥ 0. This conclude the proof. � Remark 1: From Theorem 3.2 we deduced that the solutions of system (4) are defined for all t ≥ 0. Hereafter we will assume that f̂ < 2µ(α+µe) α . The equilibria of system (4) are obtained by the solutions of the following algebraic equations   qαE −µAM = 0 (1 −q)αE −µAF = 0( 1 − AM + AF K ) f̂ (AF + AM ) − (α + µe) E = 0 (6) So, the equilibria of (4) consist of one triv- ial critical point I0 = (0, 0, 0), that always ex- ists, and a unique nontrivial critical point I∗ = Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 5 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... (A∗M,A ∗ F ,E ∗) which exists if and only if the following condition is true f̂ > µ(α + µe) α (7) In this case, we have A∗M =qK ( 1 − µ (α + µe) f̂α ) , A∗F =(1 −q)K ( 1 − µ (α + µe) f̂α ) , E∗ =µK ( 1 − µ (α + µe) f̂α ) . (8) A. Local stability of equilibrium points In this subsection we shall discuss the local stability properties of the equilibria I0 and I∗. Theorem 3.3: The equilibrium point I0 is lo- cally asymptotically stable if f̂ < µ(α+µe) α Proof: The Jacobian matrix of system (4) about the equilibrium point I0 is given by J(I0) =   −µ 0 qα 0 −µ (1 −q) α f̂ f̂ −α−µe   (9) and the characteristic polynomial associated to (9) is P(λ) = a0λ 3 + a1λ 2 + a2λ + a3, (10) where a0 = 1, a1 = α + µe + 2µ, a2 = (−f̂α + µ(α + µe)) + µ(α + µ + µe) and a3 = µ(−αf̂ + µ(α + µe)). Since a1, a2 and a3 are positive and a1a2 −a3 = (α + µe + 2µ)(−f̂α + µ(α + µe)) + µ(α + µ + µe)−µ(−αf̂ + µ(α + µe)) > 0, so by the Routh-Hurwitz criterion (see [16]) we have that the equilibrium point I0 is locally asymptotically stable. � Theorem 3.4: The equilibrium point I∗ is local asymptotically stable if and only if f̂ > µ(α+µe) α . Proof: The Jacobian matrix of system (4) about the equilibrium point I∗ is given by J(I∗) =   −µ 0 qα 0 −µ (1 −q) α γ∗ γ∗ −α−µe   , (11) where γ∗ = 2µ(α + µe) α − f̂. After a simple calculation, one finds the char- acteristic polynomial is P(l) = a0λ 3 + a1λ 2 + a2λ + a3, (12) where a0 = 1, a1 = µe + α + 2µ, a2 = f̂α + µ2, and a3 = µ(f̂α−µ(α + µe)). Since a0, a1 and a3 are positive and a1a2−a3 = (µe + α + 2µ)(f̂α + µ 2) −µ(f̂α−µ(α + µe)) > 0, then we have by the Routh-Hurwitz criterion (see [16]) that the equilibrium point I∗ is locally asymptotically stable for system (4). Therefore I∗ is locally asymptotically stable if and only if it exist, i.e. if and only if f̂ > µ(α+µe) α . � B. Global stability analysis Now, we are ready to prove that under some conditions the non trivial equilibrium point is globally asymptotically stable, this is done trans- forming system (4) into a planar system and ap- plying Bendixon’s criterion first, and next Poincare Bendixon’s Theorem. After that, we go back to system (4) and prove the global asymptotic stabil- ity of the three dimensional equilibrium point. Theorem 3.5: If f̂ > µ(α+µe) α , then the equi- librium point I∗ is global asymptotically stable in Ω; otherwise I0 is global asymptotically stable. Proof: Let (AM (0),AF (0),E(0)) be an arbi- trary initial data on Ω. If we make the change of Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 6 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... variable A = AM + AF , then system (4) can be reduced to the following bidimensional system   dA dt = αE −µA := G1(A,E) dE dt = f̂ ( 1 − A K ) A − (α + µe) E := G2(A,E). (13) Since the vector field of system (13) is directed to the interior of the region Ω̃ = { (A,E) ∈ R2 : 0 0 such that E(t) > E∗/2 for t > τ; hence∫ t 0 qαE(s)eµsds= ∫ τ 0 qαE(s)eµsds+ ∫ t τ qαE(s)eµsds ≥ ∫ τ 0 qαE(s)eµsds+ ∫ t τ qα E∗ 2 eµsds−→∞ as t→∞. Therefore, by using the L’Hospital rule, we have that lim t→∞ AM (t) = lim t→∞ AM (0)e −µt + lim t→∞ e−µt ∫ t 0 qαE(s)eµsds qα µ E∗= A∗M. Analogously, we obtain that lim t→∞ AF (t) = (1 −q)α µ E∗ = A∗F . So, the equilibrium point I∗ = (A∗M,A ∗ F ,E ∗) is globally asymptotically stable in Ω. If µ ≥ f̂α α+µe , then I0 is the unique point of equilibrium of system (4). Following an analogous reasoning to the one above, we have that I0 is globally asymptotically stable. This conclude the proof. � IV. NUMERICAL SIMULATION In this section is devoted to show numerical examples that illustrate our results. In this sense we will use Table I, and several initial conditions. If let us pick f̂ = 0.25 and K = 1, then the nontrivial critical point I∗ = (0.24385, 0.24385, 0.025361) is a global attractor (see Figure 5) Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 7 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... TABLE I PARAMETERS DESCRIPTION AND VALUES CONSIDERING THE GALAPAGOS ISLANDS AND ECUATORIAN MAINLAND GREEN SEA TURTLE POPULATIONS. Parameter Description Estimated values Reference α Maturity rate of green turtle 0.54 ±0.45 [20] µ Per capita death rate for adults 0.052 ± 0.005 [21] µe Per capita death rate for eggs 0.79 ± 0.19 [22] q Proportion of eggs that become male 0.5 [23] f̂ Fertility rate 0.35 ± 0.15 (Estimated) Fig. 5. Solutions of system (4) with f̂ = 0.25 and K = 1 If let us pick f̂ = 0.01 and K = 1, then the trivial critical point I0 is the global attractor (see Figure 6). Fig. 6. Solutions of system (4) with f̂ = 0.01 and K = 1 V. CONCLUSION AND FINAL REMARK The model developed in this paper is character- ized by studying the population dynamics through a sex-structured continuous time model. Unlike conventional discrete age-structured models, in this paper we seek, through the distinction between sexes, to determine how the different parameters described influence population dynamics and its stability. From the above results of the qualitative analysis of the model, we infer that the condition established determine the persistence and stability of the green turtle population or its extinction. In terms of the biological significance of the proposed condition, we can deduce that the number of offspring a female individual produce over time should be biger than the outflow rates of the pop- ulation stages. This model can be used as a basis for a more complex models and include a series of parameters that can offer us a better appreciation of population dynamics and how it is affected or distorted by the influence of new parameters, such as temperature. Since the proportion of males and females is directly related to the incubation temperature of the eggs, modeling the population dynamics by considering the proportion of male eggs as a function of temperature q(T) and adding details about other stages of the turtle life cycle could reveal new conditions, which help to envi- sion the population situation and take actions for the conservation of the species. The methodology applied here to study the dynamics of the popula- tion of green turtles on the Galapagos island, can be applied to study the dynamics of other types of Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 8 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... reptiles whose reproduction also depends on the environment and temperature changes. In all these cases, the region of biological interest will be of the polyhedral type, which allows us to study the vector field on each face of the region, and the invariance of this region can be achieved in the same way we did it with the prism in this work. After that, we can study the local stability of the equilibrium points by linearizing the differential equation around them, and looking at the sign of the eigenvalues; then reduce the system to a two- dimensional one and apply Bendixon’s criterion, and finish it with the Poincare-Bendixon theory to conclude that the non-trivial equilibrium point is globally asymptotically stable. From the foregoing comments, it is clear that our method can be used to study a broad class of similar problems. ACKNOWLEDGMENT We would like to thank the anonymous reviewer who carefully read the paper, his/her comments and suggestions allowed us to improve the final presentation of this manuscript. REFERENCES [1] S.J. Morreale, G.J. Ruiz and E.A. Standora, Temperature-dependent sex determination: current practices threaten conservation of sea turtles. Science, 216 (4551), 1245-1247, 1982. [2] V. Lewis-Winokur and R.M. Winokur, Incubation tem- perature affects sexual differentiation, incubation time, and posthatching survival in desert tortoises (Gopherus agassizi), Canadian Journal of Zoology, 73 (11), 2091- 2097, 1995. [3] J.J. Bull, Temperature-sensitive periods of sex determi- nation in a lizard: Similarities with turtles and crocodil- ians, Journal of Experimental Zoology, 241 (1), 143- 148, 1987. [4] D.C. Deeming and M.W. Ferguson, The mechanism of temperature dependent sex determination in crocodil- ians: a hypothesis, American Zoologist, 29 (3), 973-985, 1989. [5] M.M Fuentes, M. Hamann, M. and C.J. Limpus, Past, current and future thermal profiles of green turtle nest- ing grounds: Implications from climate change, Journal of Experimental Marine Biology and Ecology, 383 (1), 56-64, 2010. [6] Y. Matsumoto, and D. Crews, Molecular mechanisms of temperature-dependent sex determination in the context of ecological developmental biology, Molecular and cellular endocrinology, 354 (1-2), 103-110, 2012. [7] S.Y. Özdilek, B.E. Sönmez, and Y. Kaska, Sex ratio estimations of Chelonia mydas hatchlings at Samandağ Beach, Turkey, Turkish Journal of Zoology, 40 (4), 552- 560, 2016. [8] L.I. Wright, K.L. Stokes, W.J. Fuller, B.J. Godley, A. McGowan, R. Snape and A.C. Broderick, Turtle mating patterns buffer against disruptive effects of climate change, Proceedings of the Royal Society B: Biological Sciences, 279 (1736), 2122-2127, 2012. [9] N. Mrosovsky and C.L. Yntema, Temperature depen- dence of sexual differentiation in sea turtles: implica- tions for conservation practices, Biological Conserva- tion, 18 (4), 271-280, 1980. [10] M.P. Jensen, C.D. Allen, T. Eguchi, I.P. Bell, E.L. LaCasella, W.A. Hilton and P.H. Dutton, Environmental warming and feminization of one of the largest sea turtle populations in the world, Current Biology, 28 (1), 154- 159, 2018. [11] S. Yamaguchi and Y. Iwasa, Temperature-dependent sex determination, realized by hormonal dynamics with enzymatic reactions sensitive to ambient temperature, Journal of theoretical biology, 453, 146-155, 2018. [12] R. King, W.H. Cheng, C.T. Tseng, H. Chen and I.J. Cheng, Estimating the sex ratio of green sea turtles (Chelonia mydas) in Taiwan by the nest temperature and histological methods, Journal of experimental marine biology and ecology, 445, 140-147, 2013. [13] J.R. Spotila, E.A. Standora, S.J. Morreale and G.J. Ruiz, Temperature dependent sex determination in the green turtle (Chelonia mydas): effects on the sex ratio on a natural nesting beach, Herpetologica, 74-81, 1987. [14] A.B. Bolten, and G.H. Balazs, Biology of the early pelagic stage—the ‘lost year.’. Biology and Conservation of Sea Turtles, Revised edition, Smithsonian Institute Press, Washington, DC, 579, 1995. [15] A. Carr, Notes on the behavioral ecology of sea turtles. Biology and conservation of sea turtles, Smithsonian Institution press Washington, DC, 19-23, 1982. [16] F. Gantmacher, The Theory of Matrices, Vol 2, Chelsea Publishing, New York, 1974. [17] D. Green, The east Pacific green sea turtle in Galapa- gos, Noticias de Galápagos, Charles Darwin Research Station (28), 9-12, 1978 [18] H.F. Hirth and D.A. Samson, Nesting behavior of green turtles (Chelonia mydas) at Tortuguero, Costa Rica. Comportamiento de anidamiento de las tortugas verde (Chelonia mydas) en Tortuguero, Costa Rica, Caribbean Journal of Science, 23 (3/4), 374-379, 1987. [19] P.L. Lee, P. Luschi and G.C. Hays, Detecting female precise natal philopatry in green turtles using assign- ment methods, Molecular ecology, 16(1), 61-74, 2007. [20] D. Green, Growth rates of wild immature green turtles in the Galapagos Islands, Ecuador, Journal of Herpetol- ogy, 338-341, 1993. [21] J.J. Alava, P. Jiménez, M. Peñafiel, W. Aguirre and P. Amador. Sea turtle strandings and mortality in Ecuador: 1994-1999, Marine Turtle Newsletter, 108, 4-7, 2005. Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 9 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 C Herrera, C Duque, H Leiva, Qualitative analysis of a mathematical model about population of ... [22] P. Zárate, P, K.A. Bjorndal, M. Parra, P.H. Dutton, J.A. Seminoff, and A.B. Bolten, Hatching and emergence success in green turtle Chelonia mydas nests in the Galápagos Islands, Aquatic Biology, 19 (3), 217-229, 2013. [23] S. Zavala Montoya, J. Belmont, M. Hirschfeld and D. Alarcón, Estimación de la proporción de sexos de lator- tuga verde (Chelonia mydas) en áreas de alimentación en las islas Galápagos, In Simposio de Tortugas Mari- nas de Ecuador 2018. [24] L. Perko, Differential Equations and Dynamical Sys- tems, Texts in Applied Mathematics 7, Springer-Verlag, 1993. [25] J.K. Hale, Ordinary Differential Equations, Dover Pub- lications, Mineola, New York, 2009. Biomath 10 (2021), 2107293, http://dx.doi.org/10.11145/j.biomath.2021.07.293 Page 10 of 10 http://dx.doi.org/10.11145/j.biomath.2021.07.293 Introduction Model Formulation Qualitative Analysis of the Model Local stability of equilibrium points Global stability analysis Numerical Simulation Conclusion and Final Remark References