Chaos Control and Synchronization in Nonlinear Prey-Predator Model of Type 1 Diabetes: A Modern IT Perspective

Maysoon M. Aziz¹, Qusay W. Habash²

¹Department of Mathematics, College of Computer Sciences and Mathematics, University of Mosul, Mosul, Iraq.
²Department of Mathematics, College of Computer Sciences and Mathematics, University of Tikrit, Tikrit, Iraq.

DOI: 10.4236/oalib.1113897                

Abstract

This paper investigates the behavior of a prey-predator system in continuous time obtained by non-linear differential equations arising from diabetes of type one. We study the stability of equilibrium points, dissipativity multistability, Wave Form, state space analysis Lyapunove exponents, and bifurcation. A new Lyapunove function is constructed, and the analysis of stability is consistent with other methods of stability. The analysis shows that the predator-prey system is unstable and chaotic with Kaplan-York dimension D_ky = 1.5121. A novel feature of the system has coexisting attractors and multistability for two different sets of initial conditions. Finally, the adaptive control Strategy based on Lyapunove’s method has been applied to investigate chaotic control and synchronization. Numerical simulations demonstrate that the proposed control laws successfully achieve master-Slave (drive-response) Synchronization and effective chaos Suppression. The analysis was conducted using modern programming languages, most notably MATLAB 2024 and Spss 25.

Keywords

Bifurcation, Stability, Multistability, Lyapunove Dimension, Synchronization, Control

1. Introduction

In 1926, the Italian mathematician Vito voltera first proposed simple differential equations to describe the population dynamics of two interacting species, a predator and its prey [1]. The birth of modern ecology recognized complex interactions between species and was equipped with mathematical modelling [2]. Many authors have studied the dynamics of the prey-predator differential equation model [3]-[6]. The behavior of the prey-predator model with holling type I and chaotic dynamics of discrete prey-predator with holling type II was studied [7] [8]. ElSadney [9] discussed dynamical complexity in the food chain. Suarez et al. [10] studied a prey-predator model, and age structure increases the instability of the coexistence equilibrium and leads to a coexistence attractor.

Diabetes is a disease, which can affect any human. This disease is usually of two types: the first is called type I which mainly hurts children and young people and is usually treated by insulin injection, the second type is called type II which hurts older people and is treated by drugs [10]-[12].

The paper is organized as follows:

The time evolution of a continuous-time two-dimensional prey-predator is described. Blood glucose level and the insulin doses were documented for each reading in a sample of 303 females and 178 males aged between 1 and 40 years. All individuals were diagnosed with type one diabetes and registered at Al-Wafa center of Diabetes and Endocrinology Research in Mosul, Iraq.

The existence of equilibrium points and their stability are analyzed in Section 3. In Section 4. Dissipativity of chaotic behavior is investigated by numerical simulation. Section 5 discusses lyapunove exponent and fractal dimension. Chaos treatment control and synchronization are in Section 6, Statistical analysis is in section 7. Finally, the conclusion of this paper is in section 8.

2. Model Formulation [13] [14]

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=ax\left(1-\frac{x}{k}\right)-bxy\\ \frac{\text{d}y}{\text{d}t}=-cy+dxy\end{array}$ (1)

X(t) represents a measure of blood sugar.

Y(t) is the insuline dose.

a is birth rate, c is the death rate, b is the conversion rate, e is the growth rate,

$a=\frac{\text{number of briths}}{\text{population mid year}}\times 1000$

$c=\frac{\text{number of total death}}{\text{Population mid year}}\times 1000$

$e=\frac{a-c}{\text{population mid year}}\times 1000$

$d=be$ , k is the carrying capacity

$k=\frac{e}{d}$ ,

From the set of raw data of 581 patients including (sex, age, blood sugar, insuline dose) from Alwafa-center, we get

$\begin{array}{l}a=35.5157,c=4.1869,e=31.3128\\ d=921.63,b=29.4331,k=0.8816\end{array}$ (2)

3. Equilibrium Points and It’s Stability

3.1. Equilibrium Points

The equilibrium points of system (1) particularly, we are interested in non negative or interior equilibrium points. All possible equilibrium points:

$E_0=\left(0,0\right)$ , is the trivial equilibrium point.

$E_1=\left(k,0\right)$ , is the axial equilibrium point, In the absence of predator (y=0).

$E_2=\left(\frac{c}{d},\frac{a}{b}\left(1-\frac{c}{kd}\right)\right)$ , is the interior equilibrium point.

$E_2$ to be positive, we need $\frac{c}{kd}>1$ , $a,b,c,d>0$ .

From the set of data, we get $E_0=\left(0,0\right)$ , $E_1=\left(0.8816,0\right)$ , $E_2=\left(0.00454,1.2004\right)$ .

3.2. Stability of System (1)

3.2.1. Local Stability

The local stability of system (1) is analyzed by computing the Jacobian at each equilibrium. The variation matrix of system (1) at the State variable is expressed as.

$J\left(x,y\right)=\left[\begin{array}{cc}a\left(1-\frac{2x}{k}-by\right)& -bx\\ dy& -c+dx\end{array}\right]$ (3)

Proposition 1:

The equilibrium point. E₀ locally asymptotically. unstable, if a, c > 1.

Otherwise stable fixed point.

Proof: To prove this result at E₀.

$J\left(x,y\right)=\left[\begin{array}{cc}a& 0\\ 0& -c\end{array}\right]$ , it’s characteristic equation

${\lambda }^2-\left(a-c\right)\lambda -ac=0$ (4)

From (4), the eigen values of. $J\left(x,y\right)$ .

${\lambda }_1=a$ , ${\lambda }_2=c$ , so from (2) we get E₀ unstable equilibrium point.

Proposition 2:

The equilibrium point E₁, locally unstable for a > c, otherwise it is a stable equilibrium point.

Proof:

The Jacobian matrix at E₁, given by

$J=\left[\begin{array}{cc}-a& -bk\\ 0& -c+dk\end{array}\right]$

The eigenvalues of the matrix are:

${\lambda }_1=-a$ and ${\lambda }_2=-c+kd$ . From (2) we get ${\lambda }_1=-35.5157$ , ${\lambda }_2=808.3221$ as a result of one eigen values being positive, System (1) is unstable. The local Stability of Interior equilibrium point, E₂ the Jacobian matrix (3) at E₂ has the form

$J=kd\left[\begin{array}{cc}-\frac{ac}{kd}& -\frac{bc}{d}\\ \frac{a}{b}\left(d-\frac{c}{k}\right)& 0\end{array}\right]$ (5)

It’s characteristic equation

${\lambda }^2-\text{Tr}\lambda +\text{Deter}J=0$ (6)

Where Tr is the trace and Deters is the determinant of the Jacobian matrix defined in Equation (5).

$\text{Tr}=\frac{-ac}{kd}$ , $\text{Deter}=ac-\frac{a{c}^2}{kd}$

$\text{Tr}=-0.183$ , $\text{Deter}=147.9343$

So, for the interior equilibrium point E₂ the roots of equation (6) are complex with negative real part

${\lambda }_{1,2}=-0.0915\pm 12.1625i$

Hence the interior equilibrium point E₂ is stable.

3.2.2. Routh Criteria [15] [16]

The system is asymptotically stable according to Routh stability test if all elements in the first column of Routh array are strictly positive, so from equations (4) and (2), we get

${\lambda }^2-31.3288\lambda -148.7006=0$ (7)

Since $a_0=-148.7006$ , $a_1=-31.3288$ , $a_2=1$ . The system (1) does not satisfy asymptotic Stability requirements for the first Column. Table 1 has two negative elements, therefore system (1) is unstable at E₀.

Table 1. Routh array for system (1) at E₀.

λ2	1	-148.7006
λ1	−31.3288	0
λ0	−148.7006	0

Using Routh array for the rest of Equilibrium E₁, and E₂, the result is reported in Table 2.

3.2.3. Lyapunove Function [17] [18]

The Lyapunove function of System (1)

If $v\left(x,y\right)=\frac{1}{2}\left(x^2+y^2\right)$

$\dot{v}\left(x,y\right)=\frac{\partial v}{\partial x}\frac{\text{d}x}{\text{d}t}+\frac{\partial v}{\partial y}\frac{\text{d}y}{\text{d}t}$ (8)

$\dot{v}\left(x,y\right)<0$ , exists, then system (1) is globally asymptotically stable.

By substituting system (1) into equation (8), we get

$\dot{v}\left(x,y\right)=35.5154x^2+40.4552x^3-29.4331x^{2}y-4.1869y^2+921.63x{y}^2$ (9)

At all equilibrium points using equation (9) we get system (1) unstable except at E₂, the result reported in Table 2.

3.2.4. Hwritz Criteria [19]

This criterion is employed to assess a system’s stability. If the square matrix A’s principal minors are all positive, the System is stable, if not, it is unstable from equation (4) which is

${\lambda }^2=31.32881\lambda -148.7006=0$ , we find

${\Delta }_1=a_{n-1}=a_1=-31.3288<0$

${\Delta }_2=\left|\begin{array}{cc}{a}_{n-1}& a_{n-3}\\ a_n& a_{n-2}\end{array}\right|=\left|\begin{array}{cc}{a}_1& 0\\ a_2& a_0\end{array}\right|=a_1{a}_0>0$

System is unstable for the rest of Equilibriums E₁, and E₂, the result reported in Table 2.

3.2.5. Continued Fraction Stability [20]

The criteria are applied to the characteristic equation of the continuous system (1), the equation

${\lambda }^2-31.3288\lambda -148.7006=0$

$Q_1\left(\lambda \right)={\lambda }^2-148.7006$

$Q_2\left(\lambda \right)=-31.32887\lambda$

$\frac{Q_1\left(\lambda \right)}{Q_2\left(\lambda \right)}=\frac{{\lambda }^2-148.7006}{-31.3288\lambda }=\frac{-1}{31.3288}\lambda +\frac{1}{0.2106\lambda }=h_1\lambda +\frac{1}{h_2\lambda }$

If all values of (h₁, h₂) are positive, then the roots of $Q\left(\lambda \right)=0$ , have negative real part (the System is stable), since

$h_1=-\frac{1}{31.3288}<0$ and $h_2=0.2106>0$

∴ system is unstable.

For the rest of Equilibriums E₁, E₂ the result is reported in Table 2.

3.2.6. A New Method for Lyapunove Function via Continued Fraction [21]

Lyapunov’s direct method is a cornerstone in Stability analysis; constructing an appropriate Lyapunove function remains good for nonlinear systems. In this study, we proposed a Novel approach to construct Lyapunove function using continued fraction criteria as a guiding principle.

Step 1: Using the continued fraction criterion, determine the characteristic polynomial factors.

$Q\left(\lambda \right)=a_0{\lambda }^n+a_1{\lambda }^{n-1}+a_2{\lambda }^{n-2}+\cdots +a_{n-1}\lambda +a_n$

Step 2: Construct the Lyapunove function as:

$V\left(X\right)=\frac{1}{2}{\displaystyle \sum _{j=0}^n\left|h_j\right|x_j^2}$

As we can see from paragraph (3.2.6)

$h_1=-\frac{1}{31.328}$ , $h_2=0.2106$

We assume the Lyapunove function as

$V\left(x,y\right)=\frac{1}{2}{\displaystyle \sum _{j=0}^n\left|h_j\right|x_j^2}=\frac{1}{2}\left(\left|h_1\right|x^2+\left|h_2\right|y^2\right)$

$\dot{V}\left(x,y\right)=\left|h_1\right|x\dot{x}+\left|h_2\right|y\dot{y}$

$\therefore \dot{V}\left(x,y\right)=\left|h_1\right|a{x}^2-\frac{\left|h_1\right|a{x}^3}{k}-\left|h_1\right|b{x}^{2}y-\left|h_2\right|c{y}^2+d\left|h_2\right|x{y}^2$ (9)

at E₀ using (9) we get.

$\therefore \dot{v}\left(x,y\right)=0$ , so the system (1) is unstable for the rest of Equilibrium E₁, E₂, the result reported in Table 2.

Table 2. Stability analysis.

Equilibrum Points	characteristic Equation Roots	Routh Stability	Hwruitz Stability	LyapunoveFunction	continued Fraction Stability	New Method Continued for Lyapunov
E0 (0, 0)	${\lambda }^2-31.328\lambda -148.700$	$a_0=148.7006$	${\Delta }_1=-31.32888$	0	$h_1=\frac{-1}{31.3288}\lambda$	0
	${\lambda }_1=-4.1869$	$a_1=-31.3288$	${\Delta }_2=4658.61135$		$h_2=0.2106\lambda$
	${\lambda }_2=35.5157$	$a_2=1$
		$b_1=-148.7006$
	(unstable)	(unstable)	(unstable)	(unstable)	(unstable)	(unstable)
E1 (0.8779, 0)	${\lambda }^2-773.7328\lambda -25096.5122$	$a_0=-25096.5122$	${\Delta }_1=-773.732$	−0.03112	$h_1=\frac{-1}{773.7328}\lambda$	7.942405
	${\lambda }_1=-311792$	$a_1=-773.7328$	${\Delta }_2=19417974.577$		$h_2=0.03083\lambda$
	${\lambda }_2=804.912$	$a_2=1$
		$b_1=-25096.5122$
	(unstable)	(unstable)	(unstable)	(unstable)	(unstable)	(unstable)
E2 (0.00454, 1.2006)	${\lambda }^2+0.1917\lambda -147.851202$	$a_0=147.79332$	${\Delta }_1=0.1917$	−0.0038917	$h_1=\frac{1}{0.1917}\lambda$	−2.52295
	${\lambda }_1=-0.09585-12.1566i$	$a_1=0.1917$	${\Delta }_2=0.1917$		$h_2=0.00129\lambda$
	$\begin{array}{c}{\lambda }_2=-0.09585+12.1566i\end{array}$	$a_2=1$
		$b_1=147.7931$
	(unstable)	(stable)	(stable)	(stable)	(stable)	(stable)

4. Dissipativity [22]

Conservativity and dissipativity are two properties for understanding the behavior of dynamical systems, and how they can be used for various applications.

The dissipative nature of a dynamical system ensures that it does not lose energy over time, while the conservative nature ensures that its momentum is preserved. The divergence of system (1) can be obtained as:

The Jacobian matrix (J) of system (1), can be represented as:

$J=\left[\begin{array}{cc}\frac{\partial f_1}{\partial x}& \frac{\partial f_1}{\partial y}\\ \frac{\partial f_2}{\partial x}& \frac{\partial f_2}{\partial y}\end{array}\right]$ , where $f_1=f_1\left(x,y\right)$ , $f_2=f_2\left(x,y\right)$ and $F=\left|\begin{array}{c}{f}_1\left(x,y\right)\\ f_2\left(x,y\right)\end{array}\right|$

$\nabla \cdot F=\frac{\partial f_1}{\partial x}+\frac{\partial f_2}{\partial y}=\text{Tr}\left(J\right)$ (10)

$\nabla \cdot F=a-c+x\left(d-\frac{2a}{k}\right)-by$ (11)

The divergence of system (1) depends on the state variable and parameters, hence the system’s dissipativity is state-dependent.

So, we substitute the equilibrium points Values into the divergence expression to analyze the local dissipativity

at E₀, $\nabla \cdot F=a-c>0$ .

at E₁, $\nabla \cdot F=kd-a-c>0$ .

at E₂, $\nabla \cdot F=\left(a-c\right)+0.00454\left(d-\frac{2a}{k}\right)-1.2004b>0$ .

Hence system(I), is locally unbounded (Nonuniformly conservative and Nonuniformly dissipative).

Preposition 3. system (1), is locally unbonded at all equilibrium points E₀, E₁ and E₂.

5. Numerical Simulation

This section presents selected numerical simulation results and discusses their implication, we solved the System (1) with parameter values as in (2) and initial condition, X(0) = 0.02 and y(0) = 1 using the standard fourth order Rung-kuta technique (Rk4) in MATLAB 24 with Step size h = 0.001. We present numerical evidence illustrating the qualitative dynamical behavior of system (1), Bifurcation diagrams, Lyapunove Exponent, Waveform analysis and multi stability.

5.1. Bifurcation Diagram

The bifurcation analysis was conducted over the parameter range, the numerical analysis started with initial conditions (x(0), y(0)) = (0.02, 1), while the observation time with t₀ = 0 the time step t_step = 0.5 and finishing time is tend = 300 seconds.

The bifurcation diagram of system (1), as the parameter varies in three cases.

1. Fixing parameters c, d, b, k, and varying a with respect to x, a ∈ (30,38).

2. Fixing parameters, a, d, c, k and varying b with respect to x, b ∈ (26, 32).

3. Fixing parameters, a, b, d, k, with respect to x, and varying parameter c.

The resulting, bifurcation, diagrams of system (1) are presented in Figures 1-3.

Figure 1. Bifurcation diagram of x versus a.

Figure 2. Bifurcation diagram of x versus b.

Figure 3. Bifurcation diagram of x versus c.

Figure 4. Lyapunove exponent for system (1).

5.2. Lyapunove Exponent [23] [24]

According to nonlinear chaos theory Lyapunove, exponent quantifies the exponential divergence or convergence of nearby trajectories. In the system state space, for system (1), using the parameter value in equation (2) and the initial value. [0.02, 1], the computed Lyapunove exponent is

$L_1=0.7252$ $L_2=-1.4162$

The corresponding Lyapunove dimension $D_{ky}=1+\frac{L_1}{\left|L_2\right|}=1.5121$ , as illustrated.

In Figure 4. The evolution of the Lyapunove exponents confirms, that System (1) exhibits chaotic behavior.

5.3. Wave Form Analysis

One of the fundamental features of chaotic dynamical system is the non-periodic structure of the waveform form as seen in Figure 5 and Figure 6 shows the aperiodic waveforms of x(t) and y(t) in time domain.

5.4. Phase Space Analysis

In this paragraph, Figure 7 shows the chaotic attractor in (x-y) plane for (x, y) = (0.02, 1).

5.5. Multistability

Multistability in a dynamical system occurs when multiple (two) attractors exist simultaneously for the same parameter.

Figure 5. Time in second versus X.

Figure 6. Time in second versus y.

Figure 7. Phase plot chaotic attractor of system (1), in x-y plane.

Set with each attractor depending on distinct initial conditions given in Table 3 and Figure 8 shows the Coexistence of multiple attractors. The Coexistence of multiple attractors in multistable System produces complex behavior, with the trajectory following different paths.

Table 3. Coexistence with the same parameter and different Ic’s.

Initial conditions	color	Parameter	Figure
X0 = 0.02, y0 = 1	Red	a = 35.5157 b = 29.4331 c = 4.1869 d = 921.63 k = 0.8816	Figure 8.
X0 = 0.01, y0 = 1	blue

Figure 8. View of multi stability of two attractors that match Table 3.

6. Chaos Treatment

6.1. Chaos Control

Several methods can be used for obtaining chaos control in continuous-time dynamical system. One of these methods is adaptive control.

Adaptive Control

We give an adaptive control method to stabilize chaotic orbit of system (I), with the following controlled form

$\begin{array}{l}\dot{x}=35.5157x-40.2852x^2-29.4331xy+u_1\left(t\right)\\ \dot{y}=-cy+921.63xy+u_2\left(t\right)\end{array}$ (12)

Where (12) is the designed adaptive control law with unknowing parameter c, and $u_1\left(t\right)$ , $u_2\left(t\right)$ are feedback controllers greater than zero to be designed and (x, y) are State variables. In order to ensure (12) globally converges to zero, consider the adaptive control functions are:

$\begin{array}{l}{u}_1\left(t\right)=-35.5157x+40.2855x^2-29.4331xy-M_{1}x\\ u_2\left(t\right)=\widehat{c}y-921.63xy-M_{2}y\end{array}$ (13)

Where, $\widehat{c}$ is the estimate parameter of c and M_i, (i = 1,2,3) are positive constants, substitute (13) into (12), we obtain:

$\begin{array}{l}\dot{x}=-M_{1}x\\ \dot{y}=\left(-c+en\widehat{c}\right)y-M_{2}y\end{array}$ (14)

Let the parameter error estimation

$e_c=c-\widehat{c}$ (15)

Substitute (15) into (14) we obtain

$\begin{array}{l}\dot{x}=-M_{1}x\\ \dot{y}=e_{c}y-M_{2}y\end{array}$ (16)

The Lyapunove approach, is applied to derive the update Law (16), which is used to adjust the estimation of parameter c.

Let the Lyapunov. Function

$v\left(x,y,e_c\right)=\frac{1}{2}\left(x^2+y^2+e_c^2\right)$ (17)

Which is positive definite on R Differentiating the parameter estimation error (15), with respect to time, we get

${\dot{e}}_c=-\dot{\widehat{c}}$ (18)

Differentiate (17) and substitute (18) and (16) we get

$\dot{v}\left(x,y,e_c\right)=-M_1{x}^2-M_2{y}^2+e_c\left(-\dot{\widehat{c}}-y^2\right)$ (19)

In equation (19), we up date estimated parameter by:

$\dot{\widehat{c}}=y^2+M_3{e}_c$ (20)

Where M₃ is constant greater then zero Now substitute (20) into (19) we get:

$\dot{v}\left(x,y,e_c\right)=-M_1{x}^2-M_2{y}^2-M_3{e}_c^2$ (21)

Which is negative definite on R³.

The following result is obtained by applying Lyapunove stability.

Proposition (4).

The chaotic system using lyapunove stability for three equilibrium points E₀, E₁, E₂ ∈ R² by adaptive control Law (13) and the update estimated parameter law $\widehat{c}=y^2+M_3{e}_c$ , where M₁, M₂ and M₃ are positive constants.

6.2. Adaptive Synchronization

6.2.1. Theoretical Results

In order to achieve adaptive synchronization two systems are needed with unknown parameter c, the first system is called the drive (master) system corresponding to the uncontrolled system (1) and the second is called the response system (slave) depicted in (22) and (23), respectively.

$\begin{array}{l}{\dot{x}}_1=35.5157x_1\left(1-\frac{x_1}{0.8816}\right)-29.43311x_1{x}_2\\ {\dot{x}}_2=-c{x}_1+92.63x_1{x}_2\end{array}$ (22)

For the response system, the controlled dynamics are represented as

$\begin{array}{l}{\dot{y}}_1=35.5157y_1-40.2852y_1^2-29.4331y_1{y}_2+u_1\\ {\dot{y}}_2=-c{y}_2+921.63y_1{y}_2+u_2\end{array}$ (23)

Where $y_1,y_2$ are the state variables and $u_1,u_2$ are the nonlinear controllers to be designed. The synchronization error between (22) and (23) is given by.

$e_i=y_i-x_i,\left(i=1,2\right)$ (24)

Hence

$\begin{array}{l}{\dot{e}}_1=35.5157e_1-40.28522e_1^2-29.4331\left(e_1{e}_2+y{e}_1+x{e}_2\right)+u_1\\ {\dot{e}}_2=-c{e}_2+921-63\left(e_1{e}_2+y{e}_1+x{e}_2\right)+u_2\end{array}$ (25)

Now define the adaptive control function U₁(t) and U₂(t) as

$\begin{array}{l}{u}_1=-35.51579e_1+40.2852e_1^2+29.4331\left(e_1{e}_2+y{e}_1+x{e}_2\right)-M_1{e}_1\\ u_2=\widehat{c}{e}_2-921.63\left(e_1{e}_2+y{e}_1+x{e}_2\right)-M_2{e}_2\end{array}$ (26)

Where M₁, M₂ are positive real values and $\widehat{c}$ is the estimated value of the parameter c substitute (26) into (25), we get a dynamical System of the synchronization error as:

$\begin{array}{l}{\dot{e}}_1=-M_1{e}_1\\ {\dot{e}}_2=e_c{e}_2-M_2{e}_2-M_2{e}_2\end{array}$ (27)

Where $e_c=-\left(c+\widehat{c}\right)$ .

The lyapunove approach is used to prove the stability, the quadratic lyapunove function is considered as:

$V\left(e_1,e_2,e_c\right)=\frac{1}{2}\left(e_1^2+e_2^2+e_c^2\right)$ (28)

Which is positive definite on R³.

Note that

${\dot{e}}_c=-\widehat{c}$ (29)

Differentiating equation (28) and Substitute (27) and (29), we get

$\dot{v}=-M_1{e}_1^2-M_2{e}_2^2-e_c\left[e_1{e}_2-\dot{\widehat{c}}\right]$ (30)

The estimated parameter is updated by the following law

$\dot{\widehat{c}}=e_1{e}_2+M_3{e}_c$ (31)

Where M₃ is a positive constant Substitute (31) into (30), we get $\dot{v}=-M_1{e}_1^2-M_2{e}_2^2-M_3{e}_c^2$ , which is negative on R³.

Hence according to lyapunove Stability [24], the parameter estimation error and the synchronization error both decline to zero exponentially. Therefore, the following proposition (5) is established.

Proposition 5.

The identical chaotic systems-drive system (22) and response system (23) with unknown parameter c- achieve Synchronization of any initial condition through the adaptive control technique (26), where the parameter estimate is defined by (31) and the constant M_i (i = 1,2,3) are positive.

6.2.2. Numerical Results

Runge-Kutta method of the fourth order was used to solve the dynamic. System (22) and (23) as well as the synchronization of the error dynamics described in equation (27). The initial condition for the derive system (22) was set as (x₁(0), x₂(0)) = [0.02, 1] and for the response system (23) as (y₁(0), y₂(0)) = [2, 5] and assigned a value of c =4.1860.

The convergence of the synchronization error for system (27) is shown in Figure 9.

Figure 9. Trajectory convergence of the synchronization error dynamic (27).

7. Statistical Analysis

By using spss. 25 program for the set of data at Al-Wafa center of Diabetes and Endocrinology Research in Mosul, Iraq, we get the following results.

1. Table 4 and Table 5 display the research Sample consisting of 278 male and 303 females distributed a cording to age groups as illustrated in Figure 10 and Figure 11 shows approximate equality in number between males and females a cross all age groups except for the third group (21 - 30) years Where female make up to 37% compared to 30% for males.

2. The Statistical measure for the data given in Table 6 shows that males had an average blood sugar level of about 14 with a standard deviation (sd) of 68, while females had an average of 12 with s.d of 6.4.

Male received an average insulin dose of 44.8 with s.d of 21.2 while female received an average dose of 48.8 with s.d of 24.

Table 4. Total research sample according to age groups.

T
Age	Frequency	Ratio
1 - 10	54	9.3
11 - 20	166	28.6
21 - 30	196	33.7
31 - 40	165	28.4
Total	581	100.0

Table 5. Research sample of male and female distributed according to age groups.

Age	Frequency	Ratio
1 - 10	30	10.8
11 - 20	82	29.5
21 - 30	83	29.9
31 - 40	83	29.9
Total male	278	100.0
1 - 10	24	7.9
11 - 20	84	27.7
21 - 30	113	37.3
31 - 40	82	27.0
Total female	303	100.0

Figure10. The total research sample distributed according to age groups, as given in Table 4.

Figure 11. Show approximate equality in number between males and females a cross all age groups except the third group.

Table 6. Descriptive statistics for both blood sugar concentration and insulin dose.

		N	Mean	Std. Deviation	Minimum	Maximum
Blood sugar level	male	278	14.2119	6.80848	2.20	34.00
	female	303	11.8868	6.40992	2.40	38.10
Insulin dose	male	278	44.8022	21.23397	3.00	135.00
	female	303	48.8119	24.01974	3.00	208.00

3. The data did not confirm to the properties of Normal distribution by Kolmogorov Smirnov and Shapiro-wilk tests as shown in Table 7, most of the P-values were less than 1%.

Table 7. Test for normal distribution of data.

		Kolmogorov-Smirnov			Shapiro-Wilk
		Statistic	Df	Sig.	Statistic	df	Sig.
Blood sugar concentration	male	0.073	278	0.001	0.969	278	0.000
	female	0.124	303	0.000	0.929	303	0.000
Insulin dose	mail	0.070	278	0.002	0.979	278	0.000
	female	0.088	303	0.000	0.934	303	0.000
Blood sugar concentration	1 - 10	0.149	54	0.004	0.935	54	0.006
	11 - 20	0.121	166	0.000	0.949	166	0.000
	21 - 30	0.057	196	0.200	0.963	196	0.000
	31 - 40	0.118	165	0.000	0.926	165	0.000
Insulin dose	1 - 10	0.112	54	0.086	0.908	54	0.001
	11 - 20	0.091	166	0.002	0.971	166	0.001
	21 - 30	0.066	196	0.037	0.985	196	0.040
	31 - 40	0.133	165	0.000	0.900	165	0.000

4. A correlation and regression analysis of sugar between blood sugar levels (independent variable) and insulin doses (dependent variable) showed a very weak positive Pearson correlation r = 0.084, which was significant across most subgroups. Linear regression yields a low R² value overall 0.115, indicating a poor fit. The power model provided the best fit R² = 0.938 and less standard error of the estimates (SE = 0.619) with parameter estimation optimized via 168 iterations, we get the proper Regression model, which is the power Model (which is nonlinear) whose equation is

$\text{Insuline}=b_0\times {\left[\text{blood sugar}\right]}^{b_1}$

OR

$\ln \left(y\right)=\ln b_0+b_1\ln \left(x\right)$ .

The parameters b₀ and b₁ are estimated from a set of data using SPSS version 25, shown in Table 8.

Table 8. Parameter estimation of nonlinear regression parameter estimates.

Parameter	Estimate	Std. Error	95% Confidence Interval
			Lower Bound	Upper Bound
b0	6753317488501340.000	28545044677664700.000	−49311137106564600.000	62817772083567300.000
b1	4.311	1.701	0.970	7.652

8. Conclusion

In this paper, we explored the complex dynamics of a continuous-time prey-predator model derived from type 1 diabetes through nonlinear differential equations. Our investigation revealed rich dynamical behaviors, including instability, chaos and multistability analysis, using a newly constructed Lyapunove function confirmed consistency with classical Stability methods. The system demonstrated chaotic dynamics, with Lyapunove dimension D_ky = 1.5121 indicative of sensitive dependence on initial conditions. Notably, the existence of coexisting attractors under different initial conditions highlights the system’s multistable nature.

To address the chaotic behavior, we developed and effectively implemented an adaptive control strategy based on Lyapunor theory. Numerical simulation confirmed this control.

Approach effectively achieves synchronization in the master-slave (drive-response) configuration and suppresses chaos.

An additional statistical analysis of patient data using SPSS showed that blood glucose concentration and saline dosage do not follow a normal distribution. A power-law fit with R² = 0.938 indicates a potential nonlinear or fractal structure, but this alone is not enough to confirm chaotic behavior without additional dynamical analysis.

This research contributes a novel control framework for managing chaos in biology-inspired models and offers valuable insights into the nonlinear nature of diabetes-related physiological interactions.

Future work will extend the model to incorporate real-time clinical feedback and explore control robustness under uncertainty.

Conflicts of Interest

The authors declare no conflicts of interest.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Volterra, L. (1926) Dafnie pelagiche di due laghi dell’Italia Centrale. Internationale Revue der gesamten Hydrobiologie und Hydrographie, 15, 204-239. https://doi.org/10.1002/iroh.19260150303
[2]	Lotka, A.J. (1920) Analytical Note on Certain Rhythmic Relations in Organic Systems. Proceedings of the National Academy of Sciences of the United States of America, 6, 410-415. https://doi.org/10.1073/pnas.6.7.410
[3]	Hainzl, J. (1988) Stability and Hopf Bifurcation in a Predator-Prey System with Several Parameters. SIAM Journal on Applied Mathematics, 48, 170-190. https://doi.org/10.1137/0148008
[4]	He, X. (1996) Stability and Delays in a Predator-Prey System. Journal of Mathematical Analysis and Applications, 198, 355-370. https://doi.org/10.1006/jmaa.1996.0087
[5]	Murray, J.D. (2007) Mathematical Biology: I. An Introduction (Vol. 17). Springer Science & Business Media.
[6]	Aziz, M. and Hamid, M. (2019) The Possibility of Increasing the Predictability Indices after Control of 3D-Continuous-Time System. 2019 International Conference on Computing and Information Science and Technology and Their Applications (ICCISTA), Kirkuk, 3-5 March 2019, 1-5. https://doi.org/10.1109/iccista.2019.8830650
[7]	Aziz, M.M. and Al-Nuaimi, Z.A. (2015) Delay and Non-Delay Differential Models of Diabetes Type 1. International Journal of Electronics Communication and Computer Engineering, 6, 277-289.
[8]	Agiza, H.N., ELabbasy, E.M., EL-Metwally, H. and Elsadany, A.A. (2009) Chaotic Dynamics of a Discrete Prey-Predator Model with Holling Type II. Nonlinear Analysis: Real World Applications, 10, 116-129. https://doi.org/10.1016/j.nonrwa.2007.08.029
[9]	Elsadany, A.E.A., El-Metwally, H.A., Elabbasy, E.M. and Agiza, H.N. (2012) Chaos and Bifurcation of a Nonlinear Discrete Prey-Predator System. Computational Ecology and Software, 2, 169-180.
[10]	Suarez, L., Cameron, M.K., Fagan, W.F. and Levy, D. (2025) A Predator-Prey Model with Age-Structured Role Reversal. arXiv: 2502.19748.
[11]	Beniwal, N.S. and Beniwal, R. (2009) Automatic Control Systems with Matlab Programming. Laxmi Publications, Ltd.
[12]	Aziz, M.M. and Jihad, O.M. (2021) Stability, Chaos Tests with Adaptive and Feed-back Control Methods for 3D Discrete-Time Dynamical System. International Journal of Electronics Communication and Computer Engineering, 12, 31-42.
[13]	Jyothsna, N., Ramya, A., Abhilash, K. and Johnson, B.L. (2021) Pattern of Antibiotic Resistance of Various Strains of Bacteria Causing Acute Tonsillitis. International Journal of Otorhinolaryngology and Head and Neck Surgery, 7, 994-1003. https://doi.org/10.18203/issn.2454-5929.ijohns20212122
[14]	El Moutaouakil, K.E., El Ouissari, A.E., Palade, V., Charroud, A., Olaru, A., Baïzri, H., et al. (2023) Multi-Objective Optimization for Controlling the Dynamics of the Diabetic Population. Mathematics, 11, Article 2957. https://doi.org/10.3390/math11132957
[15]	Boccara, N. (2010) Spatial Models. In: Boccara, N., Ed., Modeling Complex Systems, Springer, 219-323. https://doi.org/10.1007/978-1-4419-6562-2_6
[16]	Pinheiro, R.F., Fonseca-Pinto, R. and Colón, D. (2024) A Review of the Lurie Problem and Its Applications in the Medical and Biological Fields. AIMS Mathematics, 9, 32962-32999. https://doi.org/10.3934/math.20241577
[17]	De Gaetano, A., Gaz, C. and Panunzi, S. (2019) Consistency of Compact and Extended Models of Glucose-Insulin Homeostasis: The Role of Variable Pancreatic Reserve. PLOS ONE, 14, e0211331. https://doi.org/10.1371/journal.pone.0211331
[18]	Seborg, D.E., Edgar, T.F., Mellichamp, D.A. and Doyle III, F.J. (2016) Process Dynamics and Control. John Wiley & Sons.
[19]	Runa and Sharma, R. (2015) A Lyapunove Theory Based Adaptive Fuzzy Learning Control for Robotic Manipulator. 2015 International Conference on Recent Developments in Control, Automation and Power Engineering (RDCAPE), Noida, 12-13 March 2015, 247-252. https://doi.org/10.1109/rdcape.2015.7281404
[20]	Aziz, M.M. and Kalalf, A.A. (2023) Nonlinear 6D Dynamical System with Hidden Attractors and Its Electronic Circuit. Open Access Library Journal, 10, 1-16. https://doi.org/10.4236/oalib.1109674
[21]	Maysoon M. Aziz, and Qusay W. Habash, (2025) Stability and Chaos Control in a Novel Three-Dimensional Multistable Dynamical System with Coexisting Attractors. International Journal of Computational and Experimental Science and Engineering, 11, 3573-3581. https://doi.org/10.22399/ijcesen.1635
[22]	Heidel, J. and Zhang, F. (2007) Nonchaotic and Chaotic Behavior in Three-Dimensional Quadratic Systems: Five-One Conservative Cases. International Journal of Bifurcation and Chaos, 17, 2049-2072. https://doi.org/10.1142/s021812740701821x
[23]	Hu, J., Qi, G., Wang, Z. and Chen, G. (2021) Rare Energy-Conservative Attractors on Global Invariant Hypersurfaces and Their Multistability. International Journal of Bifurcation and Chaos, 31, Article ID: 2130007. https://doi.org/10.1142/s021812742130007x
[24]	Aziz, M.M. and Jihad, O.M. (2021) Stability and Chaos Tests of 2D Discrete Time Dynamical System with Hidden Attractors. Open Access Library Journal, 8, 1-11. https://doi.org/10.4236/oalib.1107501