https://doi.org/10.53656/math2025-1-2-ano

2025/1, стр. 18 - 34

A NOTE ON A GENERALIZED DYNAMICAL SYSTEM OCCURS IN MODELLING “THE BATTLE OF THE SEXES”: CHAOS IN SOCIOBIOLOGY

Резюме: In this paper, we propose a new system that occurs in modelling “the battle of the sexes” in evolutionary biology (Hofbauer & Sigmund 1988). The existence of a heteroclinic cycle and a continuous family of periodic orbits of the system is established; then the dynamical characteristics of a time-periodic perturbation of the system are investigated. By using the well-known Melnikov’s method, a sufficient condition is obtained for the perturbed system to have a transverse hetero-clinic cycle and hence to possess chaotic behaviour in the sense of Smale. One possible application that Melnikov functions may find in the modelling and synthesis of radiating antenna patterns is considered. We demonstrate some modules for investigating the dynamics of the proposed model. This will be included as an integral part of a planned much more general Web-based application for scientific computing. The proposed new extended model contains many free parameters (the coefficients \(a_{i}, i=1,2, \ldots, N\) ), which makes it attractive for use in the fields of biological applications, chemistry, sociology, lifetime analysis, reaction kinetics, biostatistics, population dynamic, medical research, games theory etc. Finally, a special case of subharmonic solutions is discussed.

Ключови думи: new system occurs in modelling “the battle of the sexes”; hetero-clinic orbit; Melnikov function; subharmonic Melnikov function; chaotic behavior

1. Introduction

A number of authors devote their research to the classical differential system that occurs in modelling “the battle of the sexes” in evolutionary biology (Hofbauer & Sigmund 1988; Dawkins 1976; Foster & Young 1990; Schuster & Sigmund 1981; Taylor & Jonker 1978). The publications on this topic are significant and varied. It is known that the dynamical game theoretic mating behaviour of males and females can be modelled by a planar system of autonomous ordinary differential equations. This system occurs in modelling “the battle of the sexes” in evolutionary biology. In (Christie et al. 1995) the authors presented the following model

(1)\[ \left\{\begin{array}{l} \tfrac{d x}{d t}=x(1-x)(\alpha-\beta y)+\epsilon\left(a x+a_{1} \sin (\omega t)\right) \\ \tfrac{d y}{d t}=y(1-y)(\gamma-\delta x)+\epsilon\left(a y+a_{1} \sin (\omega t)\right) \end{array}\right. \]

where \(0 \leq \epsilon \lt 1\).

The unperturbed differential system is non-Hamiltonian under general assumptions, becoming Hamiltonian only for \(\{\alpha, \beta=2 \alpha, \gamma=-\alpha, \delta=-2 \alpha\}\).

The existence of a hetero-clinic cycle and a continuous family of periodic orbits of the system (1) is established. In the serious studies (Hofbauer & Sigmund 1988; Dawkins 1976; Foster & Young 1990; Schuster & Sigmund 1981; Taylor & Jonker 1978; Christie et al. 1995) cited above, the reader can find a considerable volume of literature devoted to this classic model. The reader can find important questions related to the topic Strategies and stability: an opening in game dynamics, some aspects of sociobiology, evolutionarily stable strategies, game dynamics and asymmetric conflicts in the monograph (Hofbauer & Sigmund 1988). The genetic model for the “battle of the sexes“ is due to (Smith & Hofbauer 1987).

In this paper, we suggest a new class of extended asymmetric models of the type (1) that occurs in modelling “he battle of the sexes” in evolutionary biology. Investigations in the light of Melnikov’s approach is considered. A sufficient condition is obtained for the perturbed system to have a transverse heteroclinic cycle and hence to possess chaotic behaviour in the sense of Smale. Several simulations are composed. We demonstrate some modules for investigating the dynamics of the proposed model. This will be included as an integral part of a planned much more general Webbased application for scientific computing.

The plan of the paper is as follows. We state our model in Section 2. Investigations in the light of Melnikov’s approach is considered in Section 3. Some simulations are presented in Section 4. One possible application that Melnikov functions may find in the modelling and synthesis of radiating antenna patterns is considered in Section 5. We conclude by Section 6.

2. The new model

In this paper, we suggest a modified model of the type:

(2)\[ \left\{\begin{array}{l} \tfrac{d x}{d t}=x(1-x)(\alpha-\beta y)+\epsilon\left(a x+\sum_{j=1}^{N} a_{j} \sin (j \omega t)\right) \\ \tfrac{d y}{d t}=y(1-y)(\gamma-\delta x)+\epsilon\left(a y+\sum_{j=1}^{N} a_{j} \sin (j \omega t)\right) \end{array}\right. \]

where \(0 \leq \epsilon \lt 1, a_{i} \geq 0 ; i=1,2, \ldots, N\) and \(N\) is integer.

3. Considerations in the light of Melnikov’s approach The system (2) is of the form:

\[ \begin{aligned} \tfrac{d x}{d t} & =f_{1}(x, y)+\epsilon g_{1}(x, y, t) \\ \tfrac{d y}{d t} & =f_{2}(x, y)+\epsilon g_{2}(x, y, t) \end{aligned} \]

Let \(q_{0}(t)=\left(x_{0}(t), y_{0}(t)\right)\) be one of the hetero-clinic solutions of the first integral of (2) \((\epsilon=0)\) given by (see for more details precise considerations in (Christie et al. 1995))

(3)\[ \begin{aligned} & O A: y_{0}(t)=0 ; x_{0}(t)=\tfrac{e^{\tfrac{\alpha t}{2}}}{2 \cosh \left(\tfrac{\alpha t}{2}\right)} \\ & A B: x_{0}(t)=1 ; y_{0}(t)=\tfrac{e^{\tfrac{(\gamma-\delta) t}{2}}}{2 \cosh \left(\tfrac{(\gamma-\delta) t}{2}\right)} \\ & B C: y_{0}(t)=1 ; x_{0}(t)=\tfrac{e^{\tfrac{(\alpha-\beta) t}{2}}}{2 \cosh \left(\tfrac{(\alpha-\beta) t}{2}\right)} \\ & C O: x_{0}(t)=0 ; y_{0}(t)=\tfrac{e^{\tfrac{\gamma t}{2}}}{2 \cosh \left(\tfrac{\gamma t}{2}\right)} \end{aligned} \]

(Here \(O(0,0), A(1,0), B(1,1), C(1,1)\) are the other fixed points.)

See also fig. 1 for fixed \(\alpha=1, \beta=2, \gamma=-1, \delta=-2\).

The Melnikov function (Melnikov 1963) corresponding to the system (2) is of the form:

(4)\[ \begin{aligned} & M\left(t_{0}\right)=\int_{-\infty}^{\infty} f\left(q_{0}(t)\right) \wedge g\left(q_{0}(t), t+t_{0}\right) e^{-\int_{0}^{t} \operatorname{traceDf}\left(q_{0}(s)\right) d s} d t \\ & =\int_{-\infty}^{\infty}\left(x_{0}(t)\left(1-x_{0}(t)\right)\left(\alpha-\beta y_{0}(t)\right)\left(a y_{0}(t)+\sum_{j=1}^{N} a_{j} \sin \left(j \omega\left(t+t_{0}\right)\right)\right)\right. \\ & \left.-y_{0}(t)\left(1-y_{0}(t)\right)\left(\gamma-\delta x_{0}(t)\right)\left(a x_{0}(t)+\sum_{j=1}^{N} a_{j} \sin \left(j \omega\left(t+t_{0}\right)\right)\right)\right) B(t) d t \end{aligned} \]

where

\[ B(t)=e^{-\int_{0}^{t}\left(\left(1-2 x_{0}(s)\right)\left(\alpha-\beta y_{0}(s)\right)+\left(1-2 y_{0}(s)\right)\left(\gamma-\delta x_{0}(s)\right)\right) d s} \] The Melnikov function can be represented in another way, using the following familiar equality

\[ \sin \left(i \omega\left(t+t_{0}\right)\right)=\sin (i \omega t) \cos \left(i \omega t_{0}\right)+\cos (i \omega t) \sin \left(i \omega t_{0}\right) \]

We will not dwell on this question here.

For us, it is more important to note that such a representation is more appropriate and can be used directly by users of the corresponding specialized module implemented for example in CAS Mathematica.

If \(M\left(t_{0}\right)=0\) and \(\tfrac{M\left(t_{0}\right)}{d t_{0}} \neq 0\) for some \(t_{0}\) and some sets of parameters, then chaos occurs.

The Melnikov function gives a measure of the leading order distance between the stable and unstable manifolds when \(\epsilon \neq 0\).

3.1. The case \(N=1, \alpha=1, \beta=2, \gamma=-1, \delta=-2\)

From (4) we find:

\[ \begin{gathered} \left.M\left(t_{0}\right)\right|_{O A}=a_{1} \pi \omega \operatorname{csch}(\pi \omega) \sin \left(\omega t_{0}\right) \\ \left.M\left(t_{0}\right)\right|_{A B}=-a-a_{1} \pi \omega \operatorname{csch}(\pi \omega) \sin \left(\omega t_{0}\right) \\ \left.M\left(t_{0}\right)\right|_{B C}=-a-a_{1} \pi \omega \operatorname{csch}(\pi \omega) \sin \left(\omega t_{0}\right) \\ \left.M\left(t_{0}\right)\right|_{C O}=a_{1} \pi \omega \operatorname{csch}(\pi \omega) \sin \left(\omega t_{0}\right) \end{gathered} \]

Proposition 1. If \(N=1, \alpha=1, \beta=2, \gamma=-1, \delta=-2\), then the roots of the Melnikov function \(M\left(t_{0}\right)\) are given as solutions of the equation

\[ M\left(t_{0}\right)=a+a_{1} \pi \omega c s c h(\pi \omega) \sin \left(\omega t_{0}\right)=0 \] The Melnikov condition for existence of chaotic behavior in the sense of Smale in the system (2) for sufficiently small \(\epsilon\) is

\[ \left|\tfrac{a}{a_{1}}\right| \lt \pi \omega \operatorname{csch}(\pi \omega) . \]

Remark. This result matches the estimate obtained in (Christie et al. 1995).

For example for \(N=1, ~ \alpha=1, ~ \beta=2, ~ \gamma=-1, ~ \delta=-2\), \(\omega=0.78, a=0.11, a_{1}=0.1\) the Melnikov function \(M\left(t_{0}\right)=0\) has no roots (fig. 2).

For \(N=1, \quad \alpha=1, \quad \beta=2, \gamma=-1, \delta=-2, \quad \omega=0.68, \quad a=0.12\), \(a_{1}=0.258, M\left(t_{0}\right)=0\) has root \(t_{0} \approx 6.95\) (with multiplicity two) (fig. 3).

3.2. The case \(N=2, \alpha=1, \beta=2, \gamma=-1, \delta=-2\)

Proposition 2. If \(N=2, \alpha=1, \beta=2, \gamma=-1, \delta=-2\), then the roots of Melnikov function \(M\left(t_{0}\right)\) are given as solutions of the equation

\[ M\left(t_{0}\right)=a+a_{1} \pi \omega c s c h(\pi \omega) \sin \left(\omega t_{0}\right)+2 a_{2} \pi \omega c s c h(2 \pi \omega) \sin \left(2 \omega t_{0}\right)=0 . \] For example for \(N=2, \alpha=1, \beta=2, \gamma=-1, \delta=-2, \omega=0.38, a=0.09\), \(a_{1}=0.137 ; a_{2}=0.09\) the Melnikov function \(M\left(t_{0}\right)=0\) is depicted in fig. 4.

For \(N=2, \quad \alpha=1, \quad \beta=2, \quad \gamma=-1, \quad \delta=-2, \quad \omega=0.62, \quad a=0.061\), \(a_{1}=0.05, a_{2}=0.07\) the Melnikov function \(M\left(t_{n}\right)=0\) has no roots (fig. 5).

3.3. General case

The research conducted above gives us a reason to formulate the main result.

Proposition A. For some \(N\) and fixed \(\alpha=1, \beta=2, \gamma=-1, \delta=-2\), the roots of the Melnikov function \(M\left(t_{0}\right)\) are given as solutions of the equation

(5)\[ M\left(t_{0}\right)=a+\pi \omega \sum_{i=1}^{N} i a_{i} \tfrac{\sin \left(i \omega t_{0}\right)}{\sinh (i \pi \omega)}=0 . \]

For example for \(N=3, \alpha=1, \beta=2, \gamma=-1, \delta=-2, \omega=0.28, a=\) \(0.09, a_{1}=0.137 ; a_{2}=0.09, a_{3}=0.2\) the Melnikov function \(M\left(t_{0}\right)=0\) is depicted in fig. 6.

From the above examples, the reader can himself formulate the corresponding Melnikov criterion for the occurrence of chaos in the considered dynamic system for some \(N\).

4. Some simulations

We will look at some interesting simulations on model (2):

Example 1. For given \(N=1, \omega=0.1, a=0.1, a_{1}=0.05, \epsilon=0.03\), \(\alpha=1, \beta=2, \gamma=-1, \delta=-2\) the simulations on the system (2) for \(x_{0}=0.9\);

\(y_{0}=0.1\) are depicted on fig. 7.

Example 2. For given \(N=3, \omega=0.1, a=0.1, a_{1}=0.05, a_{2}=0.02\), \(a_{3}=0.03, \epsilon=0.03, \alpha=1, \beta=2, \gamma=-1, \delta=-2\) the simulations on the system (2) for \(x_{0}=0.9 ; y_{0}=0.1\) are depicted on fig. 8.

Example 3. For given \(N=5, \omega=0.1, a=0.1, a_{1}=0.05, a_{2}=0.1\), \(a_{3}=0.03, a_{4}=0.11, a_{5}=0.01, \epsilon=0.03, \alpha=1, \beta=2, \gamma=-1, \delta=-2\) the simulations on the system (2) for \(x_{0}=0.8 ; y_{0}=0.2\) are depicted on fig. 9.

Example 4. For given \(N=7, \omega=0.95, a=0.05, a_{1}=0.5, a_{2}=0.01\), \(a_{3}=0.1, a_{4}=0.4, a_{5}=0.2, a_{6}=0.05, a_{7}=0.3, \epsilon=0.07, \alpha=1, \beta=\) \(2, \gamma=-1, \delta=-2\) the simulations on the system (2) for \(x_{0}=0.8 ; y_{0}=0.2\) are depicted on fig. 10.

With Examples 2-4, we demonstrate simulations of generated chaos in the new model.

5. A possible application of Melnikov functions in modelling and synthesis of radiating antenna patterns

Let us now focus on \(M(t)\). In our previous publications, we discussed a possible application of the Melnikov functions in modeling and synthesis of radiation antenna diagrams. For more details see for example (Kyurkchiev \& Andreev 2014), (Apostolov et al. 2018), (Kyurkchiev 2020), (Kyurkchiev et al. 2024).

We define the hypothetical normalized antenna factor as follows: \(M^{*}(\theta)=\) \(\tfrac{1}{D}\left|M\left(K \cos \theta+k_{1}\right)\right|\) where \(\theta\) is the azimuth angle; \(K=k d ; k=\tfrac{2 \pi}{\lambda} ; \lambda\) is the wave length; \(d\) is the distance between emitters; \(k_{1}\) is the phase difference.

Example 5. For fixed \(N=6, k=8, k_{1}=0, \omega-0.145, a=0.004\), \(a_{1}=0.137, a_{2}=0.09, a_{3}=0.2, a_{4}=0.3, a_{5}=0.3, a_{6}=0.2\) the Melnikov function and Melnikov antenna factor are depicted on fig. 11.

Example 6. For fixed \(N=10, k=8, k_{1}=0, \omega-0.145, a=0.001\), \(a_{1}=0.01, a_{2}=0.002, a_{3}=0.01, a_{4}=0.01, a_{5}=0.01, a_{6}=0.01, a_{7}=0.01\), \(a_{8}=0.01, a_{9}=0.01, a_{10}=0.2\) the Melnikov function and Melnikov antenna factor are depicted on fig. 12.

Example 7. For fixed \(N=12, k=8, k_{1}=0, \omega-0.145, a=0.0001\), \(a_{1}=0.01, a_{2}=0.002, a_{3}=0.01, a_{4}=0.01, a_{5}=0.01, a_{6}=0.01, a_{7}=0.01\), \(a_{8}=0.01, a_{9}=0.01, a_{10}=0.2, a_{11}=0.1, a_{12}-0.3\) the Melnikov function and Melnikov antenna factor are depicted on fig. 13.

Of course, this relatively new idea of justification and right to exist is subject to serious research by specialists working in this scientific direction.

The issue related to noise minimization (in decibels) also remains open.

6. Concluding Remarks

1. The discussions in this article are valid subject to the following limitations: \(0 \lt \tfrac{\alpha}{\beta} \lt 1,0 \lt \tfrac{\gamma}{\delta} \lt 1, \alpha \gamma \lt 0\), and \(\omega \gt 0\).

The existence of subharmonic periodic solutions of system (2) can be proved in the manner detailed in paper (Christie et al. 1995).

The system (1) \((\epsilon=0)\) has the Hamiltonian \(H(x, y)=x y(1-x)(1-y)=h\). It is known (Christie et al. 1995) that the family of periodic orbit surrounding the centre \(\left(\tfrac{1}{2}, \tfrac{1}{2}\right)\) with assumptions \(\alpha=1, \beta=2, \gamma=-1, \delta=-2\) is of the form

\[ q_{k}(t)=\left(x_{k}(t), y_{k}(t)\right)=\left(\tfrac{1}{2}-\tfrac{k}{2} \operatorname{sn}\left(\tfrac{t}{2}, k\right), \tfrac{1}{2}+\tfrac{k}{2} \tfrac{c n\left(\tfrac{t}{2}, k\right)}{\operatorname{dn}\left(\tfrac{t}{2}, k\right)}\right) \]

where \(k=\sqrt{1-16 h}\), and \(s n(u, k), c n(u, k), d n(u, k)\) are Jacobi elliptic functions with modulus \(k\).

The orbit \(q_{k}\) has period \(T(k)=8 K(k), k \in(0,1)\) and \(K(k)\) is the complete elliptic integral of the first kind.

For relatively positive integers \(m\) and \(n\) the subharmonic Melnikov integral is defined as

\[ M^{m / n}\left(t_{0}\right)=\int_{-\tfrac{m T}{2}}^{\tfrac{m T}{2}} f\left(q_{k}(t)\right) \wedge g\left(q_{k}(t), t+t_{0}\right) d t \]

and resonance condition is \(T(k)=m T\).

For some details see (Christie et al. 1995).

We note that for our new model (2) the first subharmonic Melnikov function is:

(6)\[ \begin{aligned} & M^{m / 1}\left(t_{0}\right)= \\ & =\int_{-\tfrac{m T}{2}}^{\tfrac{m T}{2}}\left(x_{k}(t)\left(1-x_{k}(t)\right)\left(1-2 y_{k}(t)\right)\left(a y_{k}(t)+\sum_{j=1}^{N} a_{j} \sin \left(j \omega\left(t+t_{0}\right)\right)\right)\right. \\ & \left.-y_{k}(t)\left(1-y_{k}(t)\right)\left(-1+2 x_{k}(t)\right)\left(a x_{k}(t)+\sum_{j=1}^{N} a_{j} \sin \left(j \omega\left(t+t_{0}\right)\right)\right)\right) d t \end{aligned} \]

We leave the calculation of the integrals of (6) to the reader.

Similar theoretical studies and simulations can be done on a wide class of reaction-kinetic models described in the literature.

The proposed new extended model contains many free parameters (the coefficients \(a_{i} ; i=1,2, \ldots, N\) ), which makes it attractive for use in the fields of biological applications (Murray 2001), chemistry, sociology, lifetime analysis, reaction kinetics (Michaelis & Menten 1913), biostatistics, population dynamic, medical research etc.

2. For example consider the following new predator-prey dynamical system:

(7)\[ \left\{\begin{array}{l} \tfrac{d x}{d t}=x(f-b y)+\epsilon\left(a x+\sum_{j=1}^{N} a_{j} \sin (j \omega t)\right) \\ \tfrac{d y}{d t}=y(d x-c)+\epsilon\left(a y+\sum_{j=1}^{N} a_{j} \sin (j \omega t)\right) \end{array}\right. \]

where \(0 \leq \epsilon \lt 1, a_{i} \geq 0 ; i=1,2, \ldots, N\) a and \(N\) is integer.

The first integral \((\epsilon=0): F(x, y)=d x+b y-c \ln x-f \ln y\) is constant on solution curves \((x(t), y(t))\).

Using the function \(F(x, y)\), one can furthermore show that all solution curves in the positive quadrant are closed, that is, all such solutions are periodic. For some details see (Prelle & Singer 1983).

The new class of extended models of the type (7) occurs in modeling of “predator-prey” systems can be studied in detail using the structural approach applied in this article for the model (2). For some models see (Kyurkchiev & Boyadjiev 2021), (Kyurkchiev et al. 2022), (Kyurkchiev et al. 2022a).

We will look at the following simulations on model (7)

Example 8. For given \(N=3, \omega=0.1, a=0.1, a_{1}=0.05, a_{2}=0.01\), \(a_{3}=0.02, \epsilon=0.03, f=b=1, c=d=-1\) the simulations on the system (7) for \(x_{0}=0.4 ; y_{0}=0.2\) are depicted on fig. 14.

Example 9. For given \(N=5, \omega=0.1, a=0.9, a_{1}=0.2, a_{2}=0.6\), \(a_{3}=0.8, a_{4}=0.7, a_{5}=0.9, \epsilon=0.03, f=b=1, c=d=-1\) the simulations on the system (7) for \(x_{0}=0.4 ; y_{0}=0.2\) are depicted on fig. 15.

3. We will note that the dynamics of the genetic model for the "battle of the sexes" by (Smith & Hofbauer 1987) can be studied with the methodology proposed in this article.

The derived results can be used as an integral part of a much more general application for scientific computing – for some details see (Golev et al. 2024), (Kyurkchiev et al. 2024), (Kyurkchiev et al. 2024a), (Kyurkchiev et al. 2024b), (Kyurkchiev & Zaevski 2023), (Kyurkchiev et al. 2023), (Kyurkchiev & Iliev 2022), (Kyurkchiev & Andreev 2014), (Apostolov et al. 2018), (Kyurkchiev 2020), (Kyurkchiev et al. 2024c), (Kyurkchiev et al. 2024d), (Kyurkchiev et al. 2024e), (Vasileva et al. 2024).

4. Nonstandard numerical methods connected to the investigation of the roots of equation \(M\left(t_{0}\right)=0\) can be found in ( (Proinov \& Vasileva 2020), (Ivanov 2024).

Acknowledgment

The first, second, forth, fifth and sixth authors are supported by the European Union-NextGenerationEU, through the National Plan for Recovery and Resilience of the Republic Bulgaria, project No BG-RP-2.004-0001-C01.

REFERENCES

APOSTOLOV, P., STEFANOV, A., APOSTOLOV, S., 2018. A study of filters selectivity with maximally flat responses with respect to Hausdorff distance. In Proceedings of the 2018 IX National Conference with International Participation (Electronica), Sofia, Bulgaria, 17-18 May 2018.

CHRISTIE, J.R., GOPALSAMY, K., LI, J., 1995. Chaos in sociobiology. Bulletin of the Australian Mathematical Society, vol. 51, pp. 439 – 451.

DAWKINS, R., 1976. The selfish gene. Oxford: Oxford University Press.

FOSTER, D., YOUNG, P., 1990. Stochastic evolutionary game dynamics. Theoretical Population Biology, vol. 38, pp. 219 – 232.

GOLEV, A., TERZIEVA, T., ILIEV, A., RAHNEV, A., KYURKCHIEV, N., 2024. Simulation on a generalized oscillator model: Web-based application. Comptes rendus de l’Academie bulgare des Sciences, vol. 77, pp. 230 – 237.

HOFBAUER, J. SIGMUND, K., 1988. The theory of evolution and \(d y\)namical systems. Cambridge: Cambridge University Press.

IVANOV, S., 2024. Families of high–order simultaneous methods with several corrections. Numerical Algorithms, vol. 97, no. 2, pp. \(945-\) 958.

KYURKCHIEV, N., 1963. Some Intrinsic Properties of Tadmor-Tanner Functions: Related Problems and Possible Applications. Mathematics, vol. 8, no. 11, p. 1963.

KYURKCHIEV, N., ANDREEV, A., 2014. Approximation and Antenna and Filters Synthesis. Some Moduli in Programming Environment MATHEMATICA. Saarbrucken: LAP LAMBERT Academic Publishing.

KYURKCHIEV, N., BOYADJIEV, G., 2021. Dynamics of Modified LotkaVolterra Model with Polynomial Intervention Factors. Methodological Aspects. III. International Journal of Differential Equations and Applications, vol. 20, pp. 121 – 132.

KYURKCHIEV, V., BOYADJIEV, G., Kyurkchiev, N., 2022. A Software Tool for Simulating the Dynamics of a New Extended Family of Lotka-Volterra Competition Model. International Journal of Differential Equations and Applications, vol. 21, pp. 33 – 46.

KYURKCHIEV, N., BOYADJIEV, G., KYURKCHIEV, V., MALINOVA, A. 2022a. A Technique for Simulating the Dynamics of Some Extended Nonlinear Models. International Electronic Journal of Pure and Applied Mathematics, vol. 16, pp. 13 – 25.

KYURKCHIEV, N., ILIEV, A., 2022. On the hypothetical oscillator model with second kind Chebyshev’s polynomial–correction: number and type of limit cycles, simulations and possible applications. Algorithms, vol. 15, p. 462.

KYURKCHIEV, V., ILIEV, A., RAHNEV, A., KYURKCHIEV, N., 2023. On a class of orthogonal polynomials as corrections in Lienard differential system. Applications. Algorithms, vol. 16, p. 297.

KYURKCHIEV, N., ZAEVSKI, T., 2023. On a hypothetical oscillator: investigations in the light of Melnikov’s approach, some simulations. International Journal of Differential Equations and Applications, vol. 22, pp. 67 – 79.

KYURKCHIEV, N., ZAEVSKI, T., ILIEV, A., KYURKCHIEV, V., RAHNEV, A., 2024. Nonlinear dynamics of a new class of micro-electromechanical oscillators – open problems. Symmetry vol. 16, p. 253.

KYURKCHIEV, N., ZAEVSKI, T., ILIEV, A., KYURKCHIEV, V., RAHNEV, A., 2024a. Modeling of Some Classes of Extended Oscillators: Simulations, Algorithms, Generating Chaos, Open Problems. Algorithms, vol. 17, p. 121.

KYURKCHIEV, N., ZAEVSKI, T., ILIEV, A., KYURKCHIEV, V., RAHNEV, A., 2024b. Generating Chaos in Dynamical Systems: Applications, Symmetry Results, and Stimulating Examples. Symmetry, vol. 16, p. 938.

KYURKCHIEV, N., ZAEVSKI, T., ILIEV, A., KYURKCHIEV, V., RAHNEV, A., 2024c. Dynamics of a new class of extended escape oscillators: Melnikov’s approach, possible application to antenna array theory. Mathematics and Informatics, vol. 67, no. 4, pp. 367 – 382.

KYURKCHIEV, N., ZAEVSKI, T., ILIEV, A., KYURKCHIEV, V., RAHNEV, A., 2024d. Notes on Modified Planar Kelvin-Stuart Models: Simulations, Applications, Probabilistic Control on the Perturbations. Axioms, vol. 13, p. 720.

KYURKCHIEV, N., ZAEVSKI, T., ILIEV, A., KYURKCHIEV, V., RAHNEV, A., 2024e. Dynamics of Some Perturbed Morse-Type Oscillators: Simulations and Applications. Mathematics, vol. 12, no. 21, p. 3368.

MELNIKOV, V. K., 1963. On the stability of a center for time–periodic perturbation. Transactions of the Moscow Mathematical Society, vol.

12, pp. 1 – 57.

MICHAELIS, L., MENTEN, M.L., 1913. Die Kinetik der Invertinwirkung. Biochemische Zeitschrift, vol. 49, pp. 333 – 369.

MURRAY, J. D., 2001. Mathematical biology. I. An introduction. 3rd ed. Berlin, Heidelberg: Springer-Verlag.

PRELLE, M., SINGER, M., 1983. Elementary first integrals of differential equations. Trans. of the Amer. Math. Soc. vol. 279, pp. 215 – 229.

PROINOV, P., VASILEVA, M., 2020. Local and Semilocal Convergence of Nourein’s Iterative Method for Finding All Zeros of a Polynomial Simultaneously. Symmetry, vol. 12, p. 1801.

SCHUSTER, P., SIGMUND, K., 1981. Coyness, philandering and stable strategies. Animal Behaviour, vol. 29, pp. 186 – 192.

SMITH, J.M., Hofbauer, J., 1987. The “battle of the sexes”: a genetic model with limit cycle behavior. Theoretical Population Biology, vol. 32, pp. 1 – 14.

TAYLOR, P.D., JONKER, L.B., 1978. Evolutionary stable strategies and game dynamics. Mathematical Biosciences, vol. 40, pp. 145 – 156.

VASILEVA, M., KYURKCHIEV, V., ILIEV, A., RAHNEV, A., ZAEVSKI, T., KYURKCHIEV, N., 2024. Some investigations and simulations on the generalized Rayleigh systems, Duffing systems with periodic parametric excitation, Mathieu and Hopf oscillators. Plovdiv: Plovdiv University Press.